《优质实用文档精选——数值分析实验报告1.docx》由会员分享,可在线阅读,更多相关《优质实用文档精选——数值分析实验报告1.docx(23页珍藏版)》请在淘文阁 - 分享文档赚钱的网站上搜索。
1、实验报告实验项目名称 插值法 实验室数学实验室 所属课程名称 数值逼近 实 验 类 型 算法设计 实 验 日 期 班 级 学 号 姓 名 成 绩 实验概述:【实验目的及要求】本次实验的目的是熟练数值分析第二章“插值法”的相关内容,掌握三种插值方法:牛顿多项式插值,三次样条插值,拉格朗日插值,并比较三种插值方法的优劣。本次试验要求编写牛顿多项式插值,三次样条插值,拉格朗日插值的程序编码,并在MATLAB软件中去实现。【实验原理】数值分析第二章“插值法”的相关内容,包括:牛顿多项式插值,三次样条插值,拉格朗日插值的相应算法和相关性质。【实验环境】(使用的软硬件)软件:MATLAB 2012a硬件:
2、电脑型号:联想 Lenovo 昭阳E46A笔记本电脑操作系统:Windows 8 专业版 处理器:Intel(R)Core(TM)i3 CPU M 350 2.27GHz 2.27GHz实验内容:【实验方案设计】第一步,将书上关于三种插值方法的内容转化成程序语言,用MATLAB实现;第二步,分别用牛顿多项式插值,三次样条插值,拉格朗日插值求解不同的问题。【实验过程】(实验步骤、记录、数据、分析)实验的主要步骤是:首先分析问题,根据分析设计MATLAB程序,利用程序算出问题答案,分析所得答案结果,再得出最后结论。实验一:已知函数在下列各点的值为xi0.20.40.6.0.81.0f(xi)0.9
3、80.920.810.640.38试用4次牛顿插值多项式P4(x)及三次样条函数S(x)(自然边界条件)对数据进行插值。用图给出(xi,yi),xi=0.2+0.08i,i=0,1, 11, 10,P4(x)及S(x)。(1)首先我们先求牛顿插值多项式,此处要用4次牛顿插值多项式处理数据。已知n次牛顿插值多项式如下:Pn=f(x0)+fx0,x1(x-x0)+ fx0,x1,x2(x-x0) (x-x1)+ fx0,x1,xn(x-x0) (x-xn-1)我们要知道牛顿插值多项式的系数,即均差表中得部分均差。在MATLAB的Editor中输入程序代码,计算牛顿插值中多项式系数的程序如下:fun
4、ction varargout=newtonliu(varargin)clear,clcx=0.2 0.4 0.6 0.8 1.0;fx=0.98 0.92 0.81 0.64 0.38;newtonchzh(x,fx);function newtonchzh(x,fx)%由此函数可得差分表n=length(x);fprintf(*差分表*n);FF=ones(n,n);FF(:,1)=fx;for i=2:n for j=i:n FF(j,i)=(FF(j,i-1)-FF(j-1,i-1)/(x(j)-x(j-i+1); endendfor i=1:n fprintf(%4.2f,x(i);
5、 for j=1:i fprintf(%10.5f,FF(i,j); end fprintf(n);end由MATLAB计算得:xi f(xi) 一阶差商二阶差商三阶差商四阶差商0.200.9800.400.920-0.300000.600.810-0.55000-0.625000.800.640-0.85000-0.75000-0.208331.000.380-1.30000-1.12500-0.62500-0.52083所以有四次插值牛顿多项式为:P4(x)=0.98-0.3(x-0.2)-0.62500 (x-0.2)(x-0.4) -0.20833 (x-0.2)(x-0.4)(x-0
6、.6)-0.52083 (x-0.2)(x-0.4)(x-0.6)(x-0.8)(2)接下来我们求三次样条插值函数。用三次样条插值函数由上题分析知,要求各点的M值:三次样条插值函数计算的程序如下:function tgsanci(n,s,t) %n代表元素数,s,t代表端点的一阶导。x=0.2 0.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.38; n=5for j=1:1:n-1 h(j)=x(j+1)-x(j);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1 u(j)=1-r(j);endfo
7、r j=1:1:n-1 f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j);end d(1)=0 d(n)=0 a=zeros(n,n);for j=1:1:n a(j,j)=2;end r(1)=0; u(n)=0;for j=1:1:n-1 a(j+1,j)=u(j+1); a(j,j+1)=r(j);endb=inv(a)m=b*dp=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:n-1 p(j,1)=m(j)/(6*h(j); p(j,2)=m(j+1)/(6
8、*h(j); p(j,3)=(y(j)-m(j)*(h(j)2/6)/h(j); p(j,4)=(y(j+1)-m(j+1)*(h(j)2/6)/h(j);end p得到m=(0 -1.6071 -1.0714 -3.1071 0)T 即M0=0 ;M1= -1.6071;M2= -1.0714; M3= -3.1071; M4=0则根据三次样条函数定义,可得:S(x)= 接着,在Command Window里输入画图的程序代码,下面是画牛顿插值以及三次样条插值图形的程序:x=0.2 0.4 0.6 0.8 1.0;y=0.98 0.92 0.81 0.64 0.38;plot(x,y)hol
9、d on for i=1:1:5y(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4) -0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endk=0 1 10 11x0=0.2+0.08*kfor i=1:1:4y0(i)=0.98-0.3*(x(i)-0.2)-0.62500*(x(i)-0.2)*(x(i)-0.4) -0.20833*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)-
10、0.52083*(x(i)-0.2)*(x(i)-0.4)*(x(i)-0.6)*(x(i)-0.8)endplot( x0,y0,o,x0,y0 )hold ony1=spline(x,y,x0)plot(x0,y1,o)hold ons=csape(x,y,variational) fnplt(s,r)hold ongtext(三次样条自然边界)gtext(原图像)gtext(4次牛顿插值)运行上述程序可知:给出的(xi,yi),xi=0.2+0.08i,i=0,1, 11, 10点,S(x)及P4(x)图形如下所示:实验二:在区间-1,1上分别取用两组等距节点对龙格函数作多项式插值及三次
11、样条插值,对每个值,分别画出插值函数即的图形。我们先求多项式插值:在MATLAB的Editor中建立一个多项式的M-file,输入如下的命令(如牛顿插值公式):function C,D=newpoly(X,Y)n=length(X);D=zeros(n,n) D(:,1)=Y for j=2:n for k=j:n D(k,j)=(D(k,j-1)- D(k-1,j-1)/(X(k)-X(k-j+1); endendC=D(n,n);for k=(n-1):-1:1 C=conv(C,poly(X(k) m=length(C); C(m)= C(m)+D(k,k);end当n=10时,我们在C
12、ommand Window中输入以下的命令:clear,clf,hold on; X=-1:0.2:1; Y=1./(1+25*X.2); C,D=newpoly(X,Y); x=-1:0.01:1; y=polyval(C,x); plot(x,y,X,Y,.); grid on; xp=-1:0.2:1; z=1./(1+25*xp.2); plot(xp,z,r)得到插值函数和f(x)图形:当n=20时,我们在Command Window中输入以下的命令:clear,clf,hold on; X=-1:0.1:1; Y=1./(1+25*X.2); C,D=newpoly(X,Y); x
13、=-1:0.01:1; y=polyval(C,x); plot(x,y,X,Y,.); grid on; xp=-1:0.1:1; z=1./(1+25*xp.2); plot(xp,z,r)得到插值函数和f(x)图形:下面再求三次样条插值函数,在MATLAB的Editor中建立一个多项式的M-file,输入下列程序代码:function S=csfit(X,Y,dx0,dxn) N=length(X)-1;H=diff(X); D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N); C=H(2:N); U=6*diff(D);B(1)=B(1)-H(1)
14、/2;U(1)=U(1)-3*(D(1);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(-D(N); for k=2:N-1 temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end M(N)=U(N-1)/B(N-1);for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2;for k=0:N-1 S(k+1,1)=(M(
15、k+2)-M(k+1)/(6*H(k+1); S(k+1,2)=M(k+1)/2; S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2)/6; S(k+1,4)=Y(k+1);end当n=10时,我们在Command Window中输入以下的命令:clear,clcX=-1:0.2:1;Y=1./(25*X.2+1); dx0= 0.0739644970414201;dxn= -0.0739644970414201;S=csfit(X,Y,dx0,dxn) x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1);x2=-0.5:0.01:0
16、;y2=polyval(S(2,:),x2-X(2); x3=0:0.01:0.5; y3=polyval(S(3,:),x3-X(3);x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4); plot(x1,y1,x2,y2,x3,y3,x4,y4, X,Y,.)结果如图:当n=20时,我们在Command Window中输入以下的命令:clear,clcX=-1:0.1:1;Y=1./(25*X.2+1); dx0= 0.0739644970414201;dxn= -0.0739644970414201;S=csfit(X,Y,dx0,dxn) x1=-1:0.0
17、1:-0.5;y1=polyval(S(1,:),x1-X(1);x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2); x3=0:0.01:0.5; y3=polyval(S(3,:),x3-X(3);x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4); plot(x1,y1,x2,y2,x3,y3,x4,y4, X,Y,.)结果如图:实验三:下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似,在区间0,64上作图。(1)用这9各点作8次多项式插值L8(x).(2)用三次样条(自然边界条件)程序求
18、S(x)。从结果看在0,64上,那个插值更精确;在区间0,1上,两种哪个更精确?L8(x)可由公式Ln(x)=得出。 三次样条可以利用自然边界条件。写成矩阵:其中j=,i=,dj=6fxj-1,xj,xj+1,n=0=0 d0=dn=0l0(x)=l1(x)= l2(x)= l3(x)= l4(x)= l5(x)= l6(x)= l7(x)= l8(x)= L8(x)= l1(x)+2 l2(x)+3 l3(x)+4 l4(x)+5 l5(x)+6 l6(x)+7 l7(x)+8 l8(x)求三次样条插值函数由MATLAB计算:可得矩阵形式的线性方程组为:在MATLAB中的Editor中输入程
19、序代码,以下是三次样条函数的程序代码:function tgsanci(n,s,t) %n代表元素数,s,t代表端点的一阶导。y=0 1 2 3 4 5 6 7 8;x=0 1 4 9 16 25 36 49 64; n=9for j=1:1:n-1 h(j)=x(j+1)-x(j);endfor j=2:1:n-1 r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1 u(j)=1-r(j);endfor j=1:1:n-1 f(j)=(y(j+1)-y(j)/h(j);endfor j=2:1:n-1 d(j)=6*(f(j)-f(j-1)/(h(j-1)+h(j
20、);end d(1)=0 d(n)=0 a=zeros(n,n);for j=1:1:n a(j,j)=2;end r(1)=0; u(n)=0;for j=1:1:n-1 a(j+1,j)=u(j+1); a(j,j+1)=r(j);endb=inv(a)m=b*dt=ap=zeros(n-1,4); %p矩阵为S(x)函数的系数矩阵for j=1:1:n-1 p(j,1)=m(j)/(6*h(j); p(j,2)=m(j+1)/(6*h(j); p(j,3)=(y(j)-m(j)*(h(j)2/6)/h(j); p(j,4)=(y(j+1)-m(j+1)*(h(j)2/6)/h(j);en
21、d p解得:M0=0;M1=-0.5209;M2=0.0558;M3=-0.0261;M4=0.0006;M5=-0.0029;M6=-0.0008;M7=-0.0009;M8=0,则三次样条函数:S(x)= 下面进行画图,在Command Window中输入画图的程序代码:%画图形比较那个插值更精确的函数:x0=0 1 4 9 16 25 36 49 64;y0=0 1 2 3 4 5 6 7 8;x=0:64;y=lagr1(x0,y0,x);plot(x0,y0,o)hold onplot(x,y,r);hold on;pp=csape(x0,y0,variational)fnplt(p
22、p,g);hold on;plot(x0,y0,:b);hold on%axis(0 2 0 1); %看0 1区间的图形时加上这条指令gtext(三次样条插值)gtext(原图像)gtext(拉格朗日插值)%下面是求拉格朗日插值的函数function y=lagr1(x0,y0,x) n=length(x0); m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s; end拉格朗日插值函
23、数与三次样条插值函数如图中所示,绿色实线条为三次样条插值曲线,蓝色虚线条为x=y2的曲线,另外一条红色线条为拉格朗日插值曲线。图3-1为0 1的曲线,图3-2为0 64区间上的曲线。图31图32由图3-1可以看出,红色的线条与蓝色虚线条几乎重合,所以可知拉格朗日插值函数的曲线更接近开平方根的函数曲线,在0,1朗格朗日插值更精确。而在区间0,64上从图3-2中可以看出蓝色虚线条和绿色线条是几乎重合的,而红色线条在30,70之间有很大的振荡,所在在区间0,64三次样条插值更精确写。【结论】(结果)单个多项式高次插值效果并不理想,有龙格现象,偏差大,没有使用价值。而分段低次插值则精确度较高,拟合效果较好,而三次样条插值具有良好的收敛性与稳定性,与分段低次插值相比较光滑度更高,而且提供的信息也相对少一些。我们可以看到,在以上的三道实验题里,我们可以从图形中看出,三次样条的拟合程度是三种插值函数里最好的。【小结】通过此次实验,我对牛顿多项式插值,三次样条插值,拉格朗日插值有了更进一步的了解,知道了三次样条的拟合程度在高次的情况下更高,在理论上和应用上都有重要意义,在利用计算机编程软件进行高次插值的时候,我们可以多考虑利用三次样条进行插值。指导教师评语及成绩:成绩: 指导教师签名: 批阅日期:第 23 页 共 23 页
限制150内