用MATLAB进行数据插值

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

三、插值与拟合的区别和联系第8页,共107页，2024年2月25日，星期天

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

节点可视为由产生,，表达式复杂,，或无封闭形式,，或未知.。

第12页,共107页，2024年2月25日，星期天

构造一个(相对简单的)函数通过全部节点,即再用计算插值，即

返回第13页,共107页，2024年2月25日，星期天

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

已知函数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页，2024年2月25日，星期天拉格朗日(Lagrange)插值特别地:两点一次(线性)插值多项式:三点二次(抛物)插值多项式:第15页,共107页，2024年2月25日，星期天

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

采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.例返回ToMatlablch(larg1)第16页,共107页，2024年2月25日，星期天第17页,共107页，2024年2月25日，星期天分段线性插值n越大，误差越小.

xjxj-1xj+1x0xnxoy第18页,共107页，2024年2月25日，星期天ToMATLABxch11，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页，2024年2月25日，星期天第20页,共107页，2024年2月25日，星期天比分段线性插值更光滑。

xyxi-1xiab

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

三次样条插值g(x)为被插值函数。第22页,共107页，2024年2月25日，星期天例用三次样条插值选取11个基点计算插值(ych)返回ToMATLABych第23页,共107页，2024年2月25日，星期天用MATLAB作插值计算一维插值函数：yi=interp1(x，y，xi，'method')插值方法被插值点插值节点xi处的插值结果‘nearest’

：最邻近插值‘linear’

：线性插值；‘spline’

：三次样条插值；‘cubic’

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

例：在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’)第28页,共107页，2024年2月25日，星期天第29页,共107页，2024年2月25日，星期天

xy机翼下轮廓线例已知飞机下轮廓线上数据如下，求x每改变0.1时的y值。ToMATLAB(plane)返回第30页,共107页，2024年2月25日，星期天第31页,共107页，2024年2月25日，星期天二维插值的定义

xyO第一种（网格节点）：第32页,共107页，2024年2月25日，星期天

已知m

n个节点其中互不相同，不妨设

构造一个二元函数通过全部已知节点,即再用计算插值，即第33页,共107页，2024年2月25日，星期天第二种（散乱节点）：

yx0第34页,共107页，2024年2月25日，星期天已知n个节点其中互不相同，

构造一个二元函数通过全部已知节点,即再用计算插值，即返回第35页,共107页，2024年2月25日，星期天

注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值x

y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O

二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。返回第36页,共107页，2024年2月25日，星期天

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

分片线性插值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第37页,共107页，2024年2月25日，星期天插值函数为：第二片(上三角形区域)：(x,y)满足插值函数为：注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足返回第38页,共107页，2024年2月25日，星期天

双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值x

y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O返回第39页,共107页，2024年2月25日，星期天

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

最邻近插值‘linear’

双线性插值‘cubic’

双三次插值缺省时,双线性插值第40页,共107页，2024年2月25日，星期天例：测得平板表面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个单位的地方进行插值.第41页,共107页，2024年2月25日，星期天第42页,共107页，2024年2月25日，星期天再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)画出插值后的温度分布曲面图.ToMATLAB(wendu)第43页,共107页，2024年2月25日，星期天第44页,共107页，2024年2月25日，星期天例>>[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页，2024年2月25日，星期天选较密的插值点，用默认的线性插值算法进行插值>>[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);>>z0=interp2(x,y,z,x1,y1);>>surf(x1,y1,z0)第46页,共107页，2024年2月25日，星期天立方和样条插值：>>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页，2024年2月25日，星期天算法误差比较>>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页，2024年2月25日，星期天山区地形地貌图已知某处山区地形选点测量坐标数据为：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第49页,共107页，2024年2月25日，星期天山区地形地貌图程序原始地貌图程序：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页，2024年2月25日，星期天山区地形地貌图结果第51页,共107页，2024年2月25日，星期天

通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。ToMATLAB(moutain)返回第52页,共107页，2024年2月25日，星期天第53页,共107页，2024年2月25日，星期天第54页,共107页，2024年2月25日，星期天第55页,共107页，2024年2月25日，星期天第56页,共107页，2024年2月25日，星期天第57页,共107页，2024年2月25日，星期天

插值函数griddata格式为:

cz

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

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

最邻近插值‘linear’

双线性插值‘cubic’

双三次插值'v4'-Matlab提供的插值方法缺省时,双线性插值第58页,共107页，2024年2月25日，星期天>>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页，2024年2月25日，星期天第60页,共107页，2024年2月25日，星期天误差分析>>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页，2024年2月25日，星期天第62页,共107页，2024年2月25日，星期天例：

在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页，2024年2月25日，星期天去除在(-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)第64页,共107页，2024年2月25日，星期天新样本点拟合曲面>>[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页，2024年2月25日，星期天误差分析>>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页，2024年2月25日，星期天命令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()函数一致。第67页,共107页，2024年2月25日，星期天例：通过函数生成一些网格型样本点，试根据样本点进行拟合，并给出拟合误差。>>[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页，2024年2月25日，星期天slice(x0,y0,z0,V1,[-0.5,0.3,0.9],[0.6,-0.1],[-1,-0.5,0.5,1])title('SlivesforFourDimFigures')第69页,共107页，2024年2月25日，星期天[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页，2024年2月25日，星期天

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

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)第75页,共107页，2024年2月25日，星期天海底曲面图第76页,共107页，2024年2月25日，星期天%程序二：插值并作出水深小于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页，2024年2月25日，星期天水深小于5的海域范围第78页,共107页，2024年2月25日，星期天拉格朗日插值法拉格朗日插值法是基于基函数的插值方法，插值多项式可表示为其中称为i次基函数：第79页,共107页，2024年2月25日，星期天在MATLAB中编程实现拉格朗日插值法函数为：Language。功能：求已知数据点的拉格朗日多项式；调用格式：f=

Language(x,y)或f=

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

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

x0为插值点的x坐标；

f为求得的拉格朗日多项式或x0处的插值。第80页,共107页，2024年2月25日，星期天

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页，2024年2月25日，星期天例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第82页,共107页，2024年2月25日，星期天利用均差的牛顿插值法函数f的零阶均差定义为，一阶定义均差为一般地，函数f的k阶均差定义为：利用均差的牛顿插值法多项式为

第83页,共107页，2024年2月25日，星期天系数的计算过程如表所示表

均差计算表格一阶均差二阶均差三阶均差……n阶均差……………第84页,共107页，2024年2月25日，星期天在MATLAB中编程实现利用均差牛顿插值法函数为：Newton。功能：求已知数据点的均差形式的牛顿插值多项式；调用格式：f=

Newton(x,y)或f=

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

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

x0为插值点的x坐标；

f为求得的牛顿插值法多项式或x0处的插值。第85页,共107页，2024年2月25日，星期天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第86页,共107页，2024年2月25日，星期天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第87页,共107页，2024年2月25日，星期天例根据下表的数据点求出其均差形式牛顿插值多项式，并计算当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第88页,共107页，2024年2月25日，星期天利用差分的牛顿插值法差分分为向前差分、后向差分和中心差分三种，它们的记法及定义如下为：

其中：代表向前差分；代表向后差分；代表向后差分。第89页,共107页，2024年2月25日，星期天

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

差分计算表格………………第90页,共107页，2024年2月25日，星期天向前牛顿插值向前牛顿插值多项式可表示如下：其中叫做步长，，且的取值范围为。第91页,共107页，2024年2月25日，星期天在MATLAB中编程实现向前牛顿插值法函数为：Newtonforward。功能：求已知数据点的向前牛顿插值法多项式；调用格式：f=

Newtonforward(x,y)或

f=

Newtonforward

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

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

x0为插值点的x坐标；

f为求得的向前牛顿插值法多项式或x0处的插值。第92页,共107页，2024年2月25日，星期天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第93页,共107页，2024年2月25日，星期天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第94页,共107页，2024年2月25日，星期天向前牛顿插值向后牛顿插值多项式可表示如下：其中叫做步长，，且的取值范围为。第95页,共107页，2024年2月25日，星期天在MATLAB中编程实现向后牛顿插值法函数为：Newtonback。功能：求已知数据点的向前牛顿插值法多项式；调用格式：f=

Newtonback

(x,y)或

f=

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

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

x0为插值点的x坐标；

f为求得的向前牛顿插值法多项式或x0处的插值。第96页,共107页，2024年2月25日，星期天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第97页,共107页，2024年2月25日，星期天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第98页,共107页，2024年2月25日，星期天例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第99页,共107页，2024年2月25日，星期天Hermite插值Hermite插值满足在节点上等于给定函数值，而且在节点上的导数值也等于给定的导数值。对于高阶导数的情况，Hermite插值多项式比较复杂，在实际中，常常遇到的是函数值与一阶导数给定的情况。在此情况下，n个节点x1，x2，…，xn的Hermite插值多项式的表达式如下：其中，，，第100页,共107页，2024年2月25日，星期天在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处的插值。第101页,共107页，2024年2月25日，星期天functionf=Hermite(x,y,y_1,x0)%求已知数据点的向后差分牛顿插值多项式%已知数据点的x坐标向量:x%已知数据点的y坐标向量:y%已知数据点的导数向量:y_1%求得的Hermite插值多项式或x0处的插值:fsymst;f=0.0;if(le