数值分析上级报告

1、数值分析上机实验报告函数插值方法课题五 函数插值方法一、 问题提出对于一元函数的n+1个节点值 j=(0,1,n)。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。数据如下：（1）0.40.550.650.800.951.050.410750.578150.696750.901.001.25382求五次Lagrange多项式，和分段三次插值多项式，计算,的值。（2）xj1234567yj0.3680.1350.0500.0180.0070.0020.001试构造Lagrange多项式，计算的值。0.165299 0.00213348二、要求 1、利用Lagrange

2、插值公式编写出插值多项式程序；2、给出插值多项式或分段三次插值多项式的表达式； 3、根据节点选取原则，对问题（2）用三点插值或二点插值，其结果如何； 4、对此插值问题用Newton插值多项式其结果如何。三、目的和意义 1、学会常用的插值方法，求函数的近似表达式，以解决其它实际问题； 2、明确插值多项式和分段插值多项式各自的优缺点； 3、熟悉插值方法的程序编制； 4、如果绘出插值函数的曲线，观察其光滑性。 四、实验步骤1、根据（1）给出的实验数据，编写MATLAB程序，求出五次Lagrange多项式L5(x)，并计算,的值,画出的函数图。%五次lagrange 插值syms L5 Xx=0.4

3、0.55 0.65 0.80 0.95 1.05;y=0.41075 0.57815 0.69675 0.90 1.00 1.25382;L5=y(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(5)*(X-x(6)/(x(1)-x(2)*(x(1)-x(3)*(x(1)-x(4)*(x(1)-x(5)*(x(1)-x(6).+y(2)*(X-x(1)*(X-x(3)*(X-x(4)*(X-x(5)*(X-x(6)/(x(2)-x(1)*(x(2)-x(3)*(x(2)-x(4)*(x(2)-x(5)*(x(2)-x(6).+y(3)*(X-x(1)*(X-x(2)*(X-x(

4、4)*(X-x(5)*(X-x(6)/(x(3)-x(1)*(x(3)-x(2)*(x(3)-x(4)*(x(3)-x(5)*(x(3)-x(6).+y(4)*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(5)*(X-x(6)/(x(4)-x(1)*(x(4)-x(2)*(x(4)-x(3)*(x(4)-x(5)*(x(4)-x(6).+y(5)*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(6)/(x(5)-x(1)*(x(5)-x(2)*(x(5)-x(3)*(x(5)-x(4)*(x(5)-x(6).+y(6)*(X-x(1)*(X-x(2)*(X

5、-x(3)*(X-x(4)*(X-x(5)/(x(6)-x(1)*(x(6)-x(2)*(x(6)-x(3)*(x(6)-x(4)*(x(6)-x(5);L5x=vpa(L5)X=0.596;f1=subs(L5x)X=0.99;f2=subs(L5x)ezplot(L5x,0,2),title(五次lagrange插值函数),xlabel(X),ylabel(L5(X), hold on plot(x,y,r*) 计算结果为：L5x =385.7933077*(X - 0.95)*(X - 0.4)*(X - 0.8)*(X - 0.65)*(X - 0.55) + 770.86666666

6、666666666666666666667*(X - 0.95)*(X - 0.4)*(X - 0.8)*(X - 0.65)*(X - 1.05) - 1548.3333333333333333333333333333*(X - 0.95)*(X - 0.4)*(X - 0.8)*(X - 0.55)*(X - 1.05) + 1600.0*(X - 0.95)*(X - 0.4)*(X - 0.65)*(X - 0.55)*(X - 1.05) - 76.596736596736596736596736596737*(X - 0.95)*(X - 0.8)*(X - 0.65)*(X -

7、0.55)*(X - 1.05) - 1010.11101*(X - 0.4)*(X - 0.8)*(X - 0.65)*(X - 0.55)*(X - 1.05)f1 =0.6257f2 =1.0542所以：f(0.596)= f1= 0.6257； f(0.99)= f2 =1.0542图中红色点为（1）中给出的实验数据，从图中可以看出，实验数据点基本分布在插值函数的曲线上。2、对实验数据（2），编写MATLAB程序，求出六次Lagrange多项式L6(x)，并计算f(1.8),f(6.15)的值，画出L6(x)的函数图。%六次lagrange 插值syms L6 Xx=1 2 3 4 5

8、 6 7;y=0.368 0.135 0.050 0.018 0.007 0.002 0.001;L6=y(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(5)*(X-x(6)*(X-x(7)/(x(1)-x(2)*(x(1)-x(3)*(x(1)-x(4)*(x(1)-x(5)*(x(1)-x(6)*(x(1)-x(7).+y(2)*(X-x(1)*(X-x(3)*(X-x(4)*(X-x(5)*(X-x(6)*(X-x(7)/(x(2)-x(1)*(x(2)-x(3)*(x(2)-x(4)*(x(2)-x(5)*(x(2)-x(6)*(x(2)-x(7).+y(3)*(X-

9、x(1)*(X-x(2)*(X-x(4)*(X-x(5)*(X-x(6)*(X-x(7)/(x(3)-x(1)*(x(3)-x(2)*(x(3)-x(4)*(x(3)-x(5)*(x(3)-x(6)*(x(3)-x(7).+y(4)*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(5)*(X-x(6)*(X-x(7)/(x(4)-x(1)*(x(4)-x(2)*(x(4)-x(3)*(x(4)-x(5)*(x(4)-x(6)*(x(4)-x(7).+y(5)*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(6)*(X-x(7)/(x(5)-x(1)*(x(5

10、)-x(2)*(x(5)-x(3)*(x(5)-x(4)*(x(5)-x(6)*(x(5)-x(7).+y(6)*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(5)*(X-x(7)/(x(6)-x(1)*(x(6)-x(2)*(x(6)-x(3)*(x(6)-x(4)*(x(6)-x(5)*(x(6)-x(7).+y(7)*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(5)*(X-x(6)/(x(7)-x(1)*(x(7)-x(2)*(x(7)-x(3)*(x(7)-x(4)*(x(7)-x(5)*(x(7)-x(6);L6x=vpa(L

11、6)X=1.8;f1yp=subs(L6x)X=6.15;f2yp=subs(L6x)ezplot(L6x,0,7),title(六次lagrange插值函数),xlabel(X),ylabel(L6(X),hold onplot(x,y,r*)计算结果为：L6x = 0.111111111111111111111111*(X - 7.0)*(X - 3.0)*(X - 6.0)*(X - 5.0)*(X - 2.0)*(X - 4.0) - 0.001125*(X - 7.0)*(X - 3.0)*(X - 6.0)*(X - 5.0)*(X - 1.0)*(X - 4.0) - 0.000

12、5*(X - 7.0)*(X - 3.0)*(X - 6.0)*(X - 5.0)*(X - 1.0)*(X - 2.0) + 0.333333333333*(X - 7.0)*(X - 3.0)*(X - 6.0)*(X - 1.0)*(X - 2.0)*(X - 4.0) - 0.6666666666666666666666667*(X - 7.0)*(X - 3.0)*(X - 5.0)*(X - 1.0)*(X - 2.0)*(X - 4.0) + 0.66666666666666666666667*(X - 7.0)*(X - 6.0)*(X - 5.0)*(X - 1.0)*(X

13、- 2.0)*(X - 4.0) + 0.88888888888889*(X - 3.0)*(X - 6.0)*(X - 5.0)*(X - 1.0)*(X - 2.0)*(X - 4.0)f1yp = 0.1648f2yp =0.0013所以：f(1.8)= f1yp= 0.1648； f(6.15)= f2yp =0.0013图中红色点为（2）中给出的实验数据分布图。3、对于（2）给出的实验数据，选取节点x2=2,y2=0.135； x4=4,y4=0.018；x6=6,y2=0.002；进行lagrange插值，并画出其函数图，与六次lagrange进行比较。%三点Lagrange插值c

14、learclcsyms L2 X1xyp=1 2 3 4 5 6 7;yyp=0.368 0.135 0.050 0.018 0.007 0.002 0.001;x=2 4 6;y=0.135 0.018 0.002;L2=y(1)*(X1-x(2)*(X1-x(3)/(x(1)-x(2)*(x(1)-x(3)+y(2)*(X1-x(1)*(X1-x(3)/(x(2)-x(1)*(x(2)-x(3).+y(3)*(X1-x(1)*(X1-x(2)/(x(3)-x(1)*(x(3)-x(2);L2x=vpa(L2)h=ezplot(L2,0,7);set(h,Color,k);set(h,lin

15、estyle,-);hold onplot(xyp,yyp,r*);hold on%六次Lagrange插值syms L6 Xx=1 2 3 4 5 6 7;y=0.368 0.135 0.050 0.018 0.007 0.002 0.001;L6=y(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(5)*(X-x(6)*(X-x(7)/(x(1)-x(2)*(x(1)-x(3)*(x(1)-x(4)*(x(1)-x(5)*(x(1)-x(6)*(x(1)-x(7).+y(2)*(X-x(1)*(X-x(3)*(X-x(4)*(X-x(5)*(X-x(6)*(X-x(7)/(

16、x(2)-x(1)*(x(2)-x(3)*(x(2)-x(4)*(x(2)-x(5)*(x(2)-x(6)*(x(2)-x(7).+y(3)*(X-x(1)*(X-x(2)*(X-x(4)*(X-x(5)*(X-x(6)*(X-x(7)/(x(3)-x(1)*(x(3)-x(2)*(x(3)-x(4)*(x(3)-x(5)*(x(3)-x(6)*(x(3)-x(7).+y(4)*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(5)*(X-x(6)*(X-x(7)/(x(4)-x(1)*(x(4)-x(2)*(x(4)-x(3)*(x(4)-x(5)*(x(4)-x(6)*(x(4)-

17、x(7).+y(5)*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(6)*(X-x(7)/(x(5)-x(1)*(x(5)-x(2)*(x(5)-x(3)*(x(5)-x(4)*(x(5)-x(6)*(x(5)-x(7).+y(6)*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(5)*(X-x(7)/(x(6)-x(1)*(x(6)-x(2)*(x(6)-x(3)*(x(6)-x(4)*(x(6)-x(5)*(x(6)-x(7).+y(7)*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(5)*(X-x(6)/(

18、x(7)-x(1)*(x(7)-x(2)*(x(7)-x(3)*(x(7)-x(4)*(x(7)-x(5)*(x(7)-x(6);L6x=vpa(L6)ezplot(L6x,0,7),title(三点Lgrange插值和六次Lgrange插值的比较);%求函数值X1=1.8;X=1.8;k0=subs(L2x)k1=subs(L6x)X1=6.15;X=6.15;g0=subs(L2x)g1=subs(L6x)计算结果为：L2x =0.016875*(X1 - 6.0)*(X1 - 4.0) - 0.0045*(X1 - 6.0)*(X1 - 2.0) + 0.00025*(X1 - 2.0)

19、*(X1 - 4.0)L6x =0.111111111111111111111111*(X - 7.0)*(X - 3.0)*(X - 6.0)*(X - 5.0)*(X - 2.0)*(X - 4.0) - 0.001125*(X - 7.0)*(X - 3.0)*(X - 6.0)*(X - 5.0)*(X - 1.0)*(X - 4.0) - 0.0005*(X - 7.0)*(X - 3.0)*(X - 6.0)*(X - 5.0)*(X - 1.0)*(X - 2.0) + 0.333333333333*(X - 7.0)*(X - 3.0)*(X - 6.0)*(X - 1.0)*

20、(X - 2.0)*(X - 4.0) - 0.6666666666666666666666667*(X - 7.0)*(X - 3.0)*(X - 5.0)*(X - 1.0)*(X - 2.0)*(X - 4.0) + 0.66666666666666666666667*(X - 7.0)*(X - 6.0)*(X - 5.0)*(X - 1.0)*(X - 2.0)*(X - 4.0) + 0.88888888888889*(X - 3.0)*(X - 6.0)*(X - 5.0)*(X - 1.0)*(X - 2.0)*(X - 4.0) X1 =1.8000X = 1.8000k0

21、= 0.1523k1 = 0.1648X1 = 6.1500X = 6.1500g0 = 0.0049g1 =0.0013所以：三点Lagrange插值：f(1.8)= k0= 0.1523； f(6.15)= k1=0.0049六次Lagrange插值：f(1.8)= g0= 0.1648； f(6.15)= g1=0.0013对比给出的数据：f(1.8)= 0.165299；f(6.15)= 0.00213348通过比较，六次lagrange能较好的拟合所给数据，取的数据点越多，插值函数越精确。4、对于（1）给出的实验数据，编写MATLAB程序，计算五次牛顿插值多项式p5(x)，画出p5(

22、x)的函数图，计算p(0.596),p(0.99)值，并与lanrange插值法进行比较。%五次牛顿插值函数clearclcsyms p5 Xx=0.4 0.55 0.65 0.80 0.95 1.05;y=0.41075 0.57815 0.69675 0.90 1.00 1.25382;f0=(y(2)-y(1)/(x(2)-x(1);f01=(y(3)-y(2)/(x(3)-x(2);f02=(y(4)-y(3)/(x(4)-x(3);f03=(y(5)-y(4)/(x(5)-x(4);f04=(y(6)-y(5)/(x(6)-x(5);f1=(f01-f0)/(x(3)-x(1);f1

23、1=(f02-f01)/(x(4)-x(2);f12=(f03-f02)/(x(5)-x(3);f13=(f04-f03)/(x(6)-x(4);f2=(f11-f1)/(x(4)-x(1);f21=(f12-f11)/(x(5)-x(2);f22=(f13-f12)/(x(6)-x(3);f3=(f21-f2)/(x(5)-x(1);f31=(f22-f21)/(x(6)-x(2);f4=(f31-f3)/(x(6)-x(1);p5=y(1)+f0*(X-x(1)+f1*(X-x(1)*(X-x(2)+f2*(X-x(1)*(X-x(2)*(X-x(3)+f3*(X-x(1)*(X-x(2)

24、*(X-x(3)*(X-x(4)+f4*(X-x(1)*(X-x(2)*(X-x(3)*(X-x(4)*(X-x(5);p5x=vpa(p5)X=0.596,p1=subs(p5x),X=0.990,p2=subs(p5x),ezplot(p5x,0,2),title(),xlabel(X),ylabel(p5(X),hold on plot(x,y,r*)计算结果为：p5x =1.116*X + (X - 0.55)*(0.28*X - 0.112) + (0.99*X - 0.396)*(X - 0.65)*(X - 0.55) - 1.0*(15.3220202*X - 6.128808

25、081)*(X - 0.8)*(X - 0.65)*(X - 0.55) + (121.626355866355837311*X - 48.65348266)*(X - 0.95)*(X - 0.8)*(X - 0.65)*(X - 0.55) - 0.03565X = 0.5960p1 = 0.6257X = 0.9900p2 = 1.0542所以：p6(0.596)=p1=0.6257， p6(0.99)= p2=1.05425、对于（2）给出的实验数据，编写MATLAB程序，计算六次牛顿插值多项式p6(x)，画出p6(x)的函数图，计算p(1.8),p(6.15)值，并与lanrange

26、插值法进行比较。%六次牛顿插值syms p6 Xx=1 2 3 4 5 6 7;y=0.368 0.135 0.050 0.018 0.007 0.002 0.001;f0=(y(2)-y(1)/(x(2)-x(1);f01=(y(3)-y(2)/(x(3)-x(2);f02=(y(4)-y(3)/(x(4)-x(3);f03=(y(5)-y(4)/(x(5)-x(4);f04=(y(6)-y(5)/(x(6)-x(5);f05=(y(7)-y(6)/(x(7)-x(6);f1=(f01-f0)/(x(3)-x(1);f11=(f02-f01)/(x(4)-x(2);f12=(f03-f02)

27、/(x(5)-x(3);f13=(f04-f03)/(x(6)-x(4);f14=(f05-f04)/(x(7)-x(5);f2=(f11-f1)/(x(4)-x(1);f21=(f12-f11)/(x(5)-x(2);f22=(f13-f12)/(x(6)-x(3);f23=(f14-f12)/(x(7)-x(4);f3=(f21-f2)/(x(5)-x(1);f31=(f22-f21)/(x(6)-x(2);f32=(f23-f22)/(x(7)-x(3);f4=(f31-f3)/(x(6)-x(1);f41=(f32-f31)/(x(7)-x(2);f5=(f41-f4)/(x(7)-x(1);p6=y(1)+f0*(X-x(1)+f1*(X-x(1)*(X-x(2)+f