MATLAB 怎样绘制拟合出的Smoothing spline曲线的导数曲线？？

Question

对大量实验测得的数据点用Smoothing spline曲线拟合，然后怎样得出拟合后的曲线的导数曲线？

Answer 1

一个简单的办法是直接用三点差分求导，我给你一个函数

function [dydx,x,y] = threePiontDerivative(x,y)
x = x(:); y = y(:);
n = length(x); m = length(y);
if(n~=m)
    error('x和y长度不一致')
end
if(n<3)
    error('x和y长度不够，长度至少要有3')
end
dx = diff(x);
h = mean(dx);
if(max(abs(h-dx))>1e-12)
    warning('x不是等间距的')
    if(~all(dx>0))
       [x,idx] = sort(x); 
       y = y(idx);
    end
    [x,ia] = unique(x);
    y = y(ia);
    dydx =  nonUniformPoint(x,y);
else
    dydx = zeros(size(y));
    dydx(2:end-1) = (y(3:end)-y(1:end-2))/(2*h);%首端点
    dydx(1) = (-3*y(1)+4*y(2)-y(3))/(2*h);%中间点
    dydx(n) = (y(n-2)-4*y(n-1)+3*y(n))/(2*h);%末端点
end
end
function dy =  nonUniformPoint(x,y)
dy = zeros(size(y));
n = length(y);
dx = diff(x);
h1 = dx(1:end-1);
h2 = dx(2:end);
h11 = h1.*h1;
h22 = h2.*h2;
h12 = h1.*h2;
a = -h2./(h11+h12);
b = (h2-h1)./h12;
c = h1./(h22+h12);
dy(2:n-1) = a.*y(1:n-2)+b.*y(2:n-1)+c.*y(3:n); 
dy(1) = -(2*h1(1)+h2(1))/(h11(1)+h12(1))*y(1) +...
    (h1(1)+h2(1))/h12(1)*y(2) +...
    -h1(1)/(h22(1)+h12(1))*y(3);
dy(n) = (h2(end)/(h11(end)+h12(end)))*y(n-2) +...
    -(h1(end)+h2(end))/h12(end)*y(n-1) +...
    (h1(end)+2*h2(end))/(h22(end)+h12(end))*y(n);
end

然后调用：

x = 0:pi/3:2*pi;
y = sin(x);
xp = 0:pi/10:2*pi;
yp = spline(x,y,xp);
figure(1);clf;
plot(x,y,'ro',xp,yp,'b-'); 
dy = threePiontDerivative(xp,yp);
hold on 
plot(xp,dy,'m--')
legend('原数据','Spline后','Spline求导')

效果：

当然我还想到了一种更简单的办法！！！， 一并实现了

x = 0:pi/3:2*pi;
y = sin(x);
xp = 0:pi/10:2*pi;
pp = spline(x,y);
yp = ppval(pp, xp); % 分段得到的样条插值值
p = pp.coefs; % 获取分段系数
q = zeros(size(p,1),size(p,2)-1);
for i = 1:size(p,1)
    q(i,:) = polyder(p(i,:)); %直接对分段多项式系数求导
end
pp.coefs = q;  pp.order = 3; % 把pp更改为求导后的系数和阶数
dy = ppval(pp, xp); % 直接用ppval获取样条插值后的导数
figure(1);clf; % 画图
plot(x,y,'ro',xp,yp,'b-');hold on;
plot(xp,dy,'m--')
legend('原数据','Spline后','Spline求导')

效果：

堪称完美！