Question

如何用matlab编写三次样条差值的程序

Answer 1

function [y0,C,x]=Spline_interp_1(X,Y,s0,sN,x0)
%X,Y是已知插值点坐标
%s0,sN是两端点的一次导数值
%x0是插值点
%y0是三次样条函数在x0处的值
%C是分段三次样条函数的系数,按升幂排列,区间[Xj,Xj+1]上的表达式为C(j,1)+C(j,2)*(x-Xj)++C(j,3)*(x-Xj)^2+C(j,4)*(x-Xj)^3.
%x是用追赶法依下标顺利求得的插值节点X,Y处的二阶导数值
N=length(X);
C=zeros(4,N-1); h=zeros(1,N-1);
mu=zeros(1,N-1); lmt=zeros(1,N-1);
d=zeros(1,N); %d表示右端函数值
h=X(1,2:N)-X(1,1:N-1);
mu(1,N-1)=1; lmt(1,1)=1;
mu(1,1:N-2)=h(1,1:N-2)./(h(1,1:N-2)+h(1,2:N-1));
lmt(1,2:N-1)=h(1,2:N-1)./(h(1,1:N-2)+h(1,2:N-1));
d(1,1)=6*((Y(1,2)-Y(1,1))/h(1,1)-s0)/h(1,1);
d(1,N)=6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1);
d(1,2:N-1)=6*((Y(1,3:N)-Y(1,2:N-1))./h(1,2:N-1)-(Y(1,2:N-1)-Y(1,1:N-2))./h(1,1:N-2))./(h(1,1:N-2)+h(1,2:N-1));
%追赶法解三对角方程组
bit=zeros(1,N-1);
bit(1,1)=lmt(1,1)/2;
for i=2:N-1
bit(1,i)=lmt(1,i)/(2-mu(1,i-1)*bit(1,i-1));
end
y=zeros(1,N);
y(1,1)=d(1,1)/2;
for i=2:N
y(1,i)=(d(1,i)-mu(1,i-1)*y(1,i-1))/(2-mu(1,i-1)*bit(1,i-1));
end
x=zeros(1,N);
x(1,N)=y(1,N);
for i=N-1:-1:1
x(1,i)=y(1,i)-bit(1,i)*x(1,i+1);
end

v=zeros(1,N-1);
v(1,1:N-1)=(Y(1,2:N)-Y(1,1:N-1))./h(1,1:N-1);
C(1,:)=Y(1,1:N-1);
C(2,:)=v-h.*(2*x(1,1:N-1)+x(1,2:N))/6;
C(3,:)=x(1,1:N-1)/2;
C(4,:)=(x(1,2:N)-x(1,1:N-1))./(6*h);
if nargin<5
y0=0;
else
for j=1:N-1
if x0>=X(1,j)& x0<=X(1,j+1)
omg=x0-X(1,j);
y0=((C(4,j)*omg+C(3,j))*omg+C(2,j))*omg+C(1,j);
end
end
end
C=C'