Non-symbolic derivative at all sample points including boundary points

Question

Suppose I have a vector t = [0 0.1 0.9 1 1.4], and a vector x = [1 3 5 2 3]. How can I compute the derivative of x with respect to time that has the same length as the original vectors?

I should not use any symbolic operations. The command diff(x)./diff(t) does not produce a vector of the same length. Should I first interpolate the x(t) function and then take its derivative?

Answer 1

Different approaches exist to calculate the derivative at the same points as your initial data:

1. Finite differences: Use a central difference scheme at your inner points and a forward/backward scheme at your first/last point

or

2. Curve fitting: Fit a curve through your points, calculate the derivative of this fitted function and sample them at the same points as the original data. Typical fitting functions are polynomials or spline functions.

Note that the curve fitting approach gives better results, but needs more tuning options and is slower (~100x).

Demonstration

As an example, I will calculate the derivative of a sine function:

t = 0:0.1:1;
y = sin(t);

Its exact derivative is well known:

dy_dt_exact = cos(t);

The derivative can approximately been calculated as:

1. Finite differences:

   dy_dt_approx = zeros(size(y));
   dy_dt_approx(1) = (y(2) - y(1))/(t(2) - t(1)); % forward difference
   dy_dt_approx(end) = (y(end) - y(end-1))/(t(end) - t(end-1)); % backward difference
   dy_dt_approx(2:end-1) = (y(3:end) - y(1:end-2))./(t(3:end) - t(1:end-2)); % central difference
   

or

2. Polynomial fitting:

   p = polyfit(t,y,5); % fit fifth order polynomial
   dp = polyder(p); % calculate derivative of polynomial
   

The results can be visualised as follows:

figure('Name', 'Derivative')
hold on
plot(t, dy_dt_exact, 'DisplayName', 'eyact');
plot(t, dy_dt_approx, 'DisplayName', 'finite difference');
plot(t, polyval(dp, t), 'DisplayName', 'polynomial');
legend show

figure('Name', 'Error')
hold on
plot(t, abs(dy_dt_approx - dy_dt_exact)/max(dy_dt_exact), 'DisplayName', 'finite difference');
plot(t, abs(polyval(dp, t) - dy_dt_exact)/max(dy_dt_exact), 'DisplayName', 'polynomial');
legend show

The first graph shows the derivatives itself and the second graph plots the relative errors made by both methods.

Discussion

One clearly sees that the curve fitting method gives better results than the finite differences, but it is ~100x slower. The curve fitting methods has a relative error of order 10^-5. Note that the finite differences approach becomes better when your data is sampled more densely or you use a higher order scheme. The disadvantage of the curve fitting approach is that one has to choose a good polynomial order. Spline functions may be better suited in general.

A 10x faster sampled dataset, i.e. t = 0:0.01:1;, results in the following graphs: