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

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.8

3、93.063.193.290.80.2算法描述：拉格朗日插值：其中是拉格朗日基函数，其表达式为：牛顿插值：其中三样条插值：所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(aX0X1Xn=j Y(i)=(Y(i)-Y(i-1)/(X(i)-X(i-j+1);else Y(i)=0;endend newt=newt,Y;end%计算牛顿插值 f=newt(1,2);for 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.mfunc

4、tion 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(1,n-2);for i=1:n%求函数的一阶差商h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i)/h(i);endfor i=2:n%求函数的二阶差商f2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);

5、d(i)=6*f2(i);endd(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(i);endA(1,2)=1;A(n+1,n)=1;for i=1:n+1 A(i,i)=2;endfor 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

6、)+(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);%三样条插值函数表达式endfor i=1:n disp(S(x)=);fprintf(%s (%d,%d)n,char(Sx(i),X(i),X(i+1);endfor i=1:nif xi=X(i)&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+

7、(M(i+1)-M(i)/(6*h(i)*(xi-X(i)3; endenddisp(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(牛顿插值结果);newton(X,Y,xi);disp(三次样条第一类边界条件插值结果);Threch1(X,Y,0.4

8、0,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=zeros(1,NUM);N=zeros(1,NUM);S=zeros(1,NUM);B=zeros(1,NUM);for i=1:NUM xi=a(i

9、); 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(a,S,r+);hold on;legend(原函数,拉格朗日插值,牛顿插值,三次样条插值,2);hold off7增加插值节点观察误差变化main4.mclear;clc;N=5;%4.5第一问Ini=zer

10、os(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(-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 on

11、plot(val,S,b)%三次样条插值函数图像str=sprintf(插值节点为%d时的插值效果,N);title(str); legend(原函数,拉格朗日插值,三次样条插值);%显示图例hold offfigureend8车门曲线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

12、,n-2);for i=1:nh(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i)/h(i);endfor i=2:nf2(i)=(f1(i)-f1(i-1)/(X(i+1)-X(i-1);d(i)=6*f2(i);endd(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(i);endA(1,2)=1;A(n+1,n)=1;for i=

13、1:n+1 A(i,i)=2;endfor i=2:n A(i,i-1)=B(i-1);A(i,i+1)=C(i-1);endM=Ad;x=zeros(1,n);S=zeros(1,n);for i=1:n x(i)=X(i)+0.5; S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x(i)-X(i)+M(i)/2*(x(i)-X(i)2+(M(i+1)-M(i)/(6*h(i)*(x(i)-X(i)3; endplot(X,Y,k); hold on;plot(x,S,o);title(三次样条插值效果图);legend(已知插值节点,三次样条插值); hold off实验结果：4.31计算插值节点处的函数值xi=0.13时Xi=0.36时2将多种插值函数即原函数图像画在同一图上4.5.1增加插值节点观察误差变化从上面三图可以看出增加插值节点并不能改善差之效果4.5.2 车门曲线