yangtiaochazhi.m

该程序能实现第一型三次样条插值

function yangtiaochazhi(Px,Py,sd0,sdn,x)
%Px为插值节点,Py为插值节点对应函数值,sd0为断点0一阶导,sdn为端点n一阶导,x为待求的插值节点
n=length(Px);
n=n-1;

[Ph,Pd,M]=M_computation(Px,Py,sd0,sdn);
flag=0;

fprintf('\n第一型3次样条插值函数是:')
for j=1:n
    if Py(j)==0
        fprintf('\nS(%d)=',j);
    else fprintf('\nS(%d)=%f',j,Py(j));
    end
        
    Ss1=Pd(j)-(M(j)/3+M(j+1)/6)*Ph(j);
    if Ss1>0
        fprintf('+%f*(x-%f)',Ss1,Px(j));
    else if Ss1<0
            fprintf('%f*(x-%f)',Ss1,Px(j));
        end
    end
    
    Ss2=M(j)/2;
    if Ss2>0
        fprintf('+%f*(x-%f)^2',Ss2,Px(j));
    else if Ss2<0
            fprintf('%f*(x-%f)^2',Ss2,Px(j));
        end
    end
    
    Ss3=(M(j+1)-M(j))/(6*Ph(j));
    if Ss3>0
        fprintf('+%f*(x-%f)^3',Ss3,Px(j));
    else if Ss3<0
            fprintf('%f*(x-%f)^3',Ss3,Px(j));
        end
    end
   
    %fprintf('\nS(%d)=%f+%f*(x-%f)+%f*(x-%f)^2+%f*(x-%f)^3',j,Py(j),Ss1,Px(j),Ss2,Px(j),Ss3,Px(j));
    fprintf('     (%f=<x=<%f)',Px(j),Px(j+1));
    
    if x>=Px(j)
        if x<=Px(j+1)
            S=Py(j)+Ss1*(x-Px(j))+Ss2*(x-Px(j))^2+Ss3*(x-Px(j))^3;
            flag=1;
        end
    end
end
 
if flag==1
    fprintf('\n当x=%f时,S(%f)=%f',x,x,S);
else
    fprintf('\n\nx=%f,该点不处于插值节点的区间内',x);
end
    
    


function [Ph,Pd,M]=M_computation(Px,Py,s0,sn)
%计算M
[Ph,Pu,Pk,Pd,Pdd]=param_computation(Px,Py);
n=length(Px);
d0=(Py(2)-Py(1))/(Px(2)-Px(1));
d0=(d0-s0)/(Px(2)-Px(1));
dn=(Py(n)-Py(n-1))/(Px(n)-Px(n-1));
dn=(sn-dn)/(Px(n)-Px(n-1));

A=zeros(n,n);
b=zeros(n,1);
M=zeros(n,1);

A(1,1)=2;
A(1,2)=1;
b(1)=6*d0;
for i=2:n-1
    A(i,i)=2;
    A(i,i-1)=Pu(i-1);
    A(i,i+1)=Pk(i-1);
    b(i)=6*Pdd(i-1);
end
A(n,n)=2;
A(n,n-1)=1;
b(n)=6*dn;

M=A\b;




function [Ph,Pu,Pk,Pd,Pdd]=param_computation(Px,Py)
%计算参数
n=length(Px);
n=n-1;
Ph=zeros(1,n);
Pu=zeros(1,n-2);
Pk=zeros(1,n-2);
Pd=zeros(1,n);
Pdd=zeros(1,n-1);

for j=1:n
    Ph(j)=Px(j+1)-Px(j);
end
for j=1:n-1
    Pu(j)=Ph(j)/(Ph(j)+Ph(j+1));
    Pk(j)=1-Pu(j);
end

Pd=chashang(Px,Py,1);
Pdd=chashang(Px,Pd,2);


function D=chashang(X,Y,i)
n=length(Y);
n=n-1;
D=zeros(1,n);

for j=1:n
    D(j)=(Y(j+1)-Y(j))/(X(j+i)-X(j));
end