数值计算课程设计报告(插值法)

1、数 值 计 算 方 法课 程 设 计 报 告课程设计名称： 数值计算方法 课程设计题目： 插值算法 年 级 专 业： 信计1302班 组员姓名学号： 高育坤 1309064043 王冬妮 1309064044 韩 建 1309064046 李 婧 1309064047 指 导 教 师： 刘丽华 完 成 时 间： 2015年6月17日 插值算法一、问题提出插值法是实用的数值方法，是函数逼近的重要方法。在生产和科学实验中，自变量x与因变量y的函数y = f(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。如何根

2、据观测点的值，构造一个比较简单的函数y=(x)，使函数在观测点的值等于已知的数值或导数值，进而用简单函数y=(x)在点x处的值来估计未知函数y=f(x)在x点的值。寻找这样的函数(x)，办法是很多的。(x)可以是一个代数多项式，或是三角多项式，也可以是有理分式；(x)可以是任意光滑（任意阶导数连续）的函数或是分段函数；函数类的不同，自然地有不同的逼近效果。二、背景分析 在许多实际问题及科学研究中，因素之间往往存在着函数关系，然而，这种关系经常很难有明显的解析表达，通常只是由观察与测试得到一些离散数值。有时，即使给出了解析表达式，却由于表达式过于复杂，不仅使用不便，而且不易于进行计算与理论分析。

3、解决这类问题的方法有两种：一种是插值法插值法，另一种是一拟合法。插值法是一种古老的数学方法，它来自生产实践，早在一千多年前，我国科学家在研究历法上就应用了线性插值与二次插值，但它的基本理论却是在微积分产生之后才逐渐完善的，其应用也日益增多，特别是在计算机软件中，许多库函数，如 ,cos,sin ex 等的计算实际上归结于它的逼近函数的计算。逼近函数一般为只含有算术运算的简单函数，如多项式、有理分式（即多项式的商）。在工程实际问题当中，我们也经常会碰到诸如此类的函数值计算问题。被计算的函数有时不容易直接计算，如表达式过于复杂或者只能通过某种手段获取该函数在某些点处的函数值信息或者导数值信息等。因

4、此，我们希望能用一个“简单函数” 逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。这种方法就叫插值逼近或者插值法。插值法要求给出函数f（x）的一个函数表，然后选定一种简单的函数形式，比如多项式、分段线性函数及三角多项式等，通过已知的函数表来确定一个简单的函数作为f（x）的近似，概括地说，就是用简单函数为离散数组建立连续模型。三、基本算法思想与实现·····已知个数据节点:构造一个(相对简单)函数(称为插值函数),通过全部结点即（j0,1,n）再用计算插值，即数学上插值方法非常多，这里介绍几种常用方法：1·插值函数插

5、值函数的基本思想:将待求的次插值多项式写成另一种表达方,式再利用插值条件 确定出插值基函由基函数条件,确定多项式系数，进而可得插值函数.（1）已知，求满足条件的插值函数。 由题可知表示过两点的直线，这个问题是我们所熟悉的，它的解可表为下列对称式 此类一次插值称为线性插值，若令（由此可得：）则有这里的可以看作是满足条件的插值多项式，这两个特殊的插值多项式称作上述问题的插值基函数。（2）求过三点的插值函数。 为了得到插值多项式先解决一个特殊的二次插值问题。 求作二次式，使满足 （2-1）这个问题是容易求解的，由式（2-1）的后两个条件知是的两个零点，因而 。再用条件确定系数c .结果得 ：类似可以

6、分别构造出满足条件的插值多项式；其表达式分别为 ，这样构造出的称作问题（2）的插值基函数。设取已知数据作为组合系数，将插值基函数组合得 验证可知，这样构造的满足已知条件，因而它就是问题（2）的解。 (3)推广到一般：已知函数在n+1个不同点上的函数值分别为求一个次数不超过n的多项式,使其满足:即个不同的点可以决定的一个次多项式。过个不同的点分别决定个次插值基函数。每个插值基多项式满足： a. 是次多项式;   b. ,而在其它个点 由于故有因子: 因其已经是n次多项式，故而仅相差一个常数因子。令:     

7、;    由,可以定出, 进而得到：  次拉格朗日型插值多项式是个次插值基本多项式的线性组合,相应的组合系数是。即:   从而是一个次数不超过n的多项式,且满足2·插值函数的构造插值法的基本思想：已知节点处的函数值或一元函数代数方程，将待求的n次插值多项式改写为具有承袭性的形式，然后根据插值条件或选取初值以求得待定系数，进而求得所要的插值函数。实践中的许多问题归结为求一元代数方程的根，如果是线性函数，则它的求根较容易；对非线性方程，只有不高于4次的代数方程有求根公式，经常需求出高于4次 

8、;的满足一定精度要求的近似解。   法的简述设是的一个近似根，把在处泰勒展开 若取前两项来近似代替，则的近似线性方程设0，设其根为，则的计算公式为=- （k=0,1,2.）这即为牛顿法，上式为牛顿迭代公式，其迭代函数为 我们知道，牛顿法是解非线性方程最著名和最有效的方法之一，在单根附近它比一般的迭代格式有较快的收速度，但也要注意它也有缺点：首先，它对迭代初值选取要求较严，初值选取不好，可能导致吧收敛；其次，它每迭代一次要计算的值，这势必增加可计算量。为回避该问题，常用一个固定的迭代若干步后再求。这就是下面要讲的简化牛顿法的基本思想。简化牛顿法和下山牛顿法简化牛顿法的公式为

9、(3-1)迭代函数 若。即 在根附近成立。则迭代法（3-1）局部收敛。此法显然化简了计算量。牛顿下山法牛顿法的收敛依赖于初值的选取，若偏离较远，则牛顿法可能发散。为防止迭代发散，我们对迭代过程在附加一项条件，即具有单调性： （3-2）保证函数值稳定下降，然后结合牛顿法加快收敛速度，即可达目的。将牛顿法的计算结果 （3-3）于前一步的近似值适当加权平均作为新的改进值 （3-4）其中称（）为下山因子，即为 （3-5）称为牛顿下山法。选择下山因子时，从开始逐次将减半进行试算。直到满足条件（3-2）为止。 3·插值法已知函数在给定个互异的节点，.上的函数值和导数值，求一个次多项式满足插值条件

10、 （）=， . k=0,1,2.n插值基本原理通常如上条件的Hermite型插值是通过构造相应的插值基函数来完成的，为方便起见以n=1为例，说明传统的求解方法，设给定的，和相应的函数值，及微商值，构造插值函数。由构造函数的办法可知：对应于和点函数值的插值函数分别为 及而对应的和点导数值的插值基函数分别为和 ，因此所要求的插值函数 （2-1）由上可发现构造插值基函数比较复杂，尤其对具有高阶导数插值条件的情况，以下将基于newton插值方法提出构造上述条件的简单格式。此时传统方法可视为这里的特例。四、具体应用实例分析1 已知,用线性插值法求的近似值.解:Matlab中有直接进行线性插值计算的命令i

11、nterp1,直接使用interp1命令即可.>> x=4 9;>> y=2 3;>> f=interp1(x,y,7,'linear')%选项使用线性插值f = 2.6000故插值计算结果为2 设,给出数据如下, 用Lagrange插值法求的近似值.0.40.50.70.8-0.916291-0.693147-0.356675-0.223144解:求解过程描述如下: format long;%输入初始数据x0=0.4 0.5 0.7 0.8;y0=-0.916291 -0.693147 -0.356675 -0.223144;x=0.6;%

12、插值点n=length(x0);s=0;%进入迭代计算过程for j=0:(n-1) t=1; for i=0:(n-1) if i=j t=t*(x-x0(i+1)/(x0(j+1)-x0(i+1); end end s=s+t*y0(j+1);ends %显示输出结果format short;程序运行结果如下:s = -0.509975500000000因此利用插值的计算结果为.3 设有如下数据, 利用插值法求的近似值.0.400.550.650.800.901.050.410750.578150.696750.888111.026521.25382解:求解程序如下clc;format l

13、ong;%显示15位x0=0.40 0.55 0.65 0.80 0.90 1.05;%x的值y0=0.41075 0.57815 0.69675 0.88811 1.02652 1.25382;%y的值x=0.596;%插值点n=max(size(x0);y=y0(1);%迭代初始值disp(y);s=1;dx=y0;for i=1:n-1%构造差商表 dx0=dx; for j=1:n-i dx(j)=(dx0(j+1)-dx0(j)/(x0(i+j)-x0(j); end df=dx(1); s=s*(x-x0(i); y=y+s*df;%计算 disp(y);end运行上述程序结果如下

14、:0.410750000000000 0.629486000000000 0.632010480000000 0.631914405504000 0.631917508079616 0.631917499231745因此插值结果为 4 给出的数据见下表,用Hermite插值多项式求的近似值,并估计其误差.0.400.500.700.80-0.916291-0.693147-0.356675-0.2231442.502.001.431.25解：先建立实现插值的M文件函数，源程序如下：function y=hermite(x0,y0,dy,x) %hermite.m%Hermite插值计算%x0为

15、输入节点的向量；y0为y的值向量，%dy为相应节点一阶倒数的函数值的向量，x为所要求的插值节点. n=length(x0); m=length(x); for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j=i h=h*(x(k)-x0(j)/(x0(i)-x0(j)2; a=1/(x0(i)-x0(j)+a; end end yy=yy+h*(x0(i)-x(k)*(2*a*y0(i)-dy(i)+y0(i); end y(k)=yy;end在Matlab命令窗口中，进行以下步骤：输入数据：>>x0=0.4 0.5 0.7

16、 0.8 x0 =0.40000 0.50000 0.70000 0.80000 >>y0=log(x0) y0 = -0.916290731874155 -0.693147180559945 -0.356674943938732 -0.223143551314210 >>dy=1./x0 %一阶导数的值dy = 2.500000000000000 2.000000000000000 1.428571428571429 1.250000000000000 >>x=0.6 %插值点x = 0.600000000000000 >>H=hermite(

17、x0,y0,dy,x) H = -0.510824121030042 由插值计算结果有，准确值为：>>log(0.6) ans = -0.510825623765991可见精度是比较高的.5 给出的数据见下表,用三次Hermite插值多项式求的近似值,精确到6位小数,并估计其误差.(精确值.)0.5000000.7071070.8660250.8660250.7071070.500000解:对于分段三次插值,Matlab里已有现成命令可用,这个命令就是interp1(),只要在选项中选用'chip'(或cubic)即可,它是用分段三次多项式插值曲线依次连接相邻样本点

18、的意思,整体上具有函数及其一阶导数连续性.我们先用这个命令来解本题.输入数据:>>x0=pi/6,pi/4,pi/3 x0 = 0.523598775598299 0.785398163397448 1.047197551196598 >>y0=sin(x0) y0 = 0.500000000000000 0.707106781186547 0.866025403784439 >>x=2*pi/9 x = 0.698131700797732 >>y=interp1(x0,y0,x,'pchip') y = 0.6438953599

19、99574 >>y1=interp1(x0,y0,x,'cubic') y1 = 0.643895359999574 可见集成命令的插值结果为: ,精确值为:>>sin(x) ans = 0.642787609686539由于题目要求用两点三次插值公式和三点三次公式求解,直接按课本公式代入运算即可.两点三次插值结果(选用点和的数据):>>x0=pi/6,pi/4,pi/3 x0 = 0.523598775598299 0.785398163397448 1.047197551196598 >>y0=sin(x0) y0 = 0.5

20、00000000000000 0.707106781186547 0.866025403784439 >>dy=cos(x0) dy = 0.866025403784439 0.707106781186548 0.500000000000000 >>x=2*pi/9 x = 0.698131700797732 >>H2=(1+2*(x-x0(1)/(x0(2)-x0(1)*(x-x0(2)/(x0(1)-x0(2)2*y0(1)+. (1+2*(x-x0(2)/(x0(1)-x0(2)*(x-x0(1)/(x0(2)-x0(1)2*y0(2)+. (x-x0

21、(1)*(x-x0(2)/(x0(1)-x0(2)2*dy(1)+. (x-x0(2)*(x-x0(1)/(x0(2)-x0(1)2*dy(2) H2 = 0.642781665744313三点三次插值结果:>>h0=(x-x0(2)/(x0(1)-x0(2)2*(x-x0(3)/(x0(1)-x0(3);>>h1=(1-(x-x0(2)*(1/(x0(2)-x0(1)+1/(x0(2)-x0(3)*(x-x0(1)*(x-x0(3)/(x0(2)-x0(1)*(x0(2)-x0(3);>>h2=(x-x0(1)/(x0(3)-x0(1)*(x-x0(2)/

22、(x0(3)-x0(2)2;>>dh=(x-x0(2)*(x-x0(1)*(x-x0(3)/(x0(2)-x0(1)*(x0(2)-x0(3);>>H3=h0*y0(1)+h1*y0(2)+h2*y0(3)+dh*dy(2) H3 = 0.6428010064027956 根据函数的数据(如下表所示),分别用两点一次插值、带导数的二次插值、两点三次插值计算的近似值，比较其精度.149161234解:先输入数据:>>x0=1 4 9 16 x0 = 1 4 9 16 >>y0=sqrt(x0) y0 = 1 2 3 4 >>dy=1./

23、2*1./sqrt(x0) dy = 0.500000000000000 0.250000000000000 0.166666666666667 0.125000000000000>>x=5 %插值点x = 5(1)两点一次插值,选取最靠近5的点和两组数据计算:%两点一次插值%找到两组最靠近插值点的数据>>for i=1:length(x0) if x0(i)<x & x0(i+1)>x break; endend>>y1=(x-x0(i+1)/(x0(i)-x0(i+1)*y0(i)+(x-x0(i)/(x0(i+1)-x0(i)*y0

24、(i+1) %计算插值点的插值结果y1 = 2.200000000000000(2)带导数的二次插值 因为,因此利用(1)中找到的两组数据即可,不同之处在于要附加利用上导数的值.>>y2=(1-(x-x0(i)/(x0(i)-x0(i+1)*(x-x0(i+1)/(x0(i)-x0(i+1)*y0(i)+power(x-x0(i)/(x0(i+1)-x0(i),2)*y0(i+1)+(x-x0(i)*(x-x0(i+1)/(x0(i)-x0(i+1)*dy(i) y2 = 2.240000000000000(3)两点三次插值,直接代入计算即可,与(2)不同在于多了一个导数条件.&g

25、t;>y3=(1-2*(x-x0(i)/(x0(i)-x0(i+1)*power(x-x0(i+1)/(x0(i)-x0(i+1),2)*y0(i)+. (1-2*(x-x0(i+1)/(x0(i+1)-x0(i)*power(x-x0(i)/(x0(i+1)-x0(i),2)*y0(i+1)+. (x-x0(i)*power(x-x0(i+1)/(x0(i)-x0(i+1),2)*dy(i)+. (x-x0(i+1)*power(x-x0(i)/(x0(i+1)-x0(i),2)*dy(i+1) y3 = 2.237333333333333此外,由于Matlab中已集成有分段三次插值函

26、数,直接利用该命令也可解答此题.>> y=interp1(4 9,2 3,5,'pchip')y = 2.200000000000000五、设计总结 总结：在条件有限情况下，构造固定的阶数的插值多项式可能会是一种简单的方案，当要反复计算逼近值时，最好用牛顿插值多项式；对于表格数据的常规插值，最好使用分段线性插值；如果插值总体平滑很重要，应该考虑运用三次样条插值或三次Hermite插值，同时表格数据构成函数的导数不存在时，最好使用三次样条插值.六、参考文献1 计算方法引论 徐萃薇,孙绳武 高等教育出版社 2001.4 2 数值分析 (美)David Kincaid &

27、amp; Ward Cheney 机械工业出版社 2005.9 (译)王国荣,俞耀明,徐兆亮 3 数值计算原理 李庆扬,关冶,白峰杉 清华大学出版社 2000.9 4 数值计算方法 郑慧娆,陈绍林,莫忠良,黄象鼎 武汉大学出版社 2002.10 5 现代数值分析 蒋耀林 国防工业出版社 2004.9 6 现代数值数学和计算 同济大学计算数学教研室 同济大学出版社 2004.7 7 数值计算方法及其应用 朱长青 科学出版社 2006.1 8 数值计算方法 曾金平 湖南大学出版社 2004.8 9 数值计算基础 沈剑华 同济大学出版社 2004.5 10 数值分析算法描述 徐士良 机械工业出版社

28、2003.4 11 数值分析 李庆扬,王能超,易大义 清华大学出版社&施普林格出版社 2001.8七、心得体会（1）韩建：通过此次课程设计我们可以知道计算机在现代生活中的应用已经如此普及，尤其是在数学计算当中，Matlab软件更是发挥了不可替代的作用Matlab以其强大的功能，方便了当今数值计算，数学教程，及工程计算等众多领域拉格朗日插值的优点：它的形式是对称的，这样很容易编程上机实现。它在理论上十分重要。 牛顿插值的优点：在计算插值多项式及求解函数近似值都比较方便且计算量相对较小。从公式中可以看出：每增加一个节点，插值多项式只增加一项，因此便于递推运算，所以其具有灵活增加节点的优点。

29、（2）高育坤： 深入了解matlab运行环境和操作环境，初步学会调试程序，运用绘图命令制作函数图象；了解常见几种插值法，以及数值分析的解决方案；懂得如何运用已有的知识更进一步了解未知的问题；独立解决和思考问题的能力有了一定的提高。（3）李婧：通过自己动手作实验学习如何用插值方法解决实际问题，提高探索和解决问题的能力。通过撰写实验报告，促使自己提炼思想，按逻辑顺序进行整理，并以他人能领会的方式表达自己思想形成的过程和理由。提高了写作、文字处理、排版等方面的能力。（4）王冬妮：各种插值法都有自己的利与弊，拉格朗日插值法运算过程相对复杂，但当和导数结合起来，组成抛物插值的时候，精度就可以提高很多。牛

30、顿插值法、拉格朗日插值法等线性插值法只能适合在已知点不多的情况下使用，当已知的坐标点很多时候应该将区间分成小段进行分段线性插值或者分段抛物插值。采用组合的思想是数值分析常用的思想和技巧之一，单个的方法得到的结果虽然不是很理想，但将多个方法按照某种方式结合在一起就能改进实验方法，我们应该触类旁通，在以后的学习中学会使用这种思想。 八、附录1分段插值的MATLAB实现 在MATLAB编辑窗口中输入 x=0:10; y=1./(1+x.2); xi=0:1:10; yi=interp1(x,y,xi,'linear'); t=0:0.1:10; z=1./(1+t.2); plot(

31、x,y,'p',xi,yi,t,z) 2 Hermite插值的MATLAB实现在MATLAB编辑窗口中输入 x=1,2; y=2,3; y1=0,-1; f=Hermite(x,y,y1,1.5) f= 2.625运行程序时调用的Hermite函数如下： function f=Hermite(x,y,y_1,x0) syms t; f=0.0; if(length(x)=length(y) if(length(y)=length(y_1) n=length(x); else disp('y和y的导数的维数不相等！'); return; end end else disp('x和y的维数不相等！); return; end for 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(i)-y_1(i)+y(i); if(i=n) if(nargin=4)