用MATLAB进行数据插值资料

我们经常会遇到大量的数据需要处理，而处理数据的关键就在于这些算法，例如数据拟合、参数估计、插值等数据处理算法。此类问题在MATLAB中有很多现成的函数可以调用，熟悉MATLAB，这些方法都能游刃有余的用好。一、概述第1页/共107页第一页，共108页。数据拟合在很多赛题中有应用，与图形处理有关的问题很多与插值和拟合有关系，例如98年美国赛A题，生物组织切片的三维插值处理，94年A题逢山开路，山体海拔高度的插值计算，2003年吵的沸沸扬扬的“非典”问题也要用到数据拟合算法，观察数据的走向进行处理，2005年的雨量预报的评价的插值计算。2001年的公交车调度拟合问题，2003年的饮酒驾车拟合问题。第2页/共107页第二页，共108页。插值问题——雨量预报的评价预测点和实测点的图形插值后的图形第3页/共107页第三页，共108页。拟合问题——饮酒驾车喝两瓶酒的拟合曲线喝1-5瓶酒的拟合曲线第4页/共107页第四页，共108页。在实际中，常常要处理由实验或测量所得到的一些离散数据。插值与拟合方法就是要通过这些数据去确定某一类已知函数的参数或寻求某个近似函数，使所得到的近似函数与已知数据有较高的拟合精度。如果要求这个近似函数（曲线或曲面）经过所已知的所有数据点，则称此类问题为插值问题。（不需要函数表达式）二、基本概念第5页/共107页第五页，共108页。如果不要求近似函数通过所有数据点，而是要求它能较好地反映数据变化规律的近似函数的方法称为数据拟合。（必须有函数表达式）近似函数不一定（曲线或曲面）通过所有的数据点。第6页/共107页第六页，共108页。1、联系都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数的方法。2、区别插值问题不一定得到近似函数的表达形式，仅通过插值方法找到未知点对应的值。数据拟合要求得到一个具体的近似函数的表达式。

三、插值与拟合的区别和联系第7页/共107页第七页，共108页。

插值当数据量不够，需要补充，且认定已有数据可信时,通常利用函数插值方法。实际问题当中碰到的函数f(x)是各种各样的，有的表达式很复杂，有的甚至给不出数学的式子，只提供了一些离散数据，警如，某些点上的函数值和导数值。第8页/共107页第八页，共108页。拉格朗日插值分段线性插值三次样条插值一维插值一、插值的定义二、插值的方法三、用Matlab解插值问题返回第9页/共107页第九页，共108页。返回二维插值一、二维插值定义二、网格节点插值法三、用Matlab解插值问题最邻近插值分片线性插值双线性插值网格节点数据的插值散点数据的插值第10页/共107页第十页，共108页。一维插值的定义已知n+1个节点其中互不相同，不妨设求任一插值点处的插值节点可视为由产生,，表达式复杂,，或无封闭形式,，或未知.。第11页/共107页第十一页，共108页。

构造一个(相对简单的)函数通过全部节点,即再用计算插值，即返回第12页/共107页第十二页，共108页。

称为拉格朗日插值基函数。

已知函数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页第十三页，共108页。拉格朗日(Lagrange)插值特别地:两点一次(线性)插值多项式:三点二次(抛物)插值多项式:第14页/共107页第十四页，共108页。

拉格朗日多项式插值的这种振荡现象叫Runge现象

采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.例返回ToMatlablch(larg1)第15页/共107页第十五页，共108页。第16页/共107页第十六页，共108页。分段线性插值n越大，误差越小.xjxj-1xj+1x0xnxoy第17页/共107页第十七页，共108页。ToMATLABxch11，xch12，xch13，xch14返回例用分段线性插值法求插值,并观察插值误差.1.在[-6,6]中平均选取5个点作插值(xch11)4.在[-6,6]中平均选取41个点作插值(xch14)2.在[-6,6]中平均选取11个点作插值(xch12)3.在[-6,6]中平均选取21个点作插值(xch13)第18页/共107页第十八页，共108页。第19页/共107页第十九页，共108页。比分段线性插值更光滑。xyxi-1xiab

在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值第20页/共107页第二十页，共108页。

三次样条插值g(x)为被插值函数。第21页/共107页第二十一页，共108页。例用三次样条插值选取11个基点计算插值(ych)返回ToMATLABych第22页/共107页第二十二页，共108页。用MATLAB作插值计算一维插值函数：yi=interp1(x，y，xi，'method')插值方法被插值点插值节点xi处的插值结果‘nearest’

：最邻近插值‘linear’

：线性插值；‘spline’

：三次样条插值；‘cubic’

：立方插值。缺省时：分段线性插值。注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。第23页/共107页第二十三页，共108页。（1）nearest方法速度最快，占用内存最小，但一般来说误差最大，插值结果最不光滑；(2)‘liner’分段线性插值：插值点处函数值由连接其最邻近的两侧点的线性函数预测，MATLAB中interp1的默认方法（3）spline三次样条插值是所有插值方法中运行耗时最长的，其插值函数以及插值函数的一阶、二阶导函数都连续，因此是最光滑的插值方法，占用内存上比cubic方法小，但当已知数据点不均匀分布时可能出现异常结果。（4）cubic三次多项式插值法中插值函数及其一阶导数都是连续的，因此其插值结果也比较光滑，运算速度比spline方法略快，但占用内存最多。在实际的使用中，应根据实际需求和运算条件选择合适的算法。第24页/共107页第二十四页，共108页。例用其他一维插值方法对以下7个离散数据点

(1,3.5)、(2,2.1)、(3,1.3)、(4.0.8)、(5,2.9)、(6,4.2)、(7,5.7)进行一维插值方法。解：在MATLAB命令窗口中输入以下命令：

>>x=[1234567];>>y=[3.52.11.30.82.94.25.7];>>xx=1:0.5:7;>>y1=interp1(x,y,xx,'nearest');>>y2=interp1(x,y,xx,'spline');>>y3=interp1(x,y,xx,'cubic');>>plot(x,y,'o',xx,y1,'-',xx,y2,'-.',xx,y3,':')第25页/共107页第二十五页，共108页。第26页/共107页第二十六页，共108页。

例：在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。ToMATLAB(temp)hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接输出数据将是很多的)plot(hours,temps,'+‘)holdonplot(h,t,hours,temps,'r:')%作图xlabel('Hour'),ylabel('DegreesCelsius’)第27页/共107页第二十七页，共108页。第28页/共107页第二十八页，共108页。xy机翼下轮廓线例已知飞机下轮廓线上数据如下，求x每改变0.1时的y值。ToMATLAB(plane)返回第29页/共107页第二十九页，共108页。第30页/共107页第三十页，共108页。二维插值的定义xyO第一种（网格节点）：第31页/共107页第三十一页，共108页。

已知mn个节点其中互不相同，不妨设

构造一个二元函数通过全部已知节点,即再用计算插值，即第32页/共107页第三十二页，共108页。第二种（散乱节点）：yx0第33页/共107页第三十三页，共108页。已知n个节点其中互不相同，

构造一个二元函数通过全部已知节点,即再用计算插值，即返回第34页/共107页第三十四页，共108页。

注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值xy(x1,y1)(x1,y2)(x2,y1)(x2,y2)O

二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。返回第35页/共107页第三十五页，共108页。

将四个插值点（矩形的四个顶点）处的函数值依次简记为：

分片线性插值xy(xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4第36页/共107页第三十六页，共108页。插值函数为：第二片(上三角形区域)：(x,y)满足插值函数为：注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足返回第37页/共107页第三十七页，共108页。

双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值xy(x1,y1)(x1,y2)(x2,y1)(x2,y2)O返回第38页/共107页第三十八页，共108页。

要求x0,y0单调；x，y可取为矩阵，或x取行向量，y取为列向量，x,y的值分别不能超出x0,y0的范围。z=interp2(x0,y0,z0,x,y,’method’)被插值点插值方法用MATLAB作网格节点数据的插值插值节点被插值点的函数值‘nearest’

最邻近插值‘linear’

双线性插值‘cubic’

双三次插值缺省时,双线性插值第39页/共107页第三十九页，共108页。例：测得平板表面3*5网格点处的温度分别为：828180828479636165818484828586试作出平板表面的温度分布曲面z=f(x,y)的图形。输入以下命令：x=1:5;y=1:3;temps=[8281808284;7963616581;8484828586];mesh(x,y,temps)1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.2．以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.第40页/共107页第四十页，共108页。第41页/共107页第四十一页，共108页。再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)画出插值后的温度分布曲面图.ToMATLAB(wendu)第42页/共107页第四十二页，共108页。第43页/共107页第四十三页，共108页。例>>[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页第四十四页，共108页。选较密的插值点，用默认的线性插值算法进行插值>>[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);>>z0=interp2(x,y,z,x1,y1);>>surf(x1,y1,z0)第45页/共107页第四十五页，共108页。立方和样条插值：>>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])第46页/共107页第四十六页，共108页。算法误差比较>>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页第四十七页，共108页。山区地形地貌图已知某处山区地形选点测量坐标数据为：x=0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5y=0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6海拔高度数据为：z=8990878592919693908782

9296989995918986848284

9698959290888584838185

8081828995969392898686

8285879899969788858283

8285899495939291868488

8892939495898786838192

9296979896939584828184

8585818280808185909395

8486819899989796958487

8081858283848790958688

8082818485868382818082

8788899899979698949287第48页/共107页第四十八页，共108页。山区地形地貌图程序原始地貌图程序：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)%加密后的地貌图第49页/共107页第四十九页，共108页。山区地形地貌图结果第50页/共107页第五十页，共108页。

通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。ToMATLAB(moutain)返回第51页/共107页第五十一页，共108页。第52页/共107页第五十二页，共108页。第53页/共107页第五十三页，共108页。第54页/共107页第五十四页，共108页。第55页/共107页第五十五页，共108页。第56页/共107页第五十六页，共108页。

插值函数griddata格式为:

cz

=griddata（x，y，z，cx，cy，‘method’）用MATLAB作散点数据的插值计算

要求cx取行向量，cy取为列向量。被插值点插值方法插值节点被插值点的函数值‘nearest’

最邻近插值‘linear’

双线性插值‘cubic’

双三次插值'v4'-Matlab提供的插值方法缺省时,双线性插值第57页/共107页第五十七页，共108页。>>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’插值法进行处理，并对误差进行比较。第58页/共107页第五十八页，共108页。第59页/共107页第五十九页，共108页。误差分析>>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])第60页/共107页第六十页，共108页。第61页/共107页第六十一页，共108页。例：

在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第62页/共107页第六十二页，共108页。去除在(-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).^2>0.5^2);%找出不在圆内的点坐标>>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页第六十三页，共108页。新样本点拟合曲面>>[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])第64页/共107页第六十四页，共108页。误差分析>>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页第六十五页，共108页。命令3

interp3

三维网格生成用meshgrid()函数,调用格式:[x,y,z]=meshgrid(x1,y1,z1)

其中x1,y1,z1为这三维所需要的分割形式，应以向量形式给出，返回x,y,z为网格的数据生成，均为三维数组。

griddata3()三维非网格形式的插值拟合命令4

interpnn维网格生成用ndgrid()函数,调用格式:[x1,x2,…,xn]=ndgrid[v1,v2,…,vn]griddatan()n维非网格形式的插值拟合interp3()、interpn()调用格式同interp2()函数一致；griddata3()、griddatan()调用格式同griddata()函数一致。第66页/共107页第六十六页，共108页。例：通过函数生成一些网格型样本点，试根据样本点进行拟合，并给出拟合误差。>>[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页第六十七页，共108页。slice(x0,y0,z0,V1,[-0.5,0.3,0.9],[0.6,-0.1],[-1,-0.5,0.5,1])title('SlivesforFourDimFigures')第68页/共107页第六十八页，共108页。[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页第六十九页，共108页。

例在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入。第70页/共107页第七十页，共108页。ToMATLABhd1返回4.作出水深小于5的海域范围,即z=5的等高线.第71页/共107页第七十一页，共108页。第72页/共107页第七十二页，共108页。第73页/共107页第七十三页，共108页。%程序一：插值并作海底曲面图

x=[129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5];y=[7.5141.523.0147.022.5137.585.5-6.5-813.056.5-66.584.0-33.5];z=[48686889988949];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)第74页/共107页第七十四页，共108页。海底曲面图第75页/共107页第七十五页，共108页。%程序二：插值并作出水深小于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页第七十六页，共108页。水深小于5的海域范围第77页/共107页第七十七页，共108页。拉格朗日插值法拉格朗日插值法是基于基函数的插值方法，插值多项式可表示为其中称为i次基函数：第78页/共107页第七十八页，共108页。在MATLAB中编程实现拉格朗日插值法函数为：Language。功能：求已知数据点的拉格朗日多项式；调用格式：f=

Language(x,y)或f=

Language(x,y,x0)。其中，x为已知数据点的x坐标向量；

y为已知数据点的y坐标向量；

x0为插值点的x坐标；

f为求得的拉格朗日多项式或x0处的插值。第79页/共107页第七十九页，共108页。

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第80页/共107页第八十页，共108页。例4-3根据下表的数据点求出其拉格朗日插值多项式，并计算当x=1.6时y的值。解：>>

x=[11.21.82.54];>>y=[0.84150.93200.97380.5985-0.7568];>>f=language(x,y)f=1.05427*t-.145485e-1*t^2-.204917*t^3+.328112e-1*t^4-.261189e-1>>f=language(x,y,1.6)f=0.9992x11.21.82.54y0.84150.93200.97380.5985-0.7568第81页/共107页第八十一页，共108页。利用均差的牛顿插值法函数f的零阶均差定义为，一阶定义均差为一般地，函数f的k阶均差定义为：利用均差的牛顿插值法多项式为

第82页/共107页第八十二页，共108页。系数的计算过程如表所示表

均差计算表格一阶均差二阶均差三阶均差……n阶均差……………第83页/共107页第八十三页，共108页。在MATLAB中编程实现利用均差牛顿插值法函数为：Newton。功能：求已知数据点的均差形式的牛顿插值多项式；调用格式：f=

Newton(x,y)或f=

Newton(x,y,x0)。其中，x为已知数据点的x坐标向量；

y为已知数据点的y坐标向量；

x0为插值点的x坐标；

f为求得的牛顿插值法多项式或x0处的插值。第84页/共107页第八十四页，共108页。functionf=Newton(x,y,x0)%求已知数据点的均差形式牛顿插值多项式%已知数据点的x坐标向量:x%已知数据点的y坐标向量:y%为插值点的x坐标:x0%求得的均差形式牛顿插值多项式或x0处的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的维数不相等！');return;end第85页/共107页第八十五页，共108页。f=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));endc(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);elsef=collect(f);%将插值多项式展开

f=vpa(f,6);endendend第86页/共107页第八十六页，共108页。例根据下表的数据点求出其均差形式牛顿插值多项式，并计算当x=2.0时y的值。解：>>x=[11.21.82.54];>>y=[11.443.246.2516];>>f=Newton(x,y)f=.182711e-14-.482154e-14*t+1.00000*t^2-.169177e-14*t^3+.211471e-15*t^4>>f=Newton(x,y,2.0)f=4x11.21.82.54y11.443.246.2516第87页/共107页第八十七页，共108页。利用差分的牛顿插值法差分分为向前差分、后向差分和中心差分三种，它们的记法及定义如下为：

其中：代表向前差分；代表向后差分；代表向后差分。第88页/共107页第八十八页，共108页。

假设。为了方便，可构造如表所示的差分表（）。表

差分计算表格………………第89页/共107页第八十九页，共108页。向前牛顿插值向前牛顿插值多项式可表示如下：其中叫做步长，，且的取值范围为。第90页/共107页第九十页，共108页。在MATLAB中编程实现向前牛顿插值法函数为：Newtonforward。功能：求已知数据点的向前牛顿插值法多项式；调用格式：f=

Newtonforward(x,y)或

f=

Newtonforward

(x,y,x0)。其中，x为已知数据点的x坐标向量；

y为已知数据点的y坐标向量；

x0为插值点的x坐标；

f为求得的向前牛顿插值法多项式或x0处的插值。第91页/共107页第九十一页，共108页。functionf=Newtonforward(x,y,x0)%求已知数据点的向前差分牛顿插值多项式%已知数据点的x坐标向量:x%已知数据点的y坐标向量:y%为插值点的x坐标:x0%求得的向前差分牛顿插值多项式或x0处的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的维数不相等！');return;end第92页/共107页第九十二页，共108页。f=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);endc(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)));elsef=collect(f);f=vpa(f,6);endendend第93页/共107页第九十三页，共108页。向前牛顿插值向后牛顿插值多项式可表示如下：其中叫做步长，，且的取值范围为。第94页/共107页第九十四页，共108页。在MATLAB中编程实现向后牛顿插值法函数为：Newtonback。功能：求已知数据点的向前牛顿插值法多项式；调用格式：f=

Newtonback

(x,y)或

f=

Newtonback(x,y,x0)。其中，x为已知数据点的x坐标向量；

y为已知数据点的y坐标向量；

x0为插值点的x坐标；

f为求得的向前牛顿插值法多项式或x0处的插值。第95页/共107页第九十五页，共108页。functionf=Newtonback(x,y,x0)%求已知数据点的向后差分牛顿插值多项式%已知数据点的x坐标向量:x%已知数据点的y坐标向量:y%为插值点的x坐标:x0%求得的向前差分牛顿插值多项式或x0处的插值:fsymst;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('x和y的维数不相等！');return;end第96页/共107页第九十六页，共108页。f=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);endc(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)));elsef=collect(f);f=vpa(f,6);endendend第97页/共107页第九十七页，共108页。例5根据下表的数据点求出其差分形式的牛顿插值多项式，并计算当x=1.55时y的值。解：>>x=[11.21.41.61.8];>>y=[0.84150.9320 0.98540.99960.9738];>>f=Newtonforward(x,y)f=.841500+.108025*t-.169042e-1*t^2-.675000e-3*t^3+.541667e-4*t^4>>f=Newtonforward(x,y,1.55)f=0.9998f=Newtonback(x,y)f=.973800-.457417e-1*t-.198042e-1*t^2+.191667e-3*t^3+.541667e-4*t^4>>f=Newtonback(x,y,1.55)f=0.9998x11.21.41.61.8y0.84150.93200.98540.99960.9738第98页/共107页第九十八页，共108页。Hermite插值Hermite插值满足在节点上等于给定函数值，而且在节点上的导数值也等于给定的导数值。对于高阶导数的情况，Hermite插值多项式比较复杂，在实际中，常常遇到的是函数值与一阶导数给定的情况。在此情况下，n个节点x1，x2，…，xn的Hermite插值多项式的表达式如下：其中，，，第99页/共107页第九十九页，共108页。在MATLAB中编程实现Hermite插值法函数为：Hermite。功能：求已知数据点的Hermite插值法多项式；调用格式：f=

Hermite

(x,y,y_1)或

f=

Hermite

(x,y,y_1,x0)。其中，x为已知数据点的x坐标向量；

y为已知数据点的y坐标向量；

y_1为已知数据点导数向量；

x0为插值点的x坐标；

f为求得的Hermite插值法多项式或x0处的插值。第100页/共107页第一百页，共108页。functionf=Hermite(x,y,y_1,x0)%求已知数据点的向后差分牛顿插值多项式%已知数据点的x坐标向量:x%已知数据点的y坐标向量:y%已知数据点的导数向量:y_1%求得的Hermite插值多项式或x0处的插值:fsymst;f=0.0;if(length(x)==length(y))if(length(y)==length(y_1))n=length(x);elsedisp('y和y的导数的维数不相等！');return;endelsedisp('x和y的维数不相等！');return;end第101页/共107页第一百零一页，共108页。fori=1:nh=1.0;a=0.0;forj=1:nif(j~=i)h=h*(t-x(j))^2/((x(i)-x(j))^2);a=a+1/(x(i)-x(j));endend

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);elsef=vpa(f,6);endendend第102页/共107页第一百零二页，共108页。例6根据