数据插值与曲线拟合

大多数数学建模问题都是从实际工程或生活中提炼出来的，往往带有大量的离散的实验观测数据，要对这类问题进行建模求解，就必须对这些数据进行处理。其目的是为了从大量的数据中寻找它们反映出来的规律。用数学语言来讲，就是要找出与这些数据相应的变量之间的近似关系。对于非确定性关系，一般用统计分析的方法来研究，如回归分析的方法。对于确定性的关系，即变量间的函数关系，一般可用数据插值与拟合的方法来研究。本讲学习数据插值与拟和的基本方法和相关的MATLAB命令。1引例拟合并不要求函数图像通过这些点，但要求在某种准则下，该函数在这些点处的函数值与给定的这些值能最接近。简单地讲，插值是对于给定的n组离散数据，寻找一个函数，使该函数的图像能严格通过这些数据对应的点。例1：对于下面给定的4组数据，求在x=175处y的值。x144169225y121315例1：对于下面给定的4组数据，求在x=175处y的值。x144169225y121315这就是一个插值问题。利用所得的函数来求x=175处y的值。我们可以先确定插值函数，再需要说明的是这3组数据事实上已经反映出x与y的的函数关系为：关系是不明显的。，当数据量较大时，这种函数也就是说，插值方法在处理数据时，不论数据本身对应的被插值函数是否已知，它都要找到一个通过这些点的插值函数，此函数是被插值函数的一个近似，从而通过插值函数来计算被插值函数在未知点处的近似值。对于所构造的插值函数要求相对简单，便于计算，一般选用多项式函数来逼近。例2：观测物体的直线运动，得以下数据，求物体的运动方程。t（秒）00.91.93.03.95.0s（米）010305080110例2：观测物体的直线运动，得以下数据，求物体的运动方程。t（秒）00.91.93.03.95.0s（米）010305080110这是一个拟合问题，其明显的特征是与数据对应的函数未知，要找到一个函数来比较准确地表述这些数据蕴藏的规律。显然，我们找出的函数不一定会通过这些点，也没有必要，因为观测数据本身并不是完全准确的。2数据插值的基本原理一般地，对于给定的n+1组数据互不相等，确定一个n次多项式使。其中称为插值函数，为插值节点，区间，为插值称为插值条件。当n=1时为线性插值。表示过两点的直线方程，即定理：满足n+1个插值节点的次数不超过n次的多多项式存在且唯一。稍加整理，即得记则它们满足：称为基函数，那么是两个基函数的线性组合,也称为Lagrange线性插值函数。当n=2时为抛物插值。表示过三点的抛物线方程，使它们满足则可表示为三个基函数的线性组合，即仿照线性插值的情形取基函数也称为Lagrange抛物插值函数。一般地，满足插值条件的n次多项式为：其中基函数满足上述多项式插值又称为n次Lagrange插值。说明：1、多项式插值的基函数仅与节点有关，而与被插值的原函数无关；2、插值多项式仅由数对确定，而与数对的排列次序无关。3、多项式插值除拉格朗日多项式插值法外，还有牛顿(Newton)插值法、埃尔米特(Hermite)插值法、三次样条插值法等，可参看有关数值分析的书籍。其中Newton插值是拉格朗日插值的一种等价变形，Hermite插值一种带导数插值条件的插值。例将[0,/2]n等分，用

g(x)=cos(x)产生n+1个节点，作Pn(x)（取n=1,2)，计算cos(/6)。解:n=1,(x0,y0)=(0,1),(x1,y1)=(/2,0),

P1(x)=1-2x/,cos(/6)=P1(/6)≈0.6667

n=2,(x0,y0)=(0,1),(x1,y1)=(/4,0.7071),(x2,y2)=(/2,0),P2(x)=8(x-/4)(x-/2)/2-16x(x-/2)0.7071/2cos(/6)=P2(/6)≈0.8508精确值：cos(/6)≈0.8660下面来求解引例1(课堂练习)。引例1：对于下面给定的4组数据，求在x=175处y的值。x144169225y121315解：用一次拉格朗日插值：所以取为插值节点，则计算得因为插值点位于和之间，，于是用二次拉格朗日插值：取，则计算得，于是（的准确值为）由上例看出，二次插值的精度明显要比一次插值要高。但对于拉格朗日多项式插值，是否插值其精度就一定越高呢？答案是：对于某些函数，适当地提高插值多项式的次数，会提高计算精度。但与此同时，多项式的次数增大可能造成插值函数的收敛性和稳定性越来越差，逼近的效果往往不理想，一个典型的例子是函数选取不同插值节点个数n+1，其中n为插值多项式的次数，使得它在结点的值与被插函数在对应结点的值相等。当n分别取2,4,6,8,10时，绘出的插值图形如下。从图中可看出，图形显示出振荡现象，在5和-5附近误差很大。这种现象叫做Runge现象。这说明，在大范围内使用高次插值，逼近效果往往并不理想。解决此问题的思路是化整为零，采用分段插值，即在小范围内使用低次多项式插值。不是去寻求整个插值区间上的一个高次多项式，也就是说插值区间划分为若干个小区间，在每个小区间上用低而是把次多项式插值，在整个插值区间上就得到一个分段插值函数。区间的划分可以是任意的，各个区间上插值多项式的次数的选取也可按具体问题选择。在分段插值中，较为简单的是分段线性插值。实际数学建模中，在光滑性要求不高的条件下，分段线性或二次插值基本可以满足需要。问题中提出的插值问题，有一些插值函数曲线要求然而实际具有较高的光滑性，如飞机机翼的下轮廓线。分段线性插值虽然简单，但插值函数在结点处的一阶导数一般不存在，光滑性不高，样条插值的提出。这就导致了三次在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。分段多项式达到较高阶光滑性的方法？是否存在较低次的就是一个很好的例子。三次样条插值3三次样条插值三次样条插值是一种非常有效的插值方法，它在实际工程中有着非常重要的应用。三次样条插值的理论推导是比较复杂的，但在数学软件MATLAB中有现成的调用程序，这样我们就可直接借助计算机来进行运算。下面简单介绍一下三次样条插值的基本原理。定义：设给定区间上的一个划分如果函数满足条件：定义：设给定区间上的一个划分如果函数满足条件：（1）在每个子区间是三次多项式；（2）在区间上连续，记作（3）对于在节点上给定的函数值满足则称为在区间上的三次样条插值函数。简单地说，已经知道函数在节点上的函数值多项式函数，现要求一个三次，使满足且。由定义可知，是区间上的分段三次插值多项式，即由于，这个函数的曲线具有二阶光滑度，看起来就很光顺了，能满足一般工程上的需要。其中是子区间插值于两点的三次多项式，即下面简单介绍一下三次样条插值函数的推导。现要求为待定系数，共4n个。已知条件：1）共n+1个方程；2）共3(n-1)

个方程。现要求4n个待定系数，但只有(n+1)+3(n-1)=4n-2个方程，故需要补充两个方程，即所谓的边界条件。通常有以下三类边界条件：3.1）给定两个端点处的导数，即3.2）给定两个端点处的导数即3.3）周期性条件，即4用MATLAB软件求解插值问题在MATLAB中提供了一个一维插值函数interp1，它的调用格式为cy=interp1(x,y,cx,‘method’)其中x、y是所给数据的横纵坐标，要求x的分量按升序或降序排列，cx是待求的插值点的横坐标，返回值cy是待求的插值点的纵坐标，method是插值方法，该函数提供了四种可选的插值方法：nearest——最邻近点插值。点和这两已知点间位置的远近来进行插值，取较近已知它根据已知两点间的插值插值点处的函数值作为未知插值点处的函数值。linear——线性插值。它将相邻的数据点用直线相连，按所生成的直线进行插值。spline——三次样条插值。它利用已知数据求出样条函数后，按样条函数进行插值。cubic——三次插值。它利用已知数据求出三次多项式函数后，按三次多项式函数进行插值。缺省时插值方法为分段线性插值。下面用该函数来求解下列插值问题。对于下面给定的4组数据，求在x=110处y的值。x100121144169y10111213输入命令：x=[100121144169];y=[10111213];cx=110;cy=interp1(x,y,cx,'linear');运行结果为cy=10.4762。由于线性插值只需要两个点，因而在上述命令中实际上只用了前两个点。若将最后一个命令中的method改为缺省、nearest、cubic和spline，运行结果为依次为cy=10.4762、cy=10、cy=10.4869、cy=10.4877通过比较，显然三次样条插值的结果最好。例：在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。程序：hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');%(直接输出数据将是很多的)plot(hours,temps,'+',h,t,'r:')%作图xlabel('Hour'),ylabel('DegreesCelsius')

xy机翼下轮廓线例已知飞机下轮廓线上数据如下，求X每改变0.1时的Y值。程序：lch(lagr1)functiony=lagr1(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;endx0=[035791112131415];y0=[01.21.72.02.12.01.81.21.01.6];x=0:0.1:15;y1=lagr1(x0,y0,x);y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,'spline');subplot(3,1,1)plot(x0,y0,'k+',x,y1,'r')gridtitle('lagrange')subplot(3,1,2)plot(x0,y0,'k+',x,y2,'r')gridtitle('piecewiselinear')subplot(3,1,3)plot(x0,y0,'k+',x,y3,'r')gridtitle('spline')曲线拟合是指：已知平面上n个点（xi,yi)i=1,…n,寻求一个函数（曲线）y=f(x),使f(x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。

+++++++++xyy=f(x)(xi,yi)

i

i

为点（xi,yi)与曲线y=f(x)的距离5曲线拟合的基本原理拟合与插值的区别函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。问题：给定一批数据点，需确定满足特定要求的曲线或曲面。解决方案：若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是曲线数据拟合，又称曲线拟合或曲面拟合。若要求所求曲线（面）通过所给所有数据点，就是插值问题；根据曲线拟合问题的定义，其关键在于准则的选取，选取的准则不同，其对应的拟合方法及其复杂程度也不相同。对于一维曲线拟合，设n个不同的离散数据点为，要寻找的拟合曲线方程为记拟合函数在处的偏差为常用的准则有：准则1：选取，使所有偏差的绝对值之和最小，即准则2：选取，使所有偏差的绝对值的最大值最小，即准则3：选取，使所有偏差的平方和最小，即相对而言，准则3最便于计算，因而通常根据准则3来选取拟合曲线。准则3又称为最小二乘准则，对应的曲线拟合方法称为最小二乘法。线性最小二乘法的基本思路第一步:先选定一组函数

r1(x),r2(x),…rm(x),m<n,令

f(x)=a1r1(x)+a2r2(x)+…+amrm(x)其中

a1,a2,…am

为待定系数。第二步:确定a1,a2,…am

的最小二乘准则：使n个点（xi,yi)与曲线y=P(x)的距离

i的平方和最小

。记

问题归结为，求

a1,a2,…am

使

J(a1,a2,…am)最小。线性最小二乘法的求解要使J(a1,a2,…am)最小，的必要条件得则由多元函数取得极值即亦即是未知量的线性方程组，称之为正定方程组。选定一组函数组解出后，就可由正规方程，于是就可得线性最小二乘拟合函数所给数据的散点图，观察数据所呈现出来的曲线的大致一般的做法是首先绘出形状，再结合该问题所在专业领域内的相关规律和结论，来确定拟合函数的形式。实际操作时可在直观判断的基础上，选几种常用的曲线分别进行拟合，比较选择拟合效果最好的曲线。面对一组数据，作线性最小二乘拟合时，恰当选定函是一个难点。数常用的曲线有直线、多项式、双曲线和指数曲线等。另外，曲线拟合又可分为线性曲线拟合和非线性曲线拟合。一般地，如果拟合函数中的系数以线性形式出现，全部如拟合函数为线性拟合，也称为多项式拟合；若拟合函数中的系数不能全部以线性形式出现，如指数拟合函数为非线性曲线拟合。实际应用中，多项式最小二乘拟合用的较多，MATLAB中也有专用函数。线性最小二乘拟合f(x)=a1r1(x)+…+amrm(x)中函数{r1(x),…rm(x)}的选取方法1.通过机理分析建立数学模型来确定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx

2.将数据(xi,yi)i=1,…n

作图，通过直观判断确定f(x)：6用MATLAB软件求解拟合问题在MATLAB中提供了一个多项式最小二乘拟合函数polyfit(x,y,n)，它的调用格式为P=polyfit(x,y,n)拟合多项式按自变量降幂排列的系数向量

输入同长度的数组X，Y拟合多项式次数下面用该函数来求解拟合问题引例2：例2：观测物体的直线运动，得以下数据，求物体的运动方程。t（秒）00.91.93.03.95.0s（米）010305080110输入命令：t=[00.91.933.95];s=[010305080110];plot(t,s,'*-')xlabel('运动时间——t（秒）')ylabel('运动位移——s（米）')gtext('物体运动的时间与位移散点图')下面显示的是物体运动的时间与位移散点图：不难看出图形近似为一条直线，因此猜测用一次多项式来拟合，输入命令：

P=polyfit(t,s,1)运行结果为：P=22.2538-7.8550即下面绘出的是拟合曲线和散点图对比图形，可以看出拟合效果并不理想。根据物理学中物体运动的方程，我们用二次曲线来拟合，输入命令：P=polyfit(t,s,2)得到拟合函数为：对比图形如下，可见曲线拟合本身就是一个猜测的过程，通常是不断地修正拟合函数，使拟合效果达到满意的程度。可以看出拟合效果有明显改善，拟合曲线与散点图基本上是吻合的，因此该物体运动的方程是7建模案例(1992年A题：农作物施肥效果分析)某地区作物生长所需要的营养元素主要有氮（N）、钾（K）、磷（P）。某作物研究所在该地区对土豆与生菜做了一定数量的实验，实验数据如下列表格所示，其中ha表示公顷，t表示吨，kg表示公斤。当一个营养素的施肥量变化时，总将另两个营养素的施肥量保持在第七个水平上，于N的施肥量做实验时，P与K的施肥量分别取为如对土豆产量关196kg/ha与372kg/ha。试分析施肥量与产量之间的关系。土豆：

NPK施肥量(kg/ha)产量(t/ha)施肥量(kg/ha)产量(t/ha)施肥量(kg/ha)产量(t/ha)0346710113520225933640447115.1821.3625.7232.2934.0339.4543.1543.4640.8330.7502449739814719624529434233.4632.4736.0637.9641.0440.0941.2642.1740.3642.730479314018627937246555865118.9827.3534.8638.5238.4437.7338.4343.8742.7746.22生菜：NPK施肥量(kg/ha)产量(t/ha)施肥量(kg/ha)产量(t/ha)施肥量(kg/ha)产量(t/ha)028568411216822428033639211.0212.7014.5616.2717.7522.5921.6319.3416.1214.11049981471962943914895876856.399.4812.4614.3817.1021.9422.6421.3422.0724.530479314018627937246555865115.7516.7616.8916.2417.5619.2017.9715.8420.1119.40模型假设：1、研究所的实验是在相同的正常实验条件(如充足的水分供应，正常的耕作程序)下进行的，产量的变化是由施肥量的改变引起的，产量与施肥量之间存在一定的规律。（此假设的目的是抓住影响产量的主要因素而剔除次要因素，使要研究的问题内部诸因素明朗化，即抓住主要矛盾）2、土壤本身已含有一定数量的氮、磷、钾等肥料，即具有一定的天然肥力。2、土壤本身已含有一定数量的氮、磷、钾等肥料，即具有一定的天然肥力。（此假设非常符合常理，而且实验数据也证明了此假设的合理性，因而此假设将实验数据中所隐藏的信息清晰化）3、每次实验是相互独立的，互不影响。（此假设澄清了在连续进行的实验中，后期实验产量与前期施肥无关）符号说明：

：农作物产量；：施肥量；N、K、P

：氮、磷、钾肥的施肥量；：农产品价格；：肥料价格；Tn、Tp、Tk：氮、磷、钾肥的价格；：常数。问题分析：1、普遍规律施肥量与产量满足下图所示关系，它分为三个不同的区段，第二区段，随着施肥量的增加，作物产量平缓上升，一定限度后，第三区段，当施肥量超过产量反而随施肥量的增加而减少。施肥量的增加而迅速增加，在第一区段，当施肥量较小时，作物产量随2、数据分析通过绘制散点图，初步得到农作物产量与施肥量间的定性认识。从散点图可以发现，氮肥施加量与农作物的产量大致呈指数关系，磷肥施加量与农作物产量大致呈分段直线关系