数值分析实习报告插值法.doc

x大学计算机与信息学院（数学类课程）课程实习报告课程名称：数值分析课程实习实习题目：插值法姓 名：系：信息与计算科学专 业：信息与计算科学年 级：2009学 号：指导教师：职 称：讲师2012 年 05 月 1 日课程实习报告评分表序号项目评分权值评定成绩1算法：正确，步骤合理0.32程序：画流程图，用语言或数学软件编写0.23质量：分析处理科学；文字通顺；计算及测试结果准确；0.44工作量、工作态度：按期完成规定的任务，工作量饱满；工作努力，遵守纪律0.1合 计1指导教师签名： 年 月 日目录1.实习的目的和任务12.实习要求13.实习地点14.主要仪器设备15.实习内容15.1问题描述15.2Lagrange插值法15.3逐次线性插值法65.4牛顿插值法85.5等距节点插值法105.6 Hermite插值法155.7分段低次插值法186问题讨论与分析197结束语20参考文献20数值分析课程实习1. 实习的目的和任务目的：通过课程实习能够应用MATLAB软来计算函数的插值；了解函数插值方法。任务：通过具体的问题，利用MATLAB软件来计算问题的结果，分析问题的结论。2. 实习要求能够从案例的自然语言描述中，抽象出其中的数学模型；能够熟练应用所学的函数插值方法；能够熟练使用MATLAB软件；对函数的插值方法有所认识，包括对不同插值法有所认识。3. 实习地点学生宿舍4. 主要仪器设备计算机Microsoft Windows 7Matlab R2010b5. 实习内容5.1 问题描述给出函数f(x)=lnx的数值表（见表5.1）用Lagrange插值法，逐次线性插值法，牛顿插值法，等距节点插值法，Hermite插值法，分段低次插值法分别求该函数0.4，0.8上的插值多项式并计算ln0.53的近似值。 表5.1x0.40.50.60.70.8lnx-0.916291-0.693147-0.510826-0.357765-0.2231445.2 Lagrange插值法拉格朗日插值法是基于基函数的插值方法，插值多项式可表示为： 其中称为i次基函数5.2.1线性插值和抛物插值线性插值选取x=0.4 x=0.5 在matlab命令窗口中输入：x=0.4 0.5;y=-0.916291 -0.693147;L0=poly(x(2)/(x(1)-x(2);L1=poly(x(1)/(x(2)-x(1);z=y(1)*L0+y(2)*L1;a=poly2sym(z) 回车结果为 a=(27893*x)/12500-1.8088670b=polyval(z,0.53) 回车结果为 b=-0.6262038在matlab命令窗口中输入：x1=0.3:0.001:0.6;y1=(27893*x1)/12500-1.8088670;x2=0.4 0.5 0.53;y2=-0.916291 -0.693147 -0.6262038;plot(x1,y1,'b-',x2,y2,'o');title('图 5-2-1 线性插值')回车后得到如图所示结果抛物插值选取x=0.4 x=0.5 x=0.6在matlab命令窗口中输入：x=0.4 0.5 0.6;y=-0.916291 -0.693147 -0.510826;L0=poly(x(2) x(3)/(x(1)-x(2)*(x(1)-x(3);L1=poly(x(1) x(3)/(x(2)-x(1)*(x(2)-x(3);L2=poly(x(1) x(2)/(x(3)-x(1)*(x(3)-x(2);z=y(1)*L0+y(2)*L1+y(3)*L2;a=poly2sym(z)回车结果为 a=-(40823*x2)/20000+4.068475*x-2.217097b=polyval(z,0.53) 回车结果为 b=-0.6341643在matlab命令窗口中输入：x=0.35:0.001:0.65;y=-(40823*x.2)/20000+4.068475*x-2.217097;x1=0.4 0.5 0.53 0.6;y1=-0.916291 -0.693147 -0.6262038 -0.510826;plot(x,y,'g+',x1,y1,'o');title('图 5-2-2 抛物插值')回车后得到如图所示结果5.2.2拉格朗日插值法的matlab实现建立M文件Language.m:function f=Language(x,y,x0)%求已知数据点的拉格朗日插值多项式%已知数据点的x坐标向量: x%已知数据点的y坐标向量: y%插值点的x坐标: x0%求得的拉格朗日插值多项式或在x0处的插值: f syms t;if(length(x)=length(y) n=length(x);else disp('x和y的维度不相等！'); return;%检错endf=0.0;for(i=1:n) l=y(i); for(j=1:i-1) l=l*(t-x(j)/(x(i)-x(j); end for(j=i+1:n) l=l*(t-x(j)/(x(i)-x(j); end; f=f+l; simplify(f);%化简 if(i=n) if(nargin=3) f=subs(f,'t',x0);%计算插值点的函数值 else f=collect(f); f=vpa(f,8);%将插值多项式的系数化成8位精度的小数 end endend在matlab工作窗口输入如下程序：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f=Lagrange(x,y)运行的结果为：f =- 0.30958333*t4 + 2.60825*t3 - 5.4860542*t2 + 5.6920925*t - 2.474362再次输入f=Lagrange(x,y,0.53)结果为 f= -0.6347048在matlab命令窗口中输入：t=0.35:0.001:0.85;y=- 0.30958333*t.4+ 2.60825*t.3 - 5.4860542*t.2 + 5.6920925*t - 2.474362;x1=0.4 0.5 0.53 0.6 0.7 0.8;y1=-0.916291 -0.693147 -0.6347048 -0.510826 -0.357765 -0.223144;plot(t,y,'r*',x1,y1,'o');title('图 5-2-3 Lagrange插值')回车后得到如图所示结果5.3逐次线性插值法5.3.1逐次线性插值法基本原理定义为，一阶定义均差为 一般地，函数f的k阶均差定义为：5.3.2逐次线性插值法的matlab实现建立M文件Aitken.mfunction f = Aitken(x,y,x0)%求已知数据点的Aitken逐次线性插值多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%为插值点的x坐标:x0%求得的Aitken插值多项式或x0处的插值:fsyms t;if(length(x) = length(y) n = length(x);else disp('x和y的维数不相等！'); return; %检错endy1(1:n)=t; for(i=1:n-1) for(j=i+1:n) y1(j) = y(j)*(t-x(i)/(x(j)-x(i)+y(i)*(t-x(j)/(x(i)-x(j); end y=y1; simplify(y1); %化简endif(nargin = 3) f = subs(y1(n),'t',x0); %计算插值点的函数值else simplify(y1(n); f = collect(y1(n); %将插值多项式展开 f = vpa(f,8); %将插值多项式的系数化成8位精度的小数end在matlab工作窗口输入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= Aitken(x,y)运行的结果为：f =- 0.30958333*t4 + 2.60825*t3 - 5.4860542*t2 + 5.6920925*t - 2.474362再次输入f= Aitken(x,y,0.53)结果为f= -0.6347048在matlab命令窗口中输入：t=0.35:0.001:0.85;y=- 0.30958333*t.4 + 2.60825*t.3 - 5.4860542*t.2 + 5.6920925*t - 2.474362x1=0.4 0.5 0.53 0.6 0.7 0.8;y1=-0.916291 -0.693147 -0.6347048 -0.510826 -0.357765 -0.223144;plot(t,y,'cd',x1,y1,'o');title('图 5-3 Aitken插值')回车后得到如图所示结果5.4牛顿插值法5.4.1牛顿插值法的基本原理零阶均差定义为，一阶定义均差为 一般地，函数f的k阶均差定义为：利用均差的牛顿插值法多项式为5.4.2牛顿插值法的matlab实现建立M文件Newton.m：function f = Newton(x,y,x0)%求已知数据点的均差形式牛顿插值多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%为插值点的x坐标:x0%求得的均差形式牛顿插值多项式或x0处的插值:fsyms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp('x和y的维数不相等！'); return; %检错endf = y(1);y1 = 0;l = 1; for(i=1:n-1) for(j=i+1:n) y1(j) = (y(j)-y(i)/(x(j)-x(i); end c(i) = y1(i+1); l = l*(t-x(i); f = f + c(i)*l; simplify(f); %化简 y = y1; if(i=n-1) if(nargin = 3) f = subs(f,'t',x0); %计算插值点的函数值 else f = collect(f); %将插值多项式展开 f = vpa(f,8); %将插值多项式的系数化成8位精度的小数 end endend在matlab工作窗口输入如下程序：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= Newton (x,y)运行的结果为f = - 0.30958333*t4 + 2.60825*t3 - 5.4860542*t2 + 5.6920925*t - 2.474362再次输入 f= Newton (x,y,0.53)结果为 f= -0.6348783在matlab命令窗口中输入：t=0.35:0.001:0.85;y=- 0.30958333*t.4 + 2.60825*t.3 - 5.4860542*t.2 + 5.6920925*t - 2.474362x1=0.4 0.5 0.53 0.6 0.7 0.8;y1=-0.916291 -0.693147 -0.6348783 -0.510826 -0.357765 -0.223144;plot(t,y,'m- ',x1,y1,'o');title('图 5-4 Newton插值')回车后得到如图所示结果5.5等距节点插值法5.5.1向前牛顿插值法的基本原理向前牛顿插值多项式可表示如下：其中h叫做步长，且t的取值范围为 。5.5.2向前牛顿插值法的matlab实现建立M文件Newtonforward.m：function f = Newtonforward(x,y,x0)%求已知数据点的向前差分牛顿插值多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%为插值点的x坐标:x0%求得的向前差分牛顿插值多项式或x0处的插值:fsyms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp('x和y的维数不相等！'); return; %检错endf = y(1);y1 = 0;xx =linspace(x(1),x(n),(x(2)-x(1);if(xx = x) disp('节点之间不是等距的！'); return;endfor(i=1:n-1) for(j=1:n-i) y1(j) = y(j+1)-y(j); end c(i) = y1(1); l = t; for(k=1:i-1) l = l*(t-k); end; f = f + c(i)*l/factorial(i); simplify(f); %化简 y = y1; if(i=n-1) if(nargin = 3) f = subs(f,'t',(x0-x(1)/(x(2)-x(1); else f = collect(f); %将插值多项式展开 f = vpa(f,8); %将插值多项式的系数化成8位精度的小数 end endend在matlab工作窗口输入如下程序：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= Newtonforward (x,y)运行的结果为f = - 0.000030958333*t4 + 0.0021129167*t3 - 0.026533542*t2 + 0.24759558*t - 0.916291再次输入 f= Newtonforward (x,y,0.53)结果为 f= -0.6348783在matlab命令窗口中输入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= Newtonforward(x,y)t=0.4:0.01:0.8;yy1=Newtonforward(x,y,t);plot(x,y,'ro',t,yy1, ':< ')title('图 5-5-1 向前牛顿插值法')回车后得到如图所示结果5.5.3向后牛顿插值法的基本原理向后牛顿插值多项式可表示如下：其中h叫做步长，且t的取值范围为。5.5.4向后牛顿插值法的matlab实现建立M文件Newtonback.m：function f = Newtonback(x,y,x0)%求已知数据点的向后差分牛顿插值多项式%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%为插值点的x坐标:x0%求得的向前差分牛顿插值多项式或x0处的插值:fsyms t;if(length(x) = length(y) n = length(x); c(1:n) = 0.0;else disp('x和y的维数不相等！'); return;endf = y(n);y1 = 0;xx =linspace(x(1),x(n),(x(2)-x(1);if(xx = x) disp('节点之间不是等距的！'); return; %检错endfor(i=1:n-1) for(j=i+1:n) y1(j) = y(j)-y(j-1); end c(i) = y1(n); l = t; for(k=1:i-1) l = l*(t+k); end; f = f + c(i)*l/factorial(i); simplify(f); %化简 y = y1;if(i=n-1) if(nargin = 3) f = subs(f,'t',(x0-x(n)/(x(2)-x(1); else f = collect(f); %将插值多项式展开 f = vpa(f,8); %将插值多项式的系数化成8位精度的小数 end endend在matlab工作窗口输入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= Newtonback (x,y)运行的结果为f = - 0.000030958333*t4 + 0.0016175833*t3 - 0.0041505417*t2 + 0.12882192*t - 0.223144再次输入 f= Newtonback (x,y,0.53)结果为 f= -0.6347048在matlab命令窗口中输入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;f= Newtonback (x,y)t=0.4:0.01:0.8;yy1=Newtonback (x,y,t);plot(x,y,'mo',t,yy1, ':x ')title('图 5-5-2 向后牛顿插值法')回车后得到如图所示结果5.6 Hermite插值法5.6.1 Hermite插值法的基本原理Hermite插值满足在节点上等于给定函数值，而且在节点上的导数值也等于给定导数值。对于高阶导数情况，Hermite插值多项式比较复杂，在实际中，常遇到的是函数值与一阶导数给定的情况。在此情况下，n个节点x1，x2，xn的Hermite插值多项式的表达式如下：其中，。5.6.2 Hermite插值法的matlab实现建立M文件Hermite.m：function f = Hermite(x,y,y_1,x0)%已知数据点的x 坐标向量:x%已知数据点的y 坐标向量:y%已知数据点的导数向量:y_1%求得的Hermite插值多项式或x0处的插值:fsyms t;f = 0.0;if(length(x) = length(y) if(length(y) = length(y_1) n = length(x); else disp('y和y的导数的维数不相等！'); return; endelse disp('x和y的维数不相等！'); return; endfor i=1:n h = 1.0; a = 0.0; for j=1:n if( j = i) h = h*(t-x(j)2/(x(i)-x(j)2); a = a + 1/(x(i)-x(j); end end f = f + h*(x(i)-t)*(2*a*y(i)-y_1(i)+y(i); if(i=n) if(nargin = 4) f = subs(f,'t',x0); else f = vpa(f,8); %将插值多项式的系数化成8位精度的小数 end endend在matlab工作窗口输入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;y_1=5/2 2 5/3 10/7 5/4;f= Hermite (x,y, y_1)运行的结果为f = 6250000.0*(1.6666667*t - 1.510826)*(t - 0.4)2*(t - 0.8)2*(t - 0.7)2*(t - 0.5)2 + 173611.11*(10.547667*t - 8.6612773)*(t - 0.4)2*(t - 0.7)2*(t - 0.6)2*(t - 0.5)2 + 2777777.8*(7.3913214*t - 5.53169)*(t - 0.4)2*(t - 0.8)2*(t - 0.6)2*(t - 0.5)2 - 173611.11*(35.678792*t - 13.355226)*(t - 0.8)2*(t - 0.7)2*(t - 0.6)2*(t - 0.5)2 - 2777777.8*(9.55245*t - 4.083078)*(t - 0.4)2*(t - 0.8)2*(t - 0.7)2*(t - 0.6)2再次输入 f= Hermite (x,y, y_1,0.53)结果为 f= -0.6349413在matlab命令窗口中输入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;y_1=5/2 2 5/3 10/7 5/4;f= Hermite (x,y,y_1)t=0.4:0.01:0.8;yy1=Hermite(x,y,y_1,t);plot(x,y,'o',t,yy1)title('图 5-6 Hermite插值')回车后得到如图所示结果5.7分段低次插值法5.7.1分段线性插值在matlab命令窗口中输入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;xi=0.53;yi=interp1(x,y,xi,'linear')结果为 yi =-0.63845070在matlab命令窗口中输入：x=0.4 0.5 0.6 0.7 0.8;y=-0.916291 -0.693147 -0.510826 -0.357765 -0.223144;xi=0.4:0.1:0.8;yi=interp1(x,y,xi,'linear');plot(x,y,'ro',x,yi);title('图 5-7 分段线性插值法')回车后得到如图所示结果6 问题讨论与分析ln0.53=-0.634878272435970 表5-8 不同插值方法的比较插值方法近似值误差线性插值-0.6262038000000000.008674472抛物插值-0.6341642850000027.1399e-004拉格朗日插值-0.6347047692624991.7350e-004逐次线性插值-0.6347047692625001.7350e-004牛顿插值-0.6347047692625001.7350e-004向前牛顿插值-0.6347047692625001.7350e-004向后牛顿插值-0.6347047692625001.7350e-004Hermite插值-0.6349412808205916.3008e-005分段线性插值-0.6384507000000000.003572427通过上表，可以看出Hermite插值效果最好，精度最高，线性插值、抛物插值、分段线性插值效果最差，误差较大。另外拉格朗日插值、逐次线性插值、牛顿插值、向前牛顿插值与向后牛顿插值效果基本相同，精度较高。其中牛顿插值法是一种利用均差构造插值多项式的方法，n次牛顿插值多项式与n次拉格朗日插值多项式是恒等关系，只是表现的形式不同而已。Hermite插值是一种带导数信息的插值方法，插值效果好，插值条件要求较高。而关于等距节点插值法，当插值点x接近数据表头时一般用向前牛顿插值法，当插值点x接近数据表头时一般用向后牛顿插值法。7结束语通过本次数值分析的课程实习，使我更加深入的了解了函数的各种插值方法，并利用matlab实现的过程。通过每天的学习与查阅资料，让我知晓到了很多书本上没有的知识。我也懂得应该如何把在书本上学到的理论知识，利用matlab去实现。但在这段时间里，我也遇到了很多难题，也许在matlab中一个小符号，小标点就让我几欲发狂，觉得自己的动手能力还是有待提高。另外，耐心的去做一件事，我相信总会收获很多。参考文献1姜健飞、胡良剑、唐俭.数值分析及其MATLAB实验上海：科学出版社，20042张德丰.matlab数值计算方法.北京:机械工业出版社，20103张圣勤.MATLAB7.0实用教程.北京:机械工业出版社，20044刘卫国.MATLAB程序设计与应用第二版.北京:高等教育出版社，2006