Suppose I have the following data and commands:
clc;clear;
t = [0:0.1:1];
t_new = [0:0.01:1];
y = [1,2,1,3,2,2,4,5,6,1,0];
p = interp1(t,y,t_new,'spline');
plot(t,y,'o',t_new,p)
You can see they work quite fine, in the sense interpolating function matches the data points at the nodes fine. But my problem is, I need to compute the exact derivative of y (i.e., p function) w.r.t. time and plot it against the t vector. How can it be done? I shall not use diff commands, because I need to make sure the derivative function has the same length as t vector. Thanks a lot.
This method calculates the actual derivative of the polynomial. If you have the curve fitting toolbox you can use:
% calculate the polynominal
pp = interp1(t,y,'spline','pp')
% take the first order derivative of it
pp_der=fnder(pp,1);
% evaluate the derivative at points t (or any other points you wish)
slopes=ppval(pp_der,t);
If you don't have the curve fitting toolbox you can replace the fnderline with:
% piece-wise polynomial
[breaks,coefs,l,k,d] = unmkpp(pp);
% get its derivative
pp_der = mkpp(breaks,repmat(k-1:-1:1,d*l,1).*coefs(:,1:k-1),d);
Source: This mathworks question. Thanks to m7913d for linking it.
Appendix:
Note that
p = interp1(t,y,t_new,'spline');
is a shortcut for
% get the polynomial
pp = interp1(t,y,'spline','pp');
% get the height of the polynomial at query points t_new
p=ppval(pp,t_new);
To get the derivative we obviously need the polynomial and can't just work with the new interpolated points. To avoid interpolating the points twice which can take quite long for a lot of data, you should replace the shortcut with the longer version. So a fully working example that includes your code example would be:
t = [0:0.1:1];
t_new = [0:0.01:1];
y = [1,2,1,3,2,2,4,5,6,1,0];
% fit a polynomial
pp = interp1(t,y,'spline','pp');
% get the height of the polynomial at query points t_new
p=ppval(pp,t_new);
% plot the new interpolated curve
plot(t,y,'o',t_new,p)
% piece-wise polynomial
[breaks,coefs,l,k,d] = unmkpp(pp);
% get its derivative
pp_der = mkpp(breaks,repmat(k-1:-1:1,d*l,1).*coefs(:,1:k-1),d);
% evaluate the derivative at points t (or any other points you wish)
slopes=ppval(pp_der,t);
A derivative of a continuous function is at its base just the difference of f(x) to f(x+infinitesimal difference) divided by said infinitesimal difference.
In matlab, eps is the smallest difference possible with a double precision. Therefore after each t_new we add a second point which is eps larger and interpolate y for the new points. Then the difference between each point and it's +eps pair divided by eps gives the derivative.
The problem is that if we work with such small differences the precision of the output derivatives is severely limited, meaning it can only have integer values. Therefore we add values slightly larger than eps to allow for higher precisions.
% how many floating points the derivatives can have
precision = 10;
% add after each t_new a second point with +eps difference
t_eps=[t_new; t_new+eps*precision];
t_eps=t_eps(:).';
% interpolate with those points and get the differences between them
differences = diff(interp1(t,y,t_eps,'spline'));
% delete all differences wich are not between t_new and t_new + eps
differences(2:2:end)=[];
% get the derivatives of each point
slopes = differences./(eps*precision);
You can of course replace t_new with t (or any other time you want to get the differential of) if you want to get the derivatives at the old points.
This method is slightly inferior to method a) in your case, as it is slower and a bit less precise. But maybe it's useful to somebody else who is in a different situation.