西北农林科技大学数值分析数值法实验报告

《西北农林科技大学数值分析数值法实验报告》由会员分享，可在线阅读，更多相关《西北农林科技大学数值分析数值法实验报告（17页珍藏版）》请在装配图网上搜索。

1、数值法实验报告专业班级：信息与计算科学 121 姓名：金辉 学号：20120142801) 实验目的 本次实验的目的是熟练数值分析第二章“插值法”的相关内容，掌握三 种插值方法：牛顿多项式插值，三次样条插值，拉格朗日插值，并比较三种插值 方法的优劣。本次试验要求编写牛顿多项式插值，三次样条插值，拉格朗日插值的程序编码， 并在MATLAB软件中去实现。2) 实验题目实验一：已知函数在下列各点的值为试用4次牛顿插值多项式P (x)及三次样条函数S (x)(自然边界条件)4据进行插值。用图给出(x，y)，x=0.2+0.08i,i ii及S (x)。i=0,1, 11, 10,P4对数x)x0.20

2、.40.6.0.81.0f (x)0.980.920.810.640.38实验二：在区间-1,1上分别取n = 10,20用两组等距节点对龙格函数f (x) =1 作多1 + 25 x 2 项式插值及三次样条插值，对每个n值，分别画出插值函数即f (X)的图形。实验三：下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似，在区间0,64上作图。(1) 用这9各点作8次多项式插值L (x).8(2) 用三次样条(自然边界条件)程序求S (x)。从结果看在0,64上,那个插值更精确；在区间0,1上，两种哪个更精确？3) 实验原理与理论基础数值分析第二章“插值法

3、”的相关内容，包括：牛顿多项式插值，三次 样条插值，拉格朗日4) 实验内容实验一：x0.20.40.6.0.81.0f (x )0.980.920.810.640.38试用4次牛顿插值多项式P (x)及三次样条函数S (x)4自然边界条件)对数1, 11, 10,P(x)4已知函数在下列各点的值为据进行插值。用图给出(x，y)，x=0.2+0.08i, i=0,i ii及S (x)。(1)首先我们先求牛顿插值多项式,此处要用4次牛顿插值多项式处理数据。 已知n次牛顿插值多项式如下：P =f(x )+fx ,x (x-x )+ fx ,x ,x (x-x ) (x-x )+ +n 0010012

4、01fx ,x(XX ) (x-x )01n 0n-1我们要知道牛顿插值多项式的系数,即均差表中得部分均差。在MATLAB的Editor中输入程序代码，计算牛顿插值中多项式系数的程序如 下：function 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);fpri ntf(*差分表 *n);FF=ones(n,n);FF(:,1)=fx;

5、for i=2:nfor j=i:n FF(j,i)=(FF(j,i-1)-FF(j-1,i-1)/(x(j)-x(j-i+1);endendfor i=1:nfprintf(%4.2f,x(i);for j=1:ifprintf(%10.5f,FF(i,j);endfprintf(n);end由MATLAB计算得:xf(x.)一阶差商二阶差商三阶差商四阶差商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

6、.52083所以有四次插值牛顿多项式为：P (x)=0.98-0.3(x-0.2)-0.62500 (x-0.2)(x-0.4) -0.208334(x-0.2)(x-0.4)(x-0.6)-0.52083 (x-0.2)(x-0.4)(x-0.6)(x-0.8)(2)接下来我们求三次样条插值函数。用三次样条插值函数由上题分析知，要求各点的M值:_ 20000 _ M _0-0 _0.50020.50000M1-3.75000.50020.5000M=-4.500000.50020.500M2-6.7500000023M40三次样条插值函数计算的程序如下:function tgsanci(n,

7、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-1h(j)=x(j+1)-x(j);endfor j=2:1:n-1r(j)=h(j)/(h(j)+h(j-1);endfor j=1:1:n-1u(j)=1-r(j);endfor j=1:1:n-1f(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);endd(1)=0d(n)=0a=zeros(n,n);for j=1:

8、1:na(j,j)=2;endr(1)=0;u(n)=0;for j=1:1:n-1a(j+1,j)=u(j+1); a(j,j+1)=r(j);endb=inv(a)m=b*dp二zeros(n-l,4); %p矩阵为S(x)函数的系数矩阵 for j=1:1:n-1p(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+l)-m(j+1)* (h(j厂2/6)/h(j); endp得到 m=(0 -1.6071-1.0714 -3.10710)t即 M=0 ;M=

9、 -1.6071;M= -1.0714; M= -3.1071; M =001234则根据三次样条函数定义，可得：0(0.4 x) 1.3393(x0.2) + 4.90(0.4 x) + 4.6536(x0.2), x w0.2,0.4S(x)二 I -1.33930.6x)3 0.8929(x-0.4) + 4.653(0.6x)+4.0857(x0.4), xw0.4,0.60.892(0.8x)3-2.589(x-0.6)3 + 4.0857(0.8x) + 3.303(x0.6), x g 0.6,0.8、-2.58931.0x)3 -0(x0.8) + 3.3036(1.0x) +

10、1.9(x0.8), x g0.8,1.0接着，在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)hold onfor 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

11、=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)-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(三次

12、样条自然边界)gtext(原图像)gtext(4次牛顿插值)运行上述程序可知：给出的 (x,y), x=0.2+0.08i, i=0, 1, 11, 10点，Si ii(x)及P (x)图形如下所示：4实验二：在区间-1,1上分别取N = 10,20用两组等距节点对龙格函数/(x)=- 作多 1 + 25 x 2项式插值及三次样条插值，对每个n值，分别画出插值函数即/(X)的图形。我们先求多项式插值：在MATLAB的Editor中建立一个多项式的M-file,输入如下的命令(如牛顿插值公式)：function C,D=newpoly(X,Y)n=length(X);D=zeros(n,n)D(

13、:,1)=Yfor j=2:nfor k=j:nD(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:1C=conv(C,poly(X(k)m=length(C);C(m)= C(m)+D(k,k);end当n=10时，我们在Command Window中输入以下的命令: clear,clf,hold on;X=-l:0.2:l;Y=l./(l+25 *X2);C,D二newpoly(X,Y);x=-l:0.01:l;y=polyval(C,x);plo t(x,y,X,Y,.);grid on;xp

14、=-l:0.2:l;z=l./(l+25 *xp.2);plot( xp,z,r)得到插值函数和f (x)图形：当n=20时，我们在Command Window中输入以下的命令: clear,clf,hold on;X=-l:0.1:l;Y=l./(l+25 *X2); C,D=newpoly(X,Y);x=l:0.01:l;y二polyval(C,x);plo t(x,y,X,Y,.);grid on;xp=-l:0.1:l;z=l./(l+25 *xp.2);plot( xp,z,r)得到插值函数和f (x)图形:下面再求三次样条插值函数，在MATLAB的Editor中建立一个多项式的M-

15、file, 输入下列程序代码：func tion S=csf it( X,Y,dxO,dxn)N=leng th(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)/2;U(1)=U(1)-3 *(D(1);B(N-l)=B(N-l)-H(N)/2;U(Nl)=U(Nl)3 *(-D(N);for k=2:N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);endM(N)=

16、U(N-1)/B(N-1);for k=N-2:-1:1M(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-1S(k+1,1)=(M(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

17、=l./(25 *X2+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;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时，我们在C

18、ommand Window中输入以下的命令: clear,clcX=-l:0.1:l;Y=l./(25 *X2+1);dx0= 0.0739644970414201;dxn二-0.0739644970414201;S=csf it( X,Y,dx0,dxn)x1=1:0.01：0.5;y1二polyval(S(1,:),x1X(1);x2=-0.5:0.01:0;y2二polyval(S(2,:),x2-X(2);x3=0:0.01:0.5; y3=polyval(S(3,:),x3X(3);x4=0.5:0.01:1;y4二polyval(S(4,:),x4-X(4);plo t(x1,y1

19、,x2,y2,x3,y3,x4,y4, X,Y,.)结果如图：实验三：下列数据点的插值x01491625364964y012345678可以得到平方根函数的近似，在区间0,64上作图。(1) 用这9各点作8次多项式插值L (x).8(2) 用三次样条(自然边界条件)程序求S (x)。从结果看在0,64上,那个插值更精确；在区间0,1上，两种哪个更精确？L (x)可由公式L (x) = E y l (x )得出。8 n k k jk=0三次样条可以利用自然边界条件。写成矩阵其中卩j= hTj-1 j2 九0卩 210 卩00200000 _ M _d _00九00Md11i2九0M=dP22九M

20、2d233330P42M4d4h九 一i= h + hj-1jdj=6fx ,x,x ，j-1 j j+1=九=0 d =d =0n 0 0 nl(x)=1(x - 0)(x 4)(x 9)(x 16)(x 25)(x 36)(x 49)(x 64) (1 0)(1 4)(1 9)(1 16)(1 25)(1 36)(1 49)(1 64)l(x)=2(x - 0)(x 1)(x 9)(x 16)(x 25)(x 36)(x 49)(x 64) (4 0)(4 1)(4 9)(4 16)(4 25)(4 36)(4 49)(4 64)l(x)=3(x - 0)(x 1)(x 4)(x 16)(

21、x 25)(x 36)(x 49)(x 64) (9 0)(9 1)(9 4)(9 16)(9 25)(9 36)(9 49)(9 64)l(x)=(x - 0)(x 1)(x 4)(x 9)(x 25)(x 36)(x 49)(x 64)4(16 0)(16 1)(16 4)(16 9)(16 25)(16 36)(16 49)(16 64)l(x)=(x - 0)(x -1)(x 4)(x 9)(x 16)(x 36)(x 49)(x 64)5(25 0)(25 1)(25 4)(25 9)(25 16)(25 36)(25 49)(25 64)l(x)=(x - 0)(x 1)(x 4)

22、(x 9)(x 16)(x 25)(x 49)(x 64)6(36 0)(36 1)(36 4)(36 9)(36 16)(36 25)(36 49)(36 64)l(x)=(x - 0)(x -1)(x 4)(x 9)(x 16)(x 25)(x 36)(x 64)7(49 0)(49 1)(49 4)(49 9)(49 16)(49 25)(49 36)(49 64)l(x)=(x - 0)(x -1)(x 一 4)(x 一 9)(x 一 16)(x 一 25)(x 一 36)(x 一 49)8(64 0)(64 1)(64 4)(64 9)(64 16)(64 25)(64 36)(64

23、 49)l (x)= (x _ 1)(x _ 4)(x _ 9)(x _ 16)(x _0(0_1)(0_4)(0_9)(0_16)(0_25)(x _25)(0 _36)( x 49)( x 64)36)(0 49)(0 64)L(x)= l(x)+2 l(x)+3 l(x)+4 l(x)+5 l(x)+6 l(x)+7 l(x)+8 l(x)812345678求三次样条插值函数由MATLAB计算:可得矩阵形式的线性方程组为：-2.000000000000 一M 一 0 -0.25002.00000.75000000000M1.000.37502.00000.6250000001M0.100

24、0.41672.00000.58330000M230.02860000.43752.00000.5625000M4=0.011900000.45002.00000.550000M50.0061000000.45832.00000.54170M60.00350000000.46342.00000.5357M70.0022_ 000000002.0000_M8_ 0 _在MATLAB中的Editor中输入程序代码， 以下是三次样条函数的程序代码： function tgsanci(n,s, t) %n代表元素数，s, t代表端点的一阶导。y=0 1 2 3 4 5 6 7 8; x=0 1 4 9

25、 16 25 36 49 64;n=9 for 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);end d(1)=0 d(n)=0 a=zeros(n,n);for j=1:1:n a(j,j)=2;endr(1)=0;u(n)=0;for j=1:1:n-1 a(j+1,j

26、)=u(j+1); a(j,j+1)=r(j);endb=inv(a)m=b*dt=ap二zeros(n-l,4); %p矩阵为S(x)函数的系数矩阵 for j=1:1:n-1p(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+l)-m(j+l) *(h(j厂2/6)/h(j); endp解得:M=0;M=-0.5209;M=0.0558;M=-0.0261;M=0.0006;M=-0.0029;M=-0.0008;M=-01234567-0.0009;M =

27、0，则三次样条函数：80(1 x)3 +0.0868(x 0)3 + 0(1 X)+ 1.0868(x -0), x g 0,1-0.0289(4 x )3 -0.0031( x-1)3 + 0.5938(4 x) + 0.6388( x 1), x g 1,4 0.0019(9x)30.0009(x-4)3 +0.3535(9x)+0.6271(x4), xg4,9S(x)=-0.0006(16 x)3 + 0(x 9) + 0.4590(16 x)-0.5708(x -9), x g 9,160 s0(25 x)-0.0001(x 16)-0.4436(25 x) + 0.5600(x -

28、16), x g 16,25 0(36 x)3 - 0(x 25)3 + 0.4599(36 x) + 0.5470(x -25), x g 25,360(48 x)3 - 0(x 36)3 + 0.4633(48 x) + 0.5404(x -36), x g 36,480(64 x)3 + 0(x 48)3 + 0.4689(64 x) + 0.5333(x -48), x g 48,64F面进行画图，在Command Window中输入画图的程序代码:%画图形比较那个插值更精确的函数：x0=0 1 4 9 16 25 36 49 64;y0=0 1 2 3 4 5 6 7 8;x=0:6

29、4;y=lagr1(x0,y0,x);plot(x0,y0,o)hold on plot(x,y,r);hold on;pp=csape(x0,y0,variational) fnplt(pp,g);hold on; plot(x0,y0,:b);hold on %axis(0 2 0 1); %看0 1区间的图形时加上这条指令 gt ext(三次样条插值)gtext(原图像)gtext(拉格朗日插值)%下面是求拉格朗日插值的函数function y=lagr1(x0,y0,x) n=length(x0); m=length(x);for i=1:mz=x(i);s=0.0;for k=1:n

30、p=1.0;for j=1:n if j=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(i)=s;end 拉格朗日插值函数与三次样条插值函数如图中所示，绿色实线条为三次样条 插值曲线，蓝色虚线条为x=y2的曲线，另外一条红色线条为拉格朗日插值曲线。 图3-1为0 1的曲线，图3-2为0 64区间上的曲线。图31图32由图3-1可以看出，红色的线条与蓝色虚线条几乎重合，所以可知拉格朗日 插值函数的曲线更接近开平方根的函数曲线，在0,1朗格朗日插值更精确。而 在区间0, 64上从图3-2中可以看出蓝色虚线条和绿色线条是几乎重合的，而 红色线条在

31、30, 70之间有很大的振荡，所在在区间0, 64三次样条插值更精 确写。5）实验结果单个多项式高次插值效果并不理想，有龙格现象，偏差大，没有使用价值。 而分段低次插值则精确度较高，拟合效果较好，而三次样条插值具有良好的收敛 性与稳定性，与分段低次插值相比较光滑度更高，而且提供的信息也相对少一些。 我们可以看到，在以上的三道实验题里，我们可以从图形中看出，三次样条的拟 合程度是三种插值函数里最好的。6）实验结果分析与小结通过此次实验，我对牛顿多项式插值，三次样条插值，拉格朗日插值有了更 进一步的了解，知道了三次样条的拟合程度在高次的情况下更高，在理论上和应 用上都有重要意义，在利用计算机编程软件进行高次插值的时候，我们可以多考 虑利用三次样条进行插值。