用matlab进行数据插值课件.ppt

关于用MATLAB进行数据插值第1页，此课件共107页哦 我们经常会遇到大量的数据需要处理，我们经常会遇到大量的数据需要处理，而处理数据的关键就在于这些算法，例如数而处理数据的关键就在于这些算法，例如数据拟合、参数估计、插值等数据处理算法。据拟合、参数估计、插值等数据处理算法。此类问题在此类问题在MATLAB中有很多现成的函数可以中有很多现成的函数可以调用，熟悉调用，熟悉MATLAB，这些方法都能游刃有余，这些方法都能游刃有余的用好。的用好。一、一、概述概述第2页，此课件共107页哦 数据拟合在很多赛题中有应用，与图形处数据拟合在很多赛题中有应用，与图形处理有关的问题很多与插值和拟合有关系，例如理有关的问题很多与插值和拟合有关系，例如9898年美国赛年美国赛A A题，生物组织切片的三维插值处理，题，生物组织切片的三维插值处理，9494年年A A题逢山开路，山体海拔高度的插值计算，题逢山开路，山体海拔高度的插值计算，20032003年吵的沸沸扬扬的年吵的沸沸扬扬的“非典非典”问题也要用到问题也要用到数据拟合算法，观察数据的走向进行处理，数据拟合算法，观察数据的走向进行处理，20052005年的雨量预报的评价的插值计算。年的雨量预报的评价的插值计算。20012001年的年的公交车调度拟合问题，公交车调度拟合问题，20032003年的饮酒驾车拟合问年的饮酒驾车拟合问题。题。第3页，此课件共107页哦预测点和实测点的图形预测点和实测点的图形插值后的图形插值后的图形第4页，此课件共107页哦喝两瓶酒的拟合曲线喝两瓶酒的拟合曲线喝喝1-5瓶酒的拟合曲线瓶酒的拟合曲线第5页，此课件共107页哦 在实际中，常常要处理由实验或测量所得在实际中，常常要处理由实验或测量所得到的一些离散数据。插值与拟合方法就是要通到的一些离散数据。插值与拟合方法就是要通过这些数据去确定某一类已知函数的参数或寻过这些数据去确定某一类已知函数的参数或寻求某个近似函数，使所得到的近似函数与已知求某个近似函数，使所得到的近似函数与已知数据有较高的拟合精度。数据有较高的拟合精度。如果要求这个近似函数（曲线或曲面）经过如果要求这个近似函数（曲线或曲面）经过所已知的所有数据点，则称此类问题为所已知的所有数据点，则称此类问题为插值问题插值问题。（不需要函数表达式）（不需要函数表达式）二、二、基本概念基本概念第6页，此课件共107页哦 如果不要求近似函数通过所有数据点，而是如果不要求近似函数通过所有数据点，而是要求它能较好地反映数据变化规律的近似函数的要求它能较好地反映数据变化规律的近似函数的方法称为方法称为数据拟合数据拟合。（必须有函数表达式）。（必须有函数表达式）近似函数不一定（曲线或曲面）通过所有近似函数不一定（曲线或曲面）通过所有的数据点。的数据点。第7页，此课件共107页哦1 1、联系、联系都是根据实际中一组已知数据来构造一个能够都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数的方法。反映数据变化规律的近似函数的方法。2 2、区别、区别插值问题插值问题不一定得到近似函数的表达形式，仅通不一定得到近似函数的表达形式，仅通过插值方法找到未知点对应的值。过插值方法找到未知点对应的值。数据拟合数据拟合要求要求得到一个具体的近似函数的表达式。得到一个具体的近似函数的表达式。三、插值与三、插值与拟合的区别和联系拟合的区别和联系第8页，此课件共107页哦 插插 值值 当数据量不够，需要补充，且认定已有数据可信时当数据量不够，需要补充，且认定已有数据可信时,通通常利用函数插值方法。常利用函数插值方法。实际问题当中碰到的函数实际问题当中碰到的函数 f(x)是各种各样的，有的表达是各种各样的，有的表达式很复杂，有的甚至给不出数学的式子，只提供了一些离散数式很复杂，有的甚至给不出数学的式子，只提供了一些离散数据，警如，某些点上的函数值和导数值。据，警如，某些点上的函数值和导数值。第9页，此课件共107页哦拉格朗日插值拉格朗日插值分段线性插值分段线性插值三次样条插值三次样条插值一一 维维 插插 值值一、插值的定义一、插值的定义二、插值的方法二、插值的方法三、用三、用Matlab解插值问题解插值问题返回返回第10页，此课件共107页哦返回返回二维插值二维插值一、二维插值定义一、二维插值定义二、网格节点插值法二、网格节点插值法三、用三、用MatlabMatlab解插值问题解插值问题最邻近插值最邻近插值分片线性插值分片线性插值双线性插值双线性插值网格节点数据的插值网格节点数据的插值散点数据的插值散点数据的插值第11页，此课件共107页哦一维插值的定义一维插值的定义已知已知 n+1个节点个节点其中其中互不相同，不妨设互不相同，不妨设求任一插值点求任一插值点处的插值处的插值节点可视为由节点可视为由产生产生,，表达式复杂表达式复杂,，或无封闭形式或无封闭形式,，或未知或未知.。第12页，此课件共107页哦 构造一个构造一个(相对简单的相对简单的)函数函数通过全部节点通过全部节点,即即再用再用计算插值，即计算插值，即返回返回第13页，此课件共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)插值插值第14页，此课件共107页哦拉格朗日拉格朗日(Lagrange)插值插值特别地特别地:两点一次两点一次(线性线性)插值多项式插值多项式:三点二次三点二次(抛物抛物)插值多项式插值多项式:第15页，此课件共107页哦 拉格朗日多项式插值的这种振荡现象叫 Runge现象现象 采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.例例返回返回To MatlabTo Matlablch(larg1)lch(larg1)第16页，此课件共107页哦第17页，此课件共107页哦分段线性插值分段线性插值n越大，误差越小.xjxj-1xj+1x0 xnxoy第18页，此课件共107页哦To MATLABxch11，xch12，xch13，xch14返回返回例例用分段线性插值法求插值用分段线性插值法求插值,并观察插值误差并观察插值误差.1.在在-6,6中平均选取中平均选取5个点作插值个点作插值(xch11)4.在在-6,6中平均选取中平均选取41个点作插值个点作插值(xch14)2.在在-6,6中平均选取中平均选取11个点作插值个点作插值(xch12)3.在在-6,6中平均选取中平均选取21个点作插值个点作插值(xch13)第19页，此课件共107页哦第20页，此课件共107页哦比分段线性插值更光滑。比分段线性插值更光滑。xyxi-1 xiab 在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值三次样条插值第21页，此课件共107页哦 三次样条插值g g(x x)为被插值函数为被插值函数。第22页，此课件共107页哦例例用三次样条插值选取用三次样条插值选取11个基点计算插值个基点计算插值(ych)返回返回To MATLABych第23页，此课件共107页哦用用MATLABMATLAB作插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linear ：线性插值；线性插值；spline ：三次样条插值；三次样条插值；cubic ：立方插值。立方插值。缺省时：缺省时：分段线性插值。分段线性插值。注意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调的，并且是单调的，并且xi不能够超过不能够超过x的范围。的范围。第24页，此课件共107页哦（1 1）nearestnearest方法速度最快，占用内存最小，但一般来说误方法速度最快，占用内存最小，但一般来说误差最大，插值结果最不光滑；差最大，插值结果最不光滑；(2)liner(2)liner分段线性插值：插值点处函数值由连接其最邻近的两侧点的线性函数预测，MATLAB中中interp1interp1的默认方的默认方法法（3 3）splinespline三次样条插值是所有插值方法中运行耗时最长的，其插值函数以及插值函数的一阶、二阶导函数都连续，因此是最光滑的插值方法，占用内存上比cubic方方法小，但当已知数据点不均匀分布时可能出现异常结果。法小，但当已知数据点不均匀分布时可能出现异常结果。（4 4）cubiccubic三次多项式插值法中插值函数及其一阶导数都三次多项式插值法中插值函数及其一阶导数都是连续的，因此其插值结果也比较光滑，运算速度比是连续的，因此其插值结果也比较光滑，运算速度比splinespline方法略快，但占用内存最多。在实际的使用中，应根方法略快，但占用内存最多。在实际的使用中，应根据实际需求和运算条件选择合适的算法。据实际需求和运算条件选择合适的算法。第25页，此课件共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)、(7,5.7)(7,5.7)进行一维插值方法。进行一维插值方法。进行一维插值方法。进行一维插值方法。解：解：在MATLABMATLAB命令窗口中输入以下命命令窗口中输入以下命令：令：x=1234567;x=1234567;y=3.52.11.30.82.94.25.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=interp1(x,y,xx,cubic);y3=interp1(x,y,xx,cubic);plot(x,y,o,xx,y1,-,xx,y2,-.,xx,y3,:)plot(x,y,o,xx,y1,-,xx,y2,-.,xx,y3,:)第26页，此课件共107页哦第27页，此课件共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 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)第28页，此课件共107页哦第29页，此课件共107页哦xy机翼下轮廓线例例 已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值。值。To MATLAB(plane)返回返回第30页，此课件共107页哦第31页，此课件共107页哦二维插值的定义二维插值的定义 xyO O第一种（网格节点）：第一种（网格节点）：第32页，此课件共107页哦 已知已知 m n个节点个节点 其中其中互不相同，不妨设互不相同，不妨设 构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，即计算插值，即第33页，此课件共107页哦第二种（散乱节点）：第二种（散乱节点）：yx0 0第34页，此课件共107页哦已知已知n个节点个节点其中其中互不相同，互不相同，构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，即计算插值，即返回返回第35页，此课件共107页哦 注意：注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值最邻近插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O 二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。返回返回第36页，此课件共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第37页，此课件共107页哦插值函数为：第二片(上三角形区域)：(x,y)满足插值函数为：注意注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足返回返回第38页，此课件共107页哦 双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值双线性插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O返回返回第39页，此课件共107页哦 要求要求x0,y0 x0,y0单调；单调；x x，y y可取可取为矩阵，或为矩阵，或x x取行向取行向量，量，y y取为列向量，取为列向量，x,yx,y的值分别不能超出的值分别不能超出x0,y0 x0,y0的范的范围。围。z=interp2(x0,y0,z0,x,y,method)被插值点插值方法用用MATLAB作网格节点数据的插值作网格节点数据的插值插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值cubiccubic 双三次插值双三次插值缺省时缺省时,双线性插值双线性插值第40页，此课件共107页哦例：测得平板表面例：测得平板表面3*53*5网格点处的温度分别为：网格点处的温度分别为：82 81 80 82 84 82 81 80 82 84 79 63 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个单位的地方进行插值.第41页，此课件共107页哦第42页，此课件共107页哦再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图.To MATLAB(wendu)第43页，此课件共107页哦第44页，此课件共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),第45页，此课件共107页哦选较密的插值点，用默认的线性插值算法进行插值x1,y1=meshgrid(-3:.2:3,-2:.2:2);z0=interp2(x,y,z,x1,y1);surf(x1,y1,z0)第46页，此课件共107页哦立方和样条插值：z1=interp2(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)第47页，此课件共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)第48页，此课件共107页哦山区地形地貌图已知某处山区地形选点测量坐标数据为：x=00.511.522.533.544.55y=00.511.522.533.544.555.56海拔高度数据为：z=8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287第49页，此课件共107页哦山区地形地貌图程序原始地貌图程序：x=0:.5:5;y=0:.5:6;xx,yy=meshgrid(x,y);z=8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287;mesh(xx,yy,z)加密后的地貌图x=0:.5:5;y=0:.5:6;z=8990878592919693908782929698999591898684828496989592908885848381858081828995969392898686828587989996978885828382858994959392918684888892939495898786838192929697989693958482818485858182808081859093958486819899989796958487808185828384879095868880828184858683828180828788899899979698949287;xi=linspace(0,5,50);%加密横坐标数据到50个yi=linspace(0,6,80);%加密纵坐标数据到60个xii,yii=meshgrid(xi,yi);%生成网格数据zii=interp2(x,y,z,xii,yii,cubic);%插值mesh(xii,yii,zii)%加密后的地貌图第50页，此课件共107页哦山区地形地貌图结果第51页，此课件共107页哦 通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。To MATLAB (moutain)返回返回第52页，此课件共107页哦第53页，此课件共107页哦第54页，此课件共107页哦第55页，此课件共107页哦第56页，此课件共107页哦第57页，此课件共107页哦 插值函数插值函数griddata格式为格式为:cz=griddata（x，y，z，cx，cy，method）用用MATLABMATLAB作散点数据的插值计算作散点数据的插值计算 要求要求cxcx取行向量，取行向量，cycy取为列向量取为列向量。被插值点插值方法插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值cubiccubic 双三次插值双三次插值v4-Matlab提供的插值方法提供的插值方法缺省时缺省时,双线性插值双线性插值第58页，此课件共107页哦x=-3+6*rand(200,1);y=-2+4*rand(200,1);z=(x.2-2*x).*exp(-x.2-y.2-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插值法进行处理，并对误差进行比较。第59页，此课件共107页哦第60页，此课件共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.15)figure;surf(x1,y1,abs(z0-z2),axis(-3,3,-2,2,0,0.15)第61页，此课件共107页哦第62页，此课件共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.5)，gridon第63页，此课件共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)第64页，此课件共107页哦新样本点拟合曲面x1,y1=meshgrid(-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)第65页，此课件共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)误差的二维等高线图第66页，此课件共107页哦命令3interp3三维网格生成用meshgrid()函数,调用格式:x,y,z=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()函函数一致。数一致。第67页，此课件共107页哦例：通过函数生成一些网格型样本点，试根据样本点进行拟合，并给出拟合误差。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第68页，此课件共107页哦slice(x0,y0,z0,V1,-0.5,0.3,0.9,0.6,-0.1,-1,-0.5,0.5,1)title(Slives for Four Dim Figures)第69页，此课件共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);第70页，此课件共107页哦 例例 在某海域测得一些点在某海域测得一些点(x,y)(x,y)处的水深处的水深z z由下表给出，由下表给出，船的吃水深度为船的吃水深度为5 5英尺，在矩形区域（英尺，在矩形区域（7575，200200）*（-50-50，150150）里的哪些地方船要避免进入。）里的哪些地方船要避免进入。第71页，此课件共107页哦To MATLAB hd1返回返回4.作出水深小于5的海域范围,即z=5的等高线.第72页，此课件共107页哦第73页，此课件共107页哦第74页，此课件共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:200;y1=-50:1:150;x1,y1=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,v4);meshc(x1,y1,z1)第75页，此课件共107页哦海底曲面图海底曲面图第76页，此课件共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)第77页，此课件共107页哦水深小于水深小于5 5的海域范围的海域范围第78页，此课件共107页哦拉格朗日插值法 拉格朗日插值法拉格朗日插值法是基于基函数的插值方法，插值多项式可表是基于基函数的插值方法，插值多项式可表示为示为其中其中称为称为i i次基函数：次基函数：第79页，此课件共107页哦在MATLAB中编程实现拉格朗日插值法拉格朗日插值法函数为：LanguageLanguage。功能：求已知数据点的拉格朗日多项式；功能：求已知数据点的拉格朗日多项式；调用格式：调用格式：f=f=Language(x,y)Language(x,y)或或f=f=Language(x,y,x0)Language(x,y,x0)。其中，其中，x x为已知数据点的为已知数据点的xx坐标向量；坐标向量；y y为已知数据点的为已知数据点的yy坐标向量；坐标向量；x0 x0为插值点的为插值点的x x坐标；f f为求得的拉格朗日多项式或为求得的拉格朗日多项式或x0处的插值。处的插值。第80页，此课件共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的维数不相等！);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第81页，此课件共107页哦例4-3根据下表的数据点求出其拉格朗日插值多项式，并计算当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-.204917*t3+.328112e-1*t4-1.05427*t-.145485e-1*t2-.204917*t3+.328112e-1*t4-.261189e-1.261189e-1f=language(x,y,1.6)f=language(x,y,1.6)f=f=0.99920.9992x x1 11.21.21.81.82.52.54 4yy0.84150.84150.93200.93200.9730.9738 80.59850.5985-0.7568-0.7568第82页，此课件共107页哦利用均差的牛顿插值法 函数函数 f f的零阶均差定义为的零阶均差定义为，一阶定义均差为，一阶定义均差为一般地，函数一般地，函数f f的的k k阶均差定义为：阶均差定义为：利用均差的牛顿插值法利用均差的牛顿插值法多项式为多项式为第83页，此课件共107页哦系数的计算过程如表所示表表 均差计算表格均差计算表格一阶均差一阶均差二阶均差二阶均差三阶均差三阶均差n n阶均差阶均差第84页，此课件共107页哦在MATLAB中编程实现利用均差牛顿插值法利用均差牛顿插值法函数为：NewtonNewton。功能：求已知数据点的均差形式的牛顿插值多项式；调用格式：调用格式：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为求得的牛顿插值法多项式或为求得的牛顿插值法多项式或x0 x0处的插值。处的插值。第85页，此课件共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)=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第86页，此课件共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*(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第87页，此课件共107页哦例根据下表的数据点求出其均差形式牛顿插值多项式，并计算当x=2.0时y的值。解：解：x=11.21.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-14*t3+.211471e-15*t4.182711e-14-.482154e-14*t+1.00000*t2-.169177e-14*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 4yy1 11.441.443.243.246.256.251616第88页，此课件共107页哦利用差分的牛顿插值法 差分分为向前差分、后向差分和中心差分三种，它们的记差分分为向前差分、后向差分和中心差分三种，它们的记法及定义如下为：法及定义如下为：其中：其中：代表向前差分；代表向前差分；代表向后差分；代表向后差分。代表向后差分。第89页，此课件共107页哦假设假设。为了方便，可构造如表所示。为了方便，可构造如表所示的差分表（的差分表（）。）。表表 差分计算表格差分计算表格第90页，此课件共107页哦向前牛顿插值 向前牛顿插值向前牛顿插值多项式可表示如下：多项式可表示如下：其中其中叫做步长，叫做步长，且且的取值范围为的取值范围为。第91