sanci.m

三次样条插值模拟曲线

%function [m,H,lambda,mu,D,A,dY,sk]=NaturalSp(X,Y)
clc;
clear;
X=[0.9	1.3	1.9	2.1	2.6	3.0	3.9	4.4	4.7	5.0	6.0 7.0	8.0	9.2	10.5 11.3 11.6 12.0 12.6 13.0 13.3];
Y=[1.3	1.5	1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1	2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25];
m=length(X);A=zeros(m,m); n=m-1;H=zeros(1,n);lambda=zeros(1,n);
mu=zeros(1,n);lambda(1)=1;A(1,1)=2;A(1,2)=lambda(1);
D=zeros(1,n);H(1)=X(2)-X(1);mu(1)=1;D(1)=3*(Y(2)-Y(1));
for k=1:n
hk=X(k+1)-X(k);H(k+1)=hk;
end
H=H(2:n+1);
for k=1:n-1
lambdak=H(k)/(H(k)+H(k+1)); lambda(k+1)=lambdak;
muk=1-lambda(k+1);mu(k)=muk;
dk=3*((mu(k).*(Y(k+1)-Y(k))./H(k))+(lambda(k+1).*(Y(k+2)-Y(k+1))./H(k+1)));
D(k+1)=dk;
end
D(m)= 3*(Y(m)-Y(m-1))/H(m-1);mu(n)=1; n;H;lambda;mu;D;
for i=1:m-1
A(i,i)=2;A(m,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i);
end
dY=A\D';
plot(X,Y,'-g');
hold on;
for k=2:21
    x=X(k-1);
    sk1=0;
    for i=1:round((X(k)-X(k-1))/0.01+1)
        sk1(i)=Y(k-1)*((H(k-1)-2*X(k-1)+2*x)*(x-X(k))^2)/(H(k-1)^3)+Y(k)*((H(k-1)+2*X(k)-2*x)*(x-X(k-1))^2)/(H(k-1)^3)+dY(k-1)*((x-X(k-1))*(x-X(k))^2)/(H(k-1)^2)+dY(k)*((x-X(k))*(x-X(k-1))^2)/(H(k-1)^2);
        x=x+0.01;
    end
    plot(X(k-1):0.01:X(k),sk1);
    hold on;
end
legend('原始数据直接模拟曲线','三次自然样条插值曲线');
title('飞行中的野鸭的顶部轮廓曲线（三次自然样条插值方法）');
xlabel('x');
ylabel('f(x)');
hold off;