cs.txt

三次样条插值的源程序

function [S,YT]=yjc(X,Y,fb0,fe0,XX,a)
％其中X,Y矩阵分别为点(x(i),y(i)]
％fb0,feo表示当a=1是为第一边界条件否则为自然边界条件
K=length(XX);％待求点的个数
N=length(X);％输入点的个数
H=diff(X); ％建立H矩阵其元素H(i)=X(i+1)-X(i)
F=diff(Y)./H;％建立F矩阵其元素为F(i)=f[x(i),x(i+1)]
for i=2:N-1
D(i)=6*(F(i)-F(i-1))/(H(i-1)+H(i));％
U(i-1)=H(i-1)/(H(i-1)+H(i));
R(i)=H(i)/(H(i-1)+H(i));
end

if (a==1)％判断是第一边界条件还是自然边界条件
D(1)=6*(F(1)-fb0)/H(1);
D(N)=6*(fe0-F(N-1))/H(N-1);
U(N-1)=1;
R(1)=1;％建立矩阵U,R,D,为用追赶法求解做准备
else
R(1)=0;
U(N-1)=0;
D(1)=2*fb0;
D(N)=2*fe0;
end
％E=2*eye(N);
％U=diag(U,-1);
％R=diag(R,1);
％EE=E+R+U;
％下面为用追赶法求解M
B(1)=R(1)/2;
for i=2:N-1
B(i)=R(i)/(2-U(i-1)*B(i-1));
end
YY(1)=D(1)/2;
for i=2:N
YY(i)=(D(i)-U(i-1)*YY(i-1))/(2-U(i-1)*B(i-1));
end
M(N)=YY(N);
for i=N-1:-1:1
M(i)=YY(i)-B(i)*M(i+1);
end
％求出系数矩阵S
for k=1:N-1
S(k,1)=(M(k+1)-M(k))/(6*H(k));
S(k,2)=M(k)/2;
S(k,3)=(Y(k+1)-Y(k))/H(k)-H(k)*(2*M(k)+M(k+1))/6;
S(k,4)=Y(k);
end
for j=1:K
for i=1:N-1
if((XX(j)>=X(i))&(XX(j)<=X(i+1)))
YT(j)=polyval(S(i,:),XX(j)-X(i));
end
end
end
％处理字符串用作画图
t4='';
for i=1:N-1
t1=['x',int2str(i),'=X(',int2str(i),'):H(',int2str(i),')/100:X(',int2str(i),'+1);y',int2str(i),'=polyval(S(',int2str(i),',:),','x',int2str(i),'-X(',int2str(i),'));'];
执行(t1);
if(i==N-1)
t4=[t4,'x',int2str(i),',','y',int2str(i)];
else
t4=[t4,'x',int2str(i),',','y',int2str(i),','];
end
end
t5=['plot','(',t4,')'];