三次样条函数

plot(px,py,'m+') %绘制插值点
[Sn]=Hermite(px,py,0,0,2) %根据插值点及其值计算三次样条函数
for i=1:10 %采用自然边界
Sx=-6+i:0.01:-5+i;
Sy=polyval(Sn(i,:),Sx)
plot(Sx,Sy,'b-.')
end %绘制函数图像
s2=[1,-X(i)]; s2=conv(s2,s2); b=[-2,H(i)+2*X(i+1)]; s2=conv(s2,b);
s3=[1,-X(i+1)]; s3=conv(s3,s3); c=[1,-X(i)]; s3=conv(s3,c);
s4=[1,-X(i)]; s4=conv(s4,s4); d=[1,-X(i+1)]; s4=conv(s4,d);
for k=1:n
H(k)=X(k+1)-X(k);
end
for k=1:m-2
u(k)=H(k+1)/(H(k+1)+H(k));
l(k)=1-u(k);
end
for i=2:m-1
A(i,i-1)=l(i);
A(i,i)=2
A(i,i+1)=u(i);
end
mY=zeros(m,1);g=mY;
for i=1:m-1
s1=[1,-X(i+1)]; s1=conv(s1,s1); a=[2,H(i)-2*X(i)]; s1=conv(s1,a);
s2=[1,-X(i)]; s2=conv(s2,s2); b=[-2,H(i)+2*X(i+1)]; s2=conv(s2,b);
s3=[1,-X(i+1)]; s3=conv(s3,s3); c=[1,-X(i)]; s3=conv(s3,c);
《数值分析》实验报告：实验二
实验名称：三次样条函数的实现
实验地点：八教
所使用的工具软件及环境：Matlab
一、实验目的：
1．利用Matlab实现三次样条函数；
2．利用Matlab绘制函数图像；
二、实验内容：
1.根据所给定的节点值和节点上的函数值,取边界条件为“自然边界条件”,编制构造三次样条函数的Matlab程序;
2.给定函数f(x)=1/(1+x^2)，取插值节点为-5:1:5，利用编制的Matlab程序求出三次样条函数，并绘制出准确函数f(x)以及所求出的三次样条函数图像，进行比较。
3.自己设定一个函数,分别构造并绘制出该函数的Lagrange插值多项式和三次样条的函数图像,进行分析和比较。
三、操作步骤：
hold on
Sn
由上图可得：在边界两端采用二导数值时的拟合度比采用自然边界的拟合度高。
任课教师：鲍亮2010年10月6日
if flag==1
A(1,1)=1;A(m,m)=1;
for k=1:m
if k==1
g(1)=f0;
elseif k==m
g(m)=fn;
else
g(k)=3*((u(k-1).*(Y(k)-Y(k-1))./H(k-1))+l(k-1).*(Y(k+1)-Y(k))./H(k));
end
end
G(k+1)=gk;
end
G(m)=3*(Y(m)-Y(m-1))/H(m-1);l(n)=1;
for i=2:(m-1)
A(i,i)=2;A(m,m)=2;
A(i+1,i)=l(i);
A(i,i+1)=u(i);
end
mY=A\G';
syms x
m=length(X);Sn=zeros(m-1,4);
hold on
Sn
[Ln] = SF_Lagrange(px,py); %根据插值点及其值计算Lagrange插值多项式
Lx = -5:0.01:5; %绘图点
Ly = polyval(Ln,Lx); %计算绘图点上的多项式函数值
plot(Lx,Ly,'g--') %绘制多项式图像legend('sin(x)','point','Ln','Location','NorthEast')
Matlab程序如下：
一、构造三次样条函数的程序：
1、
a)用2个参数（针对自然边界）
function [Sn]=SF_Hermite(X,Y)
m=length(X);n=m-1;A=zeros(m,m);
H=zeros(1,n);
u=zeros(1,n);
G=zeros(1,n);
l=zeros(1,n);
plot(Sx,Sy,'b-.')
end%绘制函数图像
hold off
Sn
3：
functionSF_ZYT
x = -5 : 0.01 : 5;
y = sin(x);
plot(x,y,'r:') %绘制准确函数图像
hold on
px = -5 : 1 : 5; %插值点
py = sin(px); %插值点上的函数值
else
g(k)=3*((u(k-1).*(Y(k)-Y(k-1))./H(k-1))+l(k-1).*(Y(k+1)-Y(k))./H(k));
end
end
end
mY=inv(A)*g;
m=length(X);Sn=zeros(m-1,4);
for i=1:m-1
s1=[1,-X(i+1)]; s1=conv(s1,s1); a=[2,H(i)-2*X(i)]; s1=conv(s1,a);
end
if flag==2
A(1,1)=2;A(m,m)=2;
A(1,2)=1;A(m,m-1)=1;
for k=1:m
if k==1
g(1)=3*(Y(2)-Y(1))/(X(2)-X(1))-H(1)/2*f0;
elseif k==m
g(m)=3*(Y(m)-Y(m-1))/(X(m)-X(m-1))+H(m-1)/2*fn;
Sn(i,:)=Y(i)*s1/(H(i)^3)+Y(i+1)*s2/(H(i)^3)+mY(i)*s3/(H(i)^2)+mY(i+1)*s4/(H(i)^2);
end
二、主程序：
2：
function SF_ZYT
X = -5 : 0.1 : 5;
Y = 1./(1+X.^2);
plot(X,Y,'r:')
hold off
Ln
四、实验结果：
图Hale Waihona Puke 绘制如下：2：从图像中可以看出：对于该函数，采用自然边界的三次样条函数图像对于标准图像有相当好的拟合度。
3：
构造的函数为f(x)=sin(x)
标准函数图像和拉格朗日插值图像比较:
标准函数图像和三次样条图像比较：
标准、拉格朗日、三次样条图像比较：
结论：由上面3个图像可以看出：在边界处采用自然边界的三次样条函数的图像与朗格朗日插值图像相比，后者对于f(x)=sin(x)的拟合度高。
u(1)=1;A(1,1)=2;A(1,2)=u(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
uk=H(k)/(H(k)+H(k+1));u(k+1)=uk;
lk=1-u(k+1);l(k)=lk;
gk=3*((l(k).*(Y(k+1)-Y(k))./H(k))+(u(k+1).*(Y(k+2)-Y(k+1))./H(k+1)));
拓展思考：
在边界两端采用二阶导数值时的拟合情况：
[Sn]=Hermite(px,py,-sin(-5),-sin(5),2) %根据插值点及其值计算三次样条函数
for i=1:10
Sx=-6+i:0.01:-5+i;
Sy=polyval(Sn(i,:),Sx)
plot(Sx,Sy,'b-.')
end %绘制函数图像
hold on
pX = -5 : 1 : 5;
pY =1./(1+pX.^2) ;
plot(pX,pY,'m+')
[Sn]=Hermite(pX,pY,0,0,2) %根据插值点及其值计算三次样条函数
for i=1:10
Sx=-6+i:0.01:-5+i;
Sy=polyval(Sn(i,:),Sx)
s4=[1,-X(i)]; s4=conv(s4,s4); d=[1,-X(i+1)]; s4=conv(s4,d);
Sn(i,:)=Y(i)*s1/(H(i)^3)+Y(i+1)*s2/(H(i)^3)+mY(i)*s3/(H(i)^2)+mY(i+1)*s4/(H(i)^2);
end
b)用5个参数(针对两端一阶、二阶导数值已知情况)：
function [Sn]=SF_Hermite(X,Y,f0,fn,flag)
if length(X)~=length(Y)
error('变量不匹配');
end
m=length(X);n=m-1;A=zeros(m,m);