三次样条插值的Matlab实现(自然边界和第一边界条件).docx

PAGE

PAGE1/4

（第一边界条件）源代码：

functiony=yt1(x0,y0,f_0,f_n,x) (1)

%第一类边界条件下三次样条插值；

%xi所求点；

%yi所求点函数值；

%x已知插值点；

%y已知插值点函数值；

%f_0左端点一次导数值；

%f_n右端点一次导数值；

n=length(x0);z=length(y0);

h=zeros(n-1,1);

k=zeros(n-2,1);

l=zeros(n-2,1);S=2*eye(n);fori=1:n-1

h(i)=x0(i+1)-x0(i);end

fori=1:n-2

k(i)=h(i+1)/(h(i+1)+h(i));l(i)=1-k(i);

end

%对于第一种边界条件：

k=[1;k]; (2)

l=[l;1]; (3)

%构建系数矩阵S：fori=1:n-1

S(i,i+1)=k(i);

S(i+1,i)=l(i);

end

%建立均差表：F=zeros(n-1,2);fori=1:n-1

F(i,1)=(y0(i+1)-y0(i))/(x0(i+1)-x0(i));

end

D=zeros(n-2,1);fori=1:n-2

F(i,2)=(F(i+1,1)-F(i,1))/(x0(i+2)-x0(i));D(i,1)=6*F(i,2);

end

%构建函数D：

d0=6*(F(1,2)-f_0)/h(1); (4)

dn=6*(f_n-F(n-1,2))/h(n-1); (5)D=[d0;D;dn]; (6)

m=S\D;

%寻找x所在位置，并求出对应插值：fori=1:length(x)

forj=1:n-1

if(x(i)=x0(j+1))(x(i)=x0(j))

y(i)=(m(j)*(x0(j+1)-x(i))^3)/(6*h(j))+...

(m(j+1)*(x(i)-x0(j))^3)/(6*h(j))+...

(y0(j)-(m(j)*h(j)^2)/6)*(x0(j+1)-x(i))/h(j)+...

(y0(j+1)-(m(j+1)*h(j)^2)/6)*(x(i)-x0(j))/h(j);break;

elsecontinue;end

endend

（2）（自然边界条件）源代码：

仅仅需要对上面部分标注的位置做如下修改:

(1):functiony=yt2(x0,y0,x)

(2):k=[0;k]

(3):l=[l;0]

(4)+(5):删除

—(6):D=[0:D:0]