数学建模插值与拟合

插值与拟合一、插值的基本原理二、拟合的基本原理三、插值与拟合的关系四、插值的MATLAB实现

我们经常会遇到大量的数据需要处理，而处理数据的关键就在于这些算法，例如数据拟合、参数估计、插值等数据处理算法。此类问题在MATLAB中有很多现成的函数可以调用，熟悉MATLAB，这些方法都能游刃有余的用好。一、概述

数据拟合在很多赛题中有应用，与图形处理有关的问题很多与插值和拟合有关系，例如98年美国赛A题，生物组织切片的三维插值处理，94年A题逢山开路，山体海拔高度的插值计算，2003年吵的沸沸扬扬的“非典”问题也要用到数据拟合算法，观察数据的走向进行处理，2005年的雨量预报的评价的插值计算。2001年的公交车调度拟合问题，2003年的饮酒驾车拟合问题。

在实际中，常常要处理由实验或测量所得到的一些离散数据。插值与拟合方法就是要通过这些数据去确定某一类已知函数的参数或寻求某个近似函数，使所得到的近似函数与已知数据有较高的拟合精度。如果要求这个近似函数（曲线或曲面）经过所已知的所有数据点，则称此类问题为插值问题。二、基本概念

如果不要求近似函数通过所有数据点，而是要求它能较好地反映数据变化规律的近似函数的方法称为数据拟合。近似函数不一定（曲线或曲面）通过所有的数据点。++++++++++f=a1+a2x+a3x2f=a1+a2x+a3x21、联系都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数的方法。2、区别插值问题不一定得到近似函数的表达形式，仅通过插值方法找到未知点对应的值。数据拟合要求得到一个具体的近似函数的表达式。

三、插值与拟合的区别和联系四、插值的使用及求解

当数据量不够，需要补充，且认定已有数据可信时,通常利用函数插值方法。实际问题当中碰到的函数f(x)

是各种各样的，有的表达式很复杂，有的甚至给不出数学的式子，只提供了一些离散数据，警如，某些点上的函数值和导数值。4.1引言

选用不同类型的插值函数，逼近的效果就不同，一般有：（1）拉格朗日插值（lagrange插值）（2）分段线性插值（3）Hermite（4）三次样条插值。4.2插值方法

插值概念与基础理论1.1.1插值问题的提法对于给定的函数表xx0x1…….xnY=f(x)y0y1……..yn(其中在[a,b]上连续，x0，

x1，…，xn

是[a,b]上的n+1个互异的点)，在某函数类{(x)

}中求一个函数(x)

，使(xi)=yi,(i=0,1,2,…,n)(2)（1）并用函数(x)

作为函数

y=f(x)的近似函数，即

y=f(x)(x),(x∈[a,b])

这类问题称为插值问题。[a,b]称为插值区间，x0,x1,...,xn

称为插值节点，（2）称为插值条件，插值条件是选择近似函数的标准，满足此条件的近似函数

(x)称为插值函数，f(x)称为被插值函数。

函数类{(x)

}有多种取法，常用的有代数多项式、三角函数和有理函数。

最简单的插值函数是代数多项式，相应的插值问题称为多项式插值。

最简单的插值函数是代数多项式，相应的插值问题称为多项式插值。

求的n次插值多项式的几何意义，就是上的若干个节点，作一条代数曲线来近似代替曲线。如图所示。

通过曲线§1.2插值多项式的求法

在前面讨论插值多项式的存在唯一性时，实际上已提供了它的一种求法，即通过求解线性方程组来确定其系数ai(i=0,1,2,…,n)

但是这种方法不仅计算量大，而且因不能获得简明的表达式而给理论和应用研究带来不便。在这里我们学习两种简便而实用的求答。1.2.1拉格朗日插值多项式

在线性代数中知道，所有次数不超过n次的多项式构成一个n+1维线性空间。其基有各种不同的取法。因此尽管满足条件（4）的n次插值多项式是唯一的，然而它的表达式可以有多种不同的形式。如果取满足条件：的一组n次多项式作为上述线性空间的基，则容易看出因此，由n+1个代数多项式线性生成的多项式（10）就是满足插值条件的n次插值多项式。（10）（9）满足条件（9）的多项式称为n+1个节点的n次基本插值多项式（或n次基函数）

显然，求拉格朗日多项式的关键是求n次插值基函数。因此，可设因为为n次多项式，且两种特殊的Lagrange插值多项式1.线性插值(两点插值)

最简单的插值是线性插值(此时n=1),这时插值问题就是求一次多项式P1(x)=a0+a1x

使它满足条件P1(x0)=y0,P1(x1)=y1,这时于是线性插值多项式为即它就是通过M0(x0,y0)和M1(x1,y1)两点的线段.2.抛物插值

线性插值仅仅用两个节点以上的信息，精确度较差。为了提高精确度，我们进一步考察以下三点的插值问题(n=2)：这时由此得到抛物插值多项式抛物插值又称三点插值.

例1

已知的函数表10111213142.30262.39792.48492.56492.6391并估计误差。分别用拉格朗日线性和抛物线插值求的近似值，%lagrange插值法的程序functiony=lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endclearx0=[1011121314];y0=[2.30262.3979,2.4849,2.56492.6391];x=10:0.1:15;y=lagrange(x0,y0,x);plot(x0,y0,’+’)Holdonplot(x,y,’-r’)10111213142.30262.39792.48492.56492.639110111213142.30262.39792.48492.56492.6391lnx.m

1901年龙格(Runge)给出一个例子:

定义在区间[-1，1]上，这是一个光滑函数，它的任意阶导数都存在，对它在[-1，1]上作等距节点插值时，插值多项式情况，见图:从图中，可见，在靠近-1或1时，余项会随n值增大而增大，如P12(0.96)=3×6!但f(0.96)=0.25

从图中，还可发现，在0附近插值效果是好的，即余项较小，另一种现象是插值多项式随节点增多而振动更多。这种插值多项式当节点增加时反而不能更好地接近被插之数的现象，称为龙格现象。

上述现象和定理，告诉我们用高次插值多项式是不妥当的，从数值计算上可解释为高次插值多项式的计算会带来舍入误差的增大，从而引起计算失真。那么如何提高插值精度呢？采用分段插值是一种办法。实践上作插值时一般只用一次、二次最多用三次插值多项式。分段线性插值的构造：

设f(x)是定义在[a,b]上的函数，在[a,b]上节点

a=x0<x1<x2<…<xn-1<xn=b,

的函数值为y0,y1,y2,…yn-1,yn。

(x)在每个子区间[xi,xi+1](i=0,1,2,,n-1)上是一次插值多项式;这种分段低次插值称为分段线性插值.在几何上就是用折线段带代替曲线,故分段线性插值又称为折线插值.1.2.2分段线性插值分段线性插值曲线图：曲线的光滑性较差在节点处有尖点但如果增加节点的数量减小步长,会改善插值效果分段二次插值即：选取跟节点x最近的三个节点xi-1,xi,xi+1进行二次插值，即在区间[xi-1,xi+1]，取：这种分段的低次插值叫分段二次插值，在几何上就是用分段抛物线代替y=f(x)，故分段二次插值又和分段抛物插值。

在许多实际问题中，不仅要求插值函数与被插值函数在节点处的函数值相同，而且还要求插值函数与被插值函数在某些节点处的导数值、甚至高阶导数值也相同。Hermite插值Hermite插值多项式Hermite插值基函数又由则可令常数2.再确定根据Lagrange基函数的定义，由1.先确定求Hermite插值基函数计算出基函数后，可组合得出Hermite插值多项式：解例3.9

应用最广泛的是三次Hermite插值多项式，即n=1时的情况。什么是样条:是指飞机或轮船等的制造过程中为描绘出光滑的外形曲线(放样)所用的工具样条本质上是一段一段的三次多项式拼合而成的曲线在拼接处,不仅函数是连续的,且一阶和二阶导数也是连续的1946年,Schoenberg将样条引入数学,即所谓的样条函数

三次样条插值------(1)定义1.三次样条插值函数xi0123yi00.521.5clearx0=[0123];y0=[00.521.5];x=0:0.1:3;pp1=csape(x0,y0,’complete’);y3=ppval(pp1,x);%计算插值函数在x处的值plot(x0,y0,’+’,x,y3,’r’)sx.mxi1245yi1342sxsec.mMatlab实现：实现分段线性插值不需要编制函数程序，它自身提供了内部的功能函数interp1(一维插值)intep2(二维)interp3(三维)intern(n维)4.3MATLAB实现插值用MATLAB作插值计算一维插值函数：yi=interp1(x，y，xi，'method')插值方法被插值点插值节点xi处的插值结果‘nearest’

最邻近插值；‘linear’

线性插值；‘spline’

三次样条插值；‘cubic’

立方插值；缺省时分段线性插值．

注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围．

例1

已知的函数表10111213142.30262.39792.48492.56492.6391分别用拉格朗日线性和抛物线插值求的近似值，clearx0=[1011121314];y0=[2.30262.3979,2.48492.56492.6391];x=10:0.1:15;y1=interp1(x0,y0,x,’linear’);yy1=interp1(x0,y0,11.5,linear’);y2=interp1(x0,y0,x,’cubic’);yy2=interp1(x0,y0,11.5,’cubic’);subplot(1,2,1)plot(x0,y0,’+’,x,y1,11.5,yy1,’rO’)title(‘Piecewiselinear’)subplot(1,2,2)plot(x0,y0,’+’,x,y2,11.5,yy2,’rO’)title(‘Piecewisecubic’)fenduanchazhi.m例：从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,'+',h,t,hours,temps,'r:')%作图xlabel('Hour'),ylabel('DegreesCelsius’)xy机翼下轮廓线例已知飞机下轮廓线上数据如下，求x每改变0.1时的y值．plane.m

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

最邻近插值；‘linear’

双线性插值；‘cubic’

双三次插值；缺省时双线性插值.二维插值函数：例：测得平板表面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个单位的地方进行插值.再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)画出插值后的温度分布曲面图.ToMATLAB(wendu)

通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较．ToMATLAB(moutain)

插值函数griddata格式为:

cz

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

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

双线性插值‘cubic’

双三次插值'v4'-MATLAB提供的插值方法缺省时,双线性插值

例在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）×（-50，150）里的哪些地方船要避免进入．ToMATLABhd14.作出水深小于5的海域范围,即z=5的等高线.2.在矩形区域(75,200)×(-50,150)进行插值。1.输入插值基点数据3.作海底曲面图%程序一：插值并作海底曲面图

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)海底曲面图%程序二：插值并作出水深小于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)

水深小于5的海域范围估计水塔的流量2、解题思路3、算法设计与编程1、问题

某居民区有一供居民用水的园柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量．通常水泵每天供水一两次，每次约两小时.水塔是一个高12.2米，直径17.4米的正园柱．按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作．表1是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量．流量估计的解题思路拟合水位~时间函数确定流量~时间函数估计一天总用水量

拟合水位~时间函数

测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段），和3个水泵不工作时段（以下称第1时段t=0到t=8.97，第2次时段t=10.95到t=20.84和第3时段t=23以后）．对第1、2时段的测量数据直接分别作多项式拟合，得到水位函数．为使拟合曲线比较光滑，多项式次数不要太高，一般在3~6．由于第3时段只有3个测量记录，无法对这一时段的水位作出较好的拟合．2、确定流量~时间函数

对于第1、2时段只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内．3、一天总用水量的估计

总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到。算法设计与编程1、拟合第1、2时段的水位，并导出流量2、拟合供水时段的流量3、估计一天总用水量4、流量及总用水量的检验1、拟合第1时段的水位，并导出流量设t，h为已输入的时刻和水位测量记录（水泵启动的4个时刻不输入），第1时段各时刻的流量可如下得：1）c1=polyfit（t（1：10），h（1：10），3）；

%用3次多项式拟合第1时段水位，c1输出3次多项式的系数2）a1=polyder（c1）；

%a1输出多项式（系数为c1）导数的系数

3）tp1=0：0.1：9；

x1=-polyval（a1，tp1）；

%x1输出多项式（系数为a1）在tp1点的函数值（取负后边为正值），即tp1时刻的流量

4）流量函数为：2、拟合第2时段的水位，并导出流量设t，h为已输入的时刻和水位测量记录（水泵启动的4个时刻不输入），第2时段各时刻的流量可如下得：1）c2=polyfit(t(10.9:21),h(10.9:21),3)；

%用3次多项式拟合第2时段水位，c2输出3次多项式的系数2）a2=polyder(c2)；

%a2输出多项式（系数为c2）导数的系数

3）tp2=10.9:0.1:21;x2=-polyval(a2,tp2);%x2输出多项式（系数为a2）在tp2点的函数值（取负后边为正值），即tp2时刻的流量4）流量函数为：

3、拟合供水时段的流量在第1供水时段（t=9~11）之前（即第1时段）和之后（即第2时段）各取几点，其流量已经得到，用它们拟合第1供水时段的流量．为使流量函数在t=9和t=11连续，我们简单地只取4个点，拟合3次多项式（即曲线必过这4个点），实现如下：

xx1=-polyval（a1，[89]）；%取第1时段在t=8,9的流量

xx2=-polyval（a2，[1112]）；%取第2时段在t=11,12的流量

xx12=[xx1xx2]；

c12=polyfit（[891112]，xx12，3）；%拟合3次多项式

tp12=9：0.1：11；

x12=polyval（c12，tp12）；%x12输出第1供水时段各时刻的流量拟合的流量函数为：

在第2供水时段之前取t=20，20.8两点的流水量，在该时刻之后（第3时段）仅有3个水位记录，我们用差分得到流量，然后用这4个数值拟合第2供水时段的流量如下：

dt3=diff（t(22：24)）；%最后3个时刻的两两之差

dh3=diff（h(22：24)）；%最后3个水位的两两之差

dht3=-dh3./dt3；%t(22)和t(23)的流量

t3=[2020.8t(22)t(23)]；

xx3=[-polyval(a2，t3(1：2)，dht3)]；%取t3各时刻的流量

c3=polyfit（t3，xx3，3）；%拟合3次多项式

t3=20.8：0.1：24；

x3=polyval（c3，tp3）；%x3输出第2供水时段（外推至t=24）各时刻的流量拟合的流量函数为：3、一天总用水量