样条插值及其（精确）导数
spline interpolation and its (exact) derivatives
假设我有以下数据和命令：
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)
您可以看到它们工作得很好，从某种意义上说，插值函数可以很好地匹配节点处的数据点。但我的问题是，我需要计算 y 的精确导数（即 p 函数）w.r.t。时间并将其绘制在 t 向量上。如何做呢？我不会使用 diff 命令，因为我需要确保导数函数与 t 向量的长度相同。非常感谢。
方法A：使用导数
此方法计算多项式的实际导数。如果你有曲线拟合工具箱，你可以使用：
% 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);
如果您没有曲线拟合工具箱，您可以将 fnder 行替换为：
% 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);
来源：This mathworks question。 感谢 m7913d 提供链接。
附录：
注意
p = interp1(t,y,t_new,'spline');
是
的快捷方式
% 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);
为了得到导数，我们显然需要多项式，不能只使用新的插值点。为避免两次插入点可能需要很长时间才能处理大量数据，您应该将快捷方式替换为更长的版本。因此，包含您的代码示例的完整工作示例将是：
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);
方法 B：使用有限差分
连续函数的导数在其基础上只是 f(x) 与 f(x+无穷小差) 除以所述无穷小差。
在 matlab 中，eps 是双精度可能的最小差异。因此，在每个 t_new 之后，我们添加一个更大 eps 的第二个点，并为新点插入 y。然后每个点和它的 +eps 对之间的差除以 eps 得到导数。
问题是，如果我们处理如此小的差异，输出导数的精度将受到严重限制，这意味着它只能有整数值。因此，我们添加的值略大于 eps 以允许更高的精度。
% 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);
如果你想得到旧点的导数，你当然可以用t替换t_new（或任何其他你想得到微分的时间）。
在您的情况下，此方法略逊于方法 a)，因为它速度较慢且不太精确。但也许它对处于不同情况的其他人有用。
假设我有以下数据和命令：
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)
您可以看到它们工作得很好，从某种意义上说，插值函数可以很好地匹配节点处的数据点。但我的问题是，我需要计算 y 的精确导数（即 p 函数）w.r.t。时间并将其绘制在 t 向量上。如何做呢？我不会使用 diff 命令，因为我需要确保导数函数与 t 向量的长度相同。非常感谢。
方法A：使用导数
此方法计算多项式的实际导数。如果你有曲线拟合工具箱，你可以使用：
% 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);
如果您没有曲线拟合工具箱，您可以将 fnder 行替换为：
% 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);
来源：This mathworks question。 感谢 m7913d 提供链接。
附录：
注意
p = interp1(t,y,t_new,'spline');
是
的快捷方式% 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);
为了得到导数，我们显然需要多项式，不能只使用新的插值点。为避免两次插入点可能需要很长时间才能处理大量数据，您应该将快捷方式替换为更长的版本。因此，包含您的代码示例的完整工作示例将是：
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);
方法 B：使用有限差分
连续函数的导数在其基础上只是 f(x) 与 f(x+无穷小差) 除以所述无穷小差。
在 matlab 中，eps 是双精度可能的最小差异。因此，在每个 t_new 之后，我们添加一个更大 eps 的第二个点，并为新点插入 y。然后每个点和它的 +eps 对之间的差除以 eps 得到导数。
问题是，如果我们处理如此小的差异，输出导数的精度将受到严重限制，这意味着它只能有整数值。因此，我们添加的值略大于 eps 以允许更高的精度。
% 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);
如果你想得到旧点的导数，你当然可以用t替换t_new（或任何其他你想得到微分的时间）。
在您的情况下，此方法略逊于方法 a)，因为它速度较慢且不太精确。但也许它对处于不同情况的其他人有用。