数值分析作业三次样条插值

1、数值计算方法作业实验名称实验4.3三次样条插值函数（P126）4.5三次样条插值函数的收敛性（P127）实验时间姓名班级学号成绩实验4.3 三次样条差值函数实验目的：掌握三次样条插值函数的三弯矩方法。实验函数:x0.00.10.20.30.4F(x)0.50000.53980.57930.61790.7554求f(0.13)和f(0.36)的近似值实验内容：（1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件；（2） 计算各插值节点的弯矩值；（3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线比较插值结果。实验4.5 三次样条差值函数的收敛性实验目的：多项式

2、插值不一定是收敛的，即插值的节点多，效果不一定好。对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。实验内容：按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。实验要求：（1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；（2） 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子，考虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一段数据如下：0123456789100.00.791.532.192.713.033.272.893.

3、063.193.290.80.2算法描述：拉格朗日插值： 其中是拉格朗日基函数，其表达式为：牛顿插值：其中三样条插值：所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a<X0<X1<Xn<b)分成的每个小区间xi-1,xi上是三次多项式，其在此区间上的表达式如下：式中Mi=.因此，只要确定了Mi的值，就确定了整个表达式，Mi的计算方法如下：令则Mi满足如下n-1个方程：常用的边界条件有如下几类：（1） 给定区间两端点的斜率m0,mn,即（2） 给定区间两端点的二阶导数M0，Mn,即（3） 假设y=f(x)是以b-a为周期的周期函数，则要求三次样条插值函数S

4、（x）也为周期函数，对S（x）加上周期条件对于第一类边界条件有对于第二类边界条件有其中那么解就可以为 对于第三类边界条件，由此推得，其中，那么解就可以为：程序代码：1拉格朗日插值函数Lang.mfunction f=lang(X,Y,xi)%X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%f求得的拉格朗日插值多项式的值 n=length(X); f=0; for i=1:n l=1; for j=1:i-1 l=l.*(xi-X(j)/(X(i)-X(j); end; for j=i+1:n l=l.*(xi-X(j)/(X(i)-X(j); end;%拉格朗日基函数 f=f

5、+l*Y(i); end fprintf('%dn',f)return2 牛顿插值函数newton.mfunction f=newton(X,Y,xi)%X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%f求得的拉格朗日插值多项式的值n=length(X);newt=X',Y'%计算差商表for j=2:n for i=n:-1:1 if i>=j Y(i)=(Y(i)-Y(i-1)/(X(i)-X(i-j+1); else Y(i)=0; end end newt=newt,Y'end %计算牛顿插值 f=newt(1,2); f

6、or i=2:n z=1; for k=1:i-1 z=(xi-X(k)*z; end f=f+newt(i-1,i)*z; end fprintf('%dn',f) return3三次样条插值第一类边界条件Threch.mfunction S=Threch1(X,Y,dy0,dyn,xi) % X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值 %dy0左端点处的一阶导数% dyn右端点处的一阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros

7、(1,n-2);for i=1:n%求函数的一阶差商 h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i)/h(i);end for i=2:n%求函数的二阶差商 f2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);d(i)=6*f2(i);end d(1)=6*(f1(1)-dy0)/h(1); d(n+1)=6*(dyn-f1(n-1)/h(n-1);%¸赋初值A=zeros(n+1,n+1); B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1 B(i)=h(i)/(h(i)+h(i+1);C(i)=1-B

8、(i);end A(1,2)=1;A(n+1,n)=1;for i=1:n+1 A(i,i)=2;end for i=2:n A(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=Ad;syms x;for i=1:n Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x-X(i). +M(i)/2*(x-X(i)2+(M(i+1)-M(i)/(6*h(i)*(x-X(i)3);digits(4); Sx(i)=vpa(Sx(i);%三样条插值函数表达式end for i=1:n disp('S(x)=');

9、fprintf('%s (%d,%d)n',char(Sx(i),X(i),X(i+1);end for i=1:n if xi>=X(i)&&xi<=X(i+1) S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(xi-X(i)+M(i)/2*(xi-X(i)2+(M(i+1)-M(i)/(6*h(i)*(xi-X(i)3; end enddisp('xi S');fprintf('%d,%dn',xi,S);return 4 三次样条插值第二类边界条件Threch2.mfunction S

10、x=Threch2(X,Y,d2y0,d2yn,xi) X为已知数据的横坐标%Y为已知数据的纵坐标%xi插值点处的横坐标%S求得的三次样条插值函数的值 %d2y0左端点处的二阶导数% d2yn右端点处的二阶导数n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:n%求一阶差商h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i)/h(i);end for i=2:n%求二阶差商f2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);d(i)=6*

11、f2(i);end d(1)=2*d2y0; d(n+1)=2*d2yn;%赋初值 A=zeros(n+1,n+1); B=zeros(1,n-1);C=zeros(1,n-1);for i=1:n-1 B(i)=h(i)/(h(i)+h(i+1);C(i)=1-B(i);end A(1,2)=0;A(n+1,n)=0;for i=1:n+1 A(i,i)=2;end for i=2:n A(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=Ad;syms x;for i=1:n Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i

12、)*(x-X(i). +M(i)/2*(x-X(i)2+(M(i+1)-M(i)/(6*h(i)*(x-X(i)3);digits(4); Sx(i)=vpa(Sx(i);end for i=1:n disp('S(x)=');fprintf('%s (%d,%d)n',char(Sx(i),X(i),X(i+1);end for i=1:n if xi>=X(i)&&xi<=X(i+1) S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(xi-X(i)+M(i)/2*(xi-X(i)2+(M(i+1)

13、-M(i)/(6*h(i)*(xi-X(i)3; end enddisp('xi S');fprintf('%d,%dn',xi,S);return5插值节点处的插值结果main3.mclearclcX=0.0,0.1,0.2,0.3,0.4;Y=0.5000,0.5398,0.5793,0.6179,0.7554;xi=0.13;%xi=0.36;disp('xi=0.13');%disp('xi=0.36');disp('拉格朗日插值结果');lang(X,Y,xi);disp('牛顿插值结果'

14、);newton(X,Y,xi);disp('三次样条第一类边界条件插值结果');Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分别为两端点处的一阶导数disp('三次样条第二类边界条件插值结果');Threch2(X,Y,0,-0.136,xi);%0,-0.136分别为两端点处的二阶导数6将多种插值函数即原函数图像画在同一张图上main2.mclearclcX=0.0,0.1,0.2,0.3,0.4;Y=0.5000,0.5398,0.5793,0.6179,0.7554;a=linspace(0,0.4,21);NUM=21;L=z

15、eros(1,NUM);N=zeros(1,NUM);S=zeros(1,NUM);B=zeros(1,NUM);for i=1:NUM xi=a(i); L(i)=lang(X,Y,xi);% 拉格朗日插值 N(i)=newton(X,Y,xi);% 牛顿插值 B(i)=normcdf(xi,0,1);%原函数 S(i)=Threch1(X,Y,0.4,0.36,xi);%三次样条函数第一类边界条件endplot(a,B,'-r');hold on;plot(a,L,'b');hold on;plot(a,N,'r');hold on;plot

16、(a,S,'r+');hold on;legend('原函数','拉格朗日插值','牛顿插值','三次样条插值',2);hold off7增加插值节点观察误差变化main4.mclear;clc;N=5;%4.5第一问Ini=zeros(1,1001);a=linspace(-1,1,1001);Ini=1./(1+25*a.2);for i=1:3 %节点数量变化次数 N=2*N; t=linspace(-1,1,N+1);%插值节点 ft=1./(1+25*t.2);%插值节点函数值 val=linspace(

17、-1,1,101); for j=1:101 L(j)=lang(t,ft,val(j); S(j)=Threch1(t,ft,0.074,-0.074,val(j);%三样条第一类边界条件插值 endplot(a,Ini,'k')%原函数图象hold onplot(val,L,'r')%拉格朗日插值函数图像hold onplot(val,S,'b')%三次样条插值函数图像str=sprintf('插值节点为%d时的插值效果',N);title(str); legend('原函数','拉格朗日插值'

18、,'三次样条插值');%显示图例hold off figureend8车门曲线main5.mclearclcX=0,1,2,3,4,5,6,7,8,9,10;Y=0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29;dy0=0.8;dyn=0.2;n=length(X)-1;d=zeros(n+1,1);h=zeros(1,n-1);f1=zeros(1,n-1);f2=zeros(1,n-2);for i=1:nh(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i)/h(i);end for i=2:nf2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);d(i)=6*f2(i);end d(1)=6*(f1(1)-dy0)/h(1); d(n+1)=6*(dyn-f