用MATLAB进行数据插值资料.pptx

《用MATLAB进行数据插值资料.pptx》由会员分享，可在线阅读，更多相关《用MATLAB进行数据插值资料.pptx（107页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 我们经常会遇到大量的数据需要处理，我们经常会遇到大量的数据需要处理，而处理数据的关键就在于这些算法，例如而处理数据的关键就在于这些算法，例如数据拟合、参数估计、插值等数据处理算数据拟合、参数估计、插值等数据处理算法。此类问题在法。此类问题在MATLAB中有很多现成的中有很多现成的函数可以调用，熟悉函数可以调用，熟悉MATLAB，这些方法，这些方法都能游刃有余的用好。都能游刃有余的用好。一、一、概述概述第1页/共107页 数据拟合在很多赛题中有应用，与图形数据拟合在很多赛题中有应用，与图形处理有关的问题很多与插值和拟合有关系，处理有关的问题很多与插值和拟合有关系，例如例如9898年美国赛年美国

2、赛A A题，生物组织切片的三维插题，生物组织切片的三维插值处理，值处理，9494年年A A题逢山开路，山体海拔高度的题逢山开路，山体海拔高度的插值计算，插值计算，20032003年吵的沸沸扬扬的年吵的沸沸扬扬的“非典非典”问题也要用到数据拟合算法，观察数据的走问题也要用到数据拟合算法，观察数据的走向进行处理，向进行处理，20052005年的雨量预报的评价的插年的雨量预报的评价的插值计算。值计算。20012001年的公交车调度拟合问题，年的公交车调度拟合问题，20032003年的饮酒驾车拟合问题。年的饮酒驾车拟合问题。第2页/共107页预测点和实测点的图形预测点和实测点的图形插值后的图形插值后的

3、图形第3页/共107页喝两瓶酒的拟合曲线喝1-5瓶酒的拟合曲线第4页/共107页 在实际中，常常要处理由实验或测量所在实际中，常常要处理由实验或测量所得到的一些离散数据。插值与拟合方法就是得到的一些离散数据。插值与拟合方法就是要通过这些数据去确定某一类已知函数的参要通过这些数据去确定某一类已知函数的参数或寻求某个近似函数，使所得到的近似函数或寻求某个近似函数，使所得到的近似函数与已知数据有较高的拟合精度。数与已知数据有较高的拟合精度。如果要求这个近似函数（曲线或曲面）如果要求这个近似函数（曲线或曲面）经过所已知的所有数据点，则称此类问题为经过所已知的所有数据点，则称此类问题为插值问题插值问题。

4、（不需要函数表达式）（不需要函数表达式）二、二、基本概念基本概念第5页/共107页 如果不要求近似函数通过所有数据点，如果不要求近似函数通过所有数据点，而是要求它能较好地反映数据变化规律的近而是要求它能较好地反映数据变化规律的近似函数的方法称为似函数的方法称为数据拟合数据拟合。（必须有函数。（必须有函数表达式）表达式）近似函数不一定（曲线或曲面）通过所近似函数不一定（曲线或曲面）通过所有的数据点。有的数据点。第6页/共107页1 1、联系、联系都是根据实际中一组已知数据来构造一个能够都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数的方法。反映数据变化规律的近似函数的方法。2

5、2、区别、区别插值问题插值问题不一定得到近似函数的表达形式，仅不一定得到近似函数的表达形式，仅通过插值方法找到未知点对应的值。通过插值方法找到未知点对应的值。数据拟合数据拟合要求得到一个具体的近似函数的表达式。要求得到一个具体的近似函数的表达式。三、插值与三、插值与拟合的区别和联系拟合的区别和联系第7页/共107页 插插 值值 当数据量不够，需要补充，且认定已有数据可信时,通常利用函数插值方法。实际问题当中碰到的函数 f(x)是各种各样的，有的表达式很复杂，有的甚至给不出数学的式子，只提供了一些离散数据，警如，某些点上的函数值和导数值。第8页/共107页拉格朗日插值分段线性插值三次样条插值一一

6、 维维 插插 值值一、一、插值的定义插值的定义二、插值的方法二、插值的方法三、用三、用Matlab解插值问题解插值问题返回第9页/共107页返回二维插值二维插值一、一、二维插值定义二维插值定义二、网格节点插值法二、网格节点插值法三、用三、用MatlabMatlab解插值问题解插值问题最邻近插值最邻近插值分片线性插值分片线性插值双线性插值双线性插值网格节点数据的插值网格节点数据的插值散点数据的插值散点数据的插值第10页/共107页一维插值的定义一维插值的定义已知已知 n+1个节点个节点其中其中互不相同，不妨设互不相同，不妨设求任一插值点求任一插值点处的插值处的插值节点可视为由节点可视为由产生产生

7、,，表达式复杂表达式复杂,，或无封闭形式或无封闭形式,，或未知或未知.。第11页/共107页 构造一个构造一个(相对简单的相对简单的)函数函数通过全部节点通过全部节点,即即再用再用计算插值，即计算插值，即返回第12页/共107页 称为拉格朗日插值基函数拉格朗日插值基函数。已知函数f(x)在n+1个点x0,x1,xn处的函数值为 y0,y1,yn。求一n次多项式函数Pn(x)，使其满足：Pn(xi)=yi,i=0,1,n.解决此问题的拉格朗日插值多项式公式如下其中Li(x)为n次多项式：拉格朗日拉格朗日(Lagrange)插值插值第13页/共107页拉格朗日拉格朗日(Lagrange)插值插值特

8、别地特别地:两点一次两点一次(线性线性)插值多项式插值多项式:三点二次三点二次(抛物抛物)插值多项式插值多项式:第14页/共107页 拉格朗日多项式插值的这种振荡现象叫 Runge现象现象 采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.例例返回To MatlabTo Matlablch(larg1)lch(larg1)第15页/共107页第16页/共107页分段线性插值分段线性插值n越大，误差越小.xjxj-1xj+1x0 xnxoy第17页/共107页To MATLABxch11，xch12，xch13，xc

9、h14返回例例用分段线性插值法求插值用分段线性插值法求插值,并观察插值误差并观察插值误差.1.在在-6,6中平均选取中平均选取5个点作插值个点作插值(xch11)4.在在-6,6中平均选取中平均选取41个点作插值个点作插值(xch14)2.在在-6,6中平均选取中平均选取11个点作插值个点作插值(xch12)3.在在-6,6中平均选取中平均选取21个点作插值个点作插值(xch13)第18页/共107页第19页/共107页比分段线性插值更光滑。比分段线性插值更光滑。xyxi-1 xiab 在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高

10、，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值三次样条插值第20页/共107页 三次样条插值g g(x x)为被插值函数为被插值函数。第21页/共107页例例用三次样条插值选取用三次样条插值选取11个基点计算插值个基点计算插值(ych)返回To MATLABych第22页/共107页用用MATLABMATLAB作插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linea

11、r ：线性插值；线性插值；spline ：三次样条三次样条插值；插值；cubic ：立方插值。立方插值。缺省时：缺省时：分段线性插值。分段线性插值。注意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调的，并且是单调的，并且xi不能够超过不能够超过x的范的范围。围。第23页/共107页 （1 1）nearestnearest方法速度最快，占用内存最小，但一般方法速度最快，占用内存最小，但一般来说误差最大，插值结果最不光滑；来说误差最大，插值结果最不光滑；(2)liner (2)liner分段线性插值：插值点处函数值由连接其分段线性插值：插值点处函数值由连接其最邻近的两侧点的线性函数预

12、测，最邻近的两侧点的线性函数预测，MATLABMATLAB中中interp1interp1的默认方法的默认方法 （3 3）splinespline三次样条插值是所有插值方法中运行耗三次样条插值是所有插值方法中运行耗时最长的，其插值函数以及插值函数的一阶、二阶导时最长的，其插值函数以及插值函数的一阶、二阶导函数都连续，因此是最光滑的插值方法，占用内存上函数都连续，因此是最光滑的插值方法，占用内存上比比cubiccubic方法小，但当已知数据点不均匀分布时可能出方法小，但当已知数据点不均匀分布时可能出现异常结果。现异常结果。（4 4）cubiccubic三次多项式插值法中插值函数及其一阶导三次多项

13、式插值法中插值函数及其一阶导数都是连续的，因此其插值结果也比较光滑，运算速数都是连续的，因此其插值结果也比较光滑，运算速度比度比splinespline方法略快，但占用内存最多。在实际的使用方法略快，但占用内存最多。在实际的使用中，应根据实际需求和运算条件选择合适的算法。中，应根据实际需求和运算条件选择合适的算法。第24页/共107页例例 用其他一维插值方法对以下用其他一维插值方法对以下7 7个离散数据点个离散数据点(1,3.5)(1,3.5)、(2,2.1)(2,2.1)、(3,1.3)(3,1.3)、(4.0.8)(4.0.8)、(5,2.9)(5,2.9)、(6,4.2)(6,4.2)、

14、(7,5.7)(7,5.7)进行一维插值方法。进行一维插值方法。解：解：在在MATLABMATLAB命令窗口中输入以下命令窗口中输入以下命令：命令：x=1 2 3 4 5 6 7;x=1 2 3 4 5 6 7;y=3.5 2.1 1.3 0.8 2.9 4.2 5.7;y=3.5 2.1 1.3 0.8 2.9 4.2 5.7;xx=1:0.5:7;xx=1:0.5:7;y1=interp1(x,y,xx,nearest);y1=interp1(x,y,xx,nearest);y2=interp1(x,y,xx,spline);y2=interp1(x,y,xx,spline);y3=int

15、erp1(x,y,xx,cubic);y3=interp1(x,y,xx,cubic);plot(x,y,o,xx,y1,-,xx,y2,-plot(x,y,o,xx,y1,-,xx,y2,-.,xx,y3,:).,xx,y3,:)第25页/共107页第26页/共107页 例：在例：在1-121-12的的1111小时内，每隔小时内，每隔1 1小时测量一次温小时测量一次温度，测得的温度依次为：度，测得的温度依次为：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424。试估计每隔。试估计每隔1/101/10小时的温小时的温度值。度值。To

16、 MATLAB(temp)hours=1:12;temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12;t=interp1(hours,temps,h,spline);(直接输出数据将是很多的)plot(hours,temps,+)hold onplot(h,t,hours,temps,r:)%作图xlabel(Hour),ylabel(Degrees Celsius)第27页/共107页第28页/共107页xy机翼下轮廓线例例 已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值。值。To MATLAB(pl

17、ane)返回第29页/共107页第30页/共107页二维插值的定义二维插值的定义 xyO O第一种（网格节点）：第一种（网格节点）：第31页/共107页 已知已知 m n个节点个节点 其中其中互不相同，不妨设互不相同，不妨设 构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，即计算插值，即第32页/共107页第二种（散乱节点）：第二种（散乱节点）：yx0 0第33页/共107页已知已知n个节点个节点其中其中互不相同，互不相同，构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，即计算插值，即返回第34页/共107页 注意

18、：注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值最邻近插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O 二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。返回第35页/共107页 将四个插值点（矩形的四个顶点）处的函数值依次简记为：分片线性插值分片线性插值xy (xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)O Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4第36页/共107页插值函数为：第二片(上三角形区域)：(x,y)满足插值函

19、数为：注意注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足返回第37页/共107页 双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值双线性插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O返回第38页/共107页 要求要求x0,y0 x0,y0单调；单调；x x，y y可取可取为矩阵，或为矩阵，或x x取取行向量，行向量，y y取为列向量，取

20、为列向量，x,yx,y的值分别不能超出的值分别不能超出x0,y0 x0,y0的范围。的范围。z=interp2(x0,y0,z0,x,y,method)被插值点插值方法用用MATLAB作网格节点数据的插值作网格节点数据的插值插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值cubiccubic 双三次插值双三次插值缺省时缺省时,双线性插值双线性插值第39页/共107页例：测得平板表面例：测得平板表面3*53*5网格点处的温度分别为：网格点处的温度分别为：82 81 80 82 84 82 81 80 82 84 79 63

21、61 65 81 79 63 61 65 81 84 84 82 85 86 84 84 82 85 86 试作出平板表面的温度分布曲面试作出平板表面的温度分布曲面z=f(x,y)z=f(x,y)的图形。的图形。输入以下命令：x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,temps)1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.2以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.第40页/共107页第41页/共107页再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi

22、=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图.To MATLAB(wendu)第42页/共107页第43页/共107页例x,y=meshgrid(-3:.6:3,-2:.4:2);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);surf(x,y,z),第44页/共107页选较密的插值点，用默认的线性插值算法进行插值x1,y1=meshgrid(-3:.2:3,-2:.2:2);z0=interp2(x,y,z,x1,y1);surf(x1,y1,z0)第45页/共107页立方和样条插值：z1=interp

23、2(x,y,z,x1,y1,cubic);z2=interp2(x,y,z,x1,y1,spline);surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5)figure;surf(x1,y1,z2),axis(-3,3,-2,2,-0.7,1.5)第46页/共107页算法误差比较z=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1);surf(x1,y1,abs(z-z1)figure;surf(x1,y1,abs(z-z2)figure;surf(x1,y1,abs(z-z0)第47页/共107页山区地形地貌图已知某处山区地形选点测量坐标数据为：

24、x=00.511.522.533.544.55y=00.511.522.533.544.555.56海拔高度数据为：z=899087859291969390878292969899959189868482849698959290888584838185808182899596939289868682858798999697888582838285899495939291868488889293949589878683819292969798969395848281848585818280808185909395848681989998979695848780818582838487909586

25、8880828184858683828180828788899899979698949287第48页/共107页山区地形地貌图 程序原始地貌图程序：x=0:.5:5;y=0:.5:6;xx,yy=meshgrid(x,y);z=899087859291969390878292969899959189868482849698959290888584838185808182899596939289868682858798999697888582838285899495939291868488889293949589878683819292969798969395848281848585818280

26、8081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287;mesh(xx,yy,z)加密后的地貌图x=0:.5:5;y=0:.5:6;z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85

27、82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87;xi=linspace(0,5,50);%加密横坐标数据到5

28、0个yi=linspace(0,6,80);%加密纵坐标数据到60个xii,yii=meshgrid(xi,yi);%生成网格数据zii=interp2(x,y,z,xii,yii,cubic);%插值mesh(xii,yii,zii)%加密后的地貌图第49页/共107页山区地形地貌图 结果第50页/共107页 通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。To MATLAB (moutain)返回第51页/共107页第52页/共107页第53页/共107页第54页/共107页第55页/共107页第56页/共107页 插值函数插值函数griddata格式为格式为:c

29、z=griddata（x，y，z，cx，cy，method）用用MATLABMATLAB作散点数据的插值计算作散点数据的插值计算 要求要求cxcx取行向量，取行向量，cycy取为列向量取为列向量。被插值点插值方法插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值cubiccubic 双三次插值双三次插值v4-Matlab提供的插值方法提供的插值方法缺省时缺省时,双线性插值双线性插值第57页/共107页x=-3+6*rand(200,1);y=-2+4*rand(200,1);z=(x.2-2*x).*exp(-x.2-y.2

30、-x.*y);x1,y1=meshgrid(-3:.2:3,-2:.2:2);z1=griddata(x,y,z,x1,y1,cubic);surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5)z2=griddata(x,y,z,x1,y1,v4);figure;surf(x1,y1,z2),axis(-3,3,-2,2,-0.7,1.5)例：在x为3，3，y为2，2矩形区域随机选择一组坐标，用 v4 与cubic插值法进行处理，并对误差进行比较。第58页/共107页第59页/共107页误差分析z0=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1)

31、;surf(x1,y1,abs(z0-z1),axis(-3,3,-2,2,0,0.15)figure;surf(x1,y1,abs(z0-z2),axis(-3,3,-2,2,0,0.15)第60页/共107页第61页/共107页例：在x为3，3，y为2，2矩形区域随机选择一组坐标中，对分布不均匀数据，进行插值分析。x=-3+6*rand(200,1);y=-2+4*rand(200,1);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);%生成已知数据plot(x,y,x)%样本点的二维分布figure,plot3(x,y,z,x),axis(-3,3,-2,2,-0.7,1

32、.5)，gridon第62页/共107页去除在(-1,-1/2)点为圆心，以0.5为半径的圆内的点。x=-3+6*rand(200,1);y=-2+4*rand(200,1);%重新生成样本点z=(x.2-2*x).*exp(-x.2-y.2-x.*y);ii=find(x+1).2+(y+0.5).20.52);%找出不在圆内的点坐标x=x(ii);y=y(ii);z=z(ii);plot(x,y,x)t=0:.1:2*pi,2*pi;x0=-1+0.5*cos(t);y0=-0.5+0.5*sin(t);line(x0,y0)第63页/共107页新样本点拟合曲面x1,y1=meshgrid

33、(-3:.2:3,-2:.2:2);z1=griddata(x,y,z,x1,y1,v4);surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5)第64页/共107页误差分析z0=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1);surf(x1,y1,abs(z0-z1),axis(-3,3,-2,2,0,0.1)contour(x1,y1,abs(z0-z1),30);holdon,plot(x,y,x);line(x0,y0)误差的二维等高线图第65页/共107页命令3interp3三维网格生成用meshgrid()函数,调用格式:x,y,z=

34、meshgrid(x1,y1,z1)其中x1,y1,z1为这三维所需要的分割形式，应以向量形式给出，返回x,y,z为网格的数据生成，均为三维数组。griddata3()三维非网格形式的插值拟合命令4interpnn维网格生成用ndgrid()函数,调用格式:x1,x2,xn=ndgridv1,v2,vngriddatan()n维非网格形式的插值拟合interp3()、interpn()调用格式同interp2()函数一致；griddata3()、griddatan()调用格式同griddata()函数一致。第66页/共107页例：通过函数生成一些网格型样本点，试根据样本点进行拟合，并给出拟合误

35、差。x,y,z=meshgrid(-1:0.2:1);x0,y0,z0=meshgrid(-1:0.05:1);V=exp(x.2.*z+y.2.*x+z.2.*y).*cos(x.2.*y.*z+z.2.*y.*x);V0=exp(x0.2.*z0+y0.2.*x0+z0.2.*y0).*cos(x0.2.*y0.*z0+z0.2.*y0.*x0);V1=interp3(x,y,z,V,x0,y0,z0,spline);err=V1-V0;max(err(:)ans=0.0419第67页/共107页slice(x0,y0,z0,V1,-0.5,0.3,0.9,0.6,-0.1,-1,-0.5

36、,0.5,1)title(Slives for Four Dim Figures)第68页/共107页x,y,z=meshgrid(-3:.1:3,-4:.2:4,-4:.4:4);v=x.*exp(-x.2-y.2-z.2);xslice=-2.2,.5,3;yslice=3,0.9;zslice=-4,1;slice(x,y,z,v,xslice,yslice,zslice);第69页/共107页 例 在某海域测得一些点(x,y)(x,y)处的水深z z由下表给出，船的吃水深度为5 5英尺，在矩形区域（7575，200200）*（-50-50，150150）里的哪些地方船要避免进入。第70

37、页/共107页To MATLAB hd1返回4.作出水深小于5的海域范围,即z=5的等高线.第71页/共107页第72页/共107页第73页/共107页%程序一：插值并作海底曲面图程序一：插值并作海底曲面图 x =129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5;y=7.5 141.5 23.0 147.0 22.5 137.5 85.5 -6.5 -81 3.0 56.5 -66.5 84.0 -33.5;z=4 8 6 8 6 8 8 9 9 8 8 9 4 9;x1=75:1:

38、200;y1=-50:1:150;x1,y1=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,v4);meshc(x1,y1,z1)第74页/共107页海底曲面图海底曲面图第75页/共107页%程序二：插值并作出水深小于程序二：插值并作出水深小于5 5的海域范围。的海域范围。x1=75:1:200;y1=-50:1:150;x1,y1=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,v4);%插值z1(z1=5)=nan;%将水深大于5的置为nan，这样绘图就不会显示出来meshc(x1,y1,z1)第76页/共107页水深小于水

39、深小于5 5的海域范围的海域范围第77页/共107页拉格朗日插值法 拉格朗日插值法拉格朗日插值法是基于基函数的插值方法，插值多是基于基函数的插值方法，插值多项式可表示为项式可表示为其中其中称为称为i i次基函数：次基函数：第78页/共107页 在在MATLABMATLAB中编程实现中编程实现拉格朗日插值法拉格朗日插值法函数为：函数为：LanguageLanguage。功能：求已知数据点的拉格朗日多项式；功能：求已知数据点的拉格朗日多项式；调用格式：调用格式：f=f=Language(x,y)Language(x,y)或或f=f=Language(x,y,x0)Language(x,y,x0)。

40、其中，其中，x x为已知数据点的为已知数据点的xx坐标向量；坐标向量；y y为已知数据点的为已知数据点的yy坐标向量；坐标向量；x0 x0为插值点的为插值点的x x坐标；坐标；f f为求得的拉格朗日多项式或为求得的拉格朗日多项式或x0 x0处的处的插值。插值。第79页/共107页functionf=Language(x,y,x0)%求已知数据点的拉格朗日多项式%已知数据点的x坐标向量:x%已知数据点的y坐标向量:y%为插值点的x坐标:x0%求得的拉格朗日多项式或x0处的插值:fsymst;if(length(x)=length(y)n=length(x);elsedisp(x和y的维数不相等！

41、);return;end%检错f=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);%计算插值点的函数值elsef=collect(f);%将插值多项式展开f=vpa(f,6);%将插值多项式的系数化成6位精度的小数endendend第80页/共107页例4-3 根据下表的数据点求出其拉格朗日插值多

42、项式，并计算当x=1.6时y的值。解：解：x=11.21.82.54;x=11.21.82.54;y=0.84150.93200.97380.5985-0.7568;y=0.84150.93200.97380.5985-0.7568;f=language(x,y)f=language(x,y)f=f=1.05427*t-.145485e-1*t2-1.05427*t-.145485e-1*t2-.204917*t3+.328112e-1*t4-.261189e-1.204917*t3+.328112e-1*t4-.261189e-1f=language(x,y,1.6)f=language(x

43、,y,1.6)f=f=0.99920.9992x x1 11.21.21.81.82.52.54 4 y y0.84150.84150.93200.93200.9730.9738 80.59850.5985-0.7568-0.7568第81页/共107页利用均差的牛顿插值法 函数函数 f f的零阶均差定义为的零阶均差定义为，一阶定义均，一阶定义均差为差为一般地，函数一般地，函数f f的的k k阶均差定义为：阶均差定义为：利用均差的牛顿插值法利用均差的牛顿插值法多项式为多项式为第82页/共107页系数的计算过程如表所示表表 均差计算表格均差计算表格一阶均差一阶均差二阶均差二阶均差三阶均差三阶均差

44、n n阶均差阶均差第83页/共107页 在在MATLABMATLAB中编程实现中编程实现利用均差牛顿插值法利用均差牛顿插值法函数为：函数为：NewtonNewton。功能：求已知数据点的均差形式的牛顿插值多项式；功能：求已知数据点的均差形式的牛顿插值多项式；调用格式：调用格式：f=f=Newton(x,y)Newton(x,y)或或f=f=Newton(x,y,x0)Newton(x,y,x0)。其中，其中，x x为已知数据点的为已知数据点的xx坐标向量；坐标向量；y y为已知数据点的为已知数据点的yy坐标向量；坐标向量；x0 x0为插值点的为插值点的x x坐标；坐标；f f为求得的牛顿插值法

45、多项式或为求得的牛顿插值法多项式或x0 x0处的插处的插值。值。第84页/共107页functionf=Newton(x,y,x0)functionf=Newton(x,y,x0)%求已知数据点的均差形式牛顿插值多项式求已知数据点的均差形式牛顿插值多项式%已知数据点的已知数据点的xx坐标向量坐标向量:x:x%已知数据点的已知数据点的yy坐标向量坐标向量:y:y%为插值点的为插值点的x x坐标坐标:x0:x0%求得的均差形式牛顿插值多项式或求得的均差形式牛顿插值多项式或x0 x0处的插值处的插值:f:fsymst;symst;if(length(x)=length(y)if(length(x)=

46、length(y)n=length(x);n=length(x);c(1:n)=0.0;c(1:n)=0.0;elseelsedisp(xdisp(x和和y y的维数不相等！的维数不相等！););return;return;endend第85页/共107页f=y(1);f=y(1);y1=0;y1=0;l=1;l=1;for(i=1:n-1)for(i=1:n-1)for(j=i+1:n)for(j=i+1:n)y1(j)=(y(j)-y(i)/(x(j)-x(i);y1(j)=(y(j)-y(i)/(x(j)-x(i);endendc(i)=y1(i+1);c(i)=y1(i+1);l=l*

47、(t-x(i);l=l*(t-x(i);f=f+c(i)*l;f=f+c(i)*l;simplify(f);simplify(f);y=y1;y=y1;if(i=n-1)if(i=n-1)if(nargin=3)if(nargin=3)f=subs(f,t,x0);f=subs(f,t,x0);elseelsef=collect(f);%f=collect(f);%将插值多项式展开将插值多项式展开f=vpa(f,6);f=vpa(f,6);endendendendendend第86页/共107页例 根据下表的数据点求出其均差形式牛顿插值多项式，并计算当x=2.0时y的值。解：解：x=11.21

48、.82.54;x=11.21.82.54;y=11.443.246.2516;y=11.443.246.2516;f=Newton(x,y)f=Newton(x,y)f=f=.182711e-14-.482154e-14*t+1.00000*t2-.169177e-.182711e-14-.482154e-14*t+1.00000*t2-.169177e-14*t3+.211471e-15*t414*t3+.211471e-15*t4f=Newton(x,y,2.0)f=Newton(x,y,2.0)f=f=44x x1 11.21.21.81.82.52.54 4 y y1 11.441.4

49、43.243.246.256.251616第87页/共107页利用差分的牛顿插值法差分分为向前差分、后向差分和中心差分三种，差分分为向前差分、后向差分和中心差分三种，它们的记法及定义如下为：它们的记法及定义如下为：其中：其中：代表向前差分；代表向前差分；代表向后差分；代表向后差分；代表向后差分。代表向后差分。第88页/共107页假设假设。为了。为了方便，可构造如表所示的差分表（方便，可构造如表所示的差分表（）。）。表表 差分计算表格差分计算表格第89页/共107页向前牛顿插值 向前牛顿插值向前牛顿插值多项式可表示如下：多项式可表示如下：其中其中叫做步长，叫做步长，且且的取的取值范围为值范围为。

50、第90页/共107页 在在MATLABMATLAB中编程实现中编程实现向前牛顿插值法向前牛顿插值法函数为：函数为：NewtonforwardNewtonforward。功能：求已知数据点的向前牛顿插值法多项式；功能：求已知数据点的向前牛顿插值法多项式；调用格式：调用格式：f=f=Newtonforward(x,y)Newtonforward(x,y)或或f=f=NewtonforwardNewtonforward(x,y,x0)(x,y,x0)。其中，其中，x x为已知数据点的为已知数据点的xx坐标向量；坐标向量；y y为已知数据点的为已知数据点的yy坐标向量；坐标向量；x0 x0为插值点的为