4

I have a problem with numerical derivative of a vector that is x: Nx1 with respect to another vector t (time) that is the same size of x. I do the following (x is chosen to be sine function as an example):

t=t0:ts:tf;
x=sin(t);
xd=diff(x)/ts;

but the answer xd is (N-1)x1 and I figured out that it does not compute derivative corresponding to the first element of x.

is there any other way to compute this derivative?

0m3r
  • 12,286
  • 15
  • 35
  • 71
milad
  • 169
  • 1
  • 2
  • 10
  • 3
    Numerically that's the way to go. It's natural that you get one value less than the original, since you are using _differences_ to approximate the derivative – Luis Mendo Aug 11 '14 at 14:14
  • I know that. is there any other way to compute that? – milad Aug 11 '14 at 14:15
  • 2
    _Numerically_ the only approach I know of is based on differences. You could also compute _symbolically_, if you know the analytical expression of the function – Luis Mendo Aug 11 '14 at 14:23
  • the analytical function is not available and x is just a series of data points. – milad Aug 11 '14 at 14:30
  • @milad - I removed my answer as it wasn't quite useful. Hope you solve your question soon! – rayryeng Aug 11 '14 at 14:34
  • Related: [Calculate the derivative of a vector](http://stackoverflow.com/questions/20613859/calculate-the-derivative-of-a-vector). – horchler Aug 11 '14 at 15:20

3 Answers3

9

You are looking for the numerical gradient I assume.

t0 = 0;
ts = pi/10;
tf = 2*pi;

t  = t0:ts:tf;
x  = sin(t);
dx = gradient(x)/ts

The purpose of this function is a different one (vector fields), but it offers what diff doesn't: input and output vector of equal length.

gradient calculates the central difference between data points. For an array, matrix, or vector with N values in each row, the ith value is defined by

enter image description here

The gradient at the end points, where i=1 and i=N, is calculated with a single-sided difference between the endpoint value and the next adjacent value within the row. If two or more outputs are specified, gradient also calculates central differences along other dimensions. Unlike the diff function, gradient returns an array with the same number of elements as the input.

Robert Seifert
  • 25,078
  • 11
  • 68
  • 113
2

I know I'm a little late to the game here, but you can also get an approximation of the numerical derivative by taking the derivatives of the polynomial (cubic) splines that runs through your data:

function dy = splineDerivative(x,y)

% the spline has continuous first and second derivatives
pp = spline(x,y); % could also use pp = pchip(x,y);

[breaks,coefs,K,r,d] = unmkpp(pp);
% pre-allocate the coefficient vector
dCoeff = zeroes(K,r-1);

% Columns are ordered from highest to lowest power. Both spline and pchip
% return 4xn matrices, ordered from 3rd to zeroth power. (Thanks to the
% anonymous person who suggested this edit).
dCoeff(:, 1) = 3 * coefs(:, 1); % d(ax^3)/dx = 3ax^2;
dCoeff(:, 2) = 2 * coefs(:, 2); % d(ax^2)/dx = 2ax;
dCoeff(:, 3) = 1 * coefs(:, 3); % d(ax^1)/dx = a;

dpp = mkpp(breaks,dCoeff,d);

dy = ppval(dpp,x);

The spline polynomial is always guaranteed to have continuous first and second derivatives at each point. I haven not tested and compared this against using pchip instead of spline, but that might be another option as it too has continuous first derivatives (but not second derivatives) at every point.

The advantage of this is that there is no requirement that the step size be even.

craigim
  • 3,884
  • 1
  • 22
  • 42
0

There are some options to work-around your issue.

First: you can make your domain larger. Instead of N, use N+1 gridpoints.

Second: depending on the end-point of interest, you can use

  1. Forward difference: F(x + dx) - F(x)
  2. Backward difference: F(x) - F(x - dx)
gire
  • 1,105
  • 1
  • 6
  • 16
  • 1
    Dear gire I have the x values as a set of data and i dont have the analytic function. so how i can extend the domain? I think forward difference is the same as ``diff`` . – milad Aug 11 '14 at 14:28
  • @milad by extending the domain I mean: if your `t = [0, 1]` and you are interested in the derivative `dF(t = 0)` and `dF(t = 1)`, then you can extend `t` to `t [-0.5, 1.5]` (where `F` is your function and `dF` represents the numerical derivative) – gire Aug 11 '14 at 14:33