数据插值与曲线拟合.ppt

《数据插值与曲线拟合.ppt》由会员分享，可在线阅读，更多相关《数据插值与曲线拟合.ppt（69页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、数据插值与曲线拟合现在学习的是第1页，共69页1 引例引例拟合并不要求函数图像通过这些点，但要求在某种准则下，该函数在这些点处的函数值与给定的这些值能最接近。简单地讲，插值是对于给定的n组离散数据，寻找一个函数，使该函数的图像能严格通过这些数据对应的点。例1：对于下面给定的4组数据，求在x=175处 y 的值。x144169225y121315现在学习的是第2页，共69页例1：对于下面给定的4组数据，求在x=175处 y 的值。x144169225y121315这就是一个插值问题。利用所得的函数来求x=175处 y 的值。我们可以先确定插值函数，再需要说明的是这3组数据事实上已经反映出 x与y

2、的 的函数关系为：关系是不明显的。，当数据量较大时，这种函数也就是说，插值方法在处理数据时，不论数据本身对应的被插值函数 是否已知，它都要找到一个通过这些点的插值函数，此函数是被 现在学习的是第3页，共69页插值函数的一个近似，从而通过插值函数来计算被插值函数在未知点处的近似值。对于所构造的插值函数要求相对简单，便于计算，一般选用多项式函数来逼近。例2：观测物体的直线运动，得以下数据，求物体的运动方程。t（秒）00.91.93.03.95.0s（米）010305080110现在学习的是第4页，共69页例2：观测物体的直线运动，得以下数据，求物体的运动方程。t（秒）00.91.93.03.95.

3、0s（米）010305080110这是一个拟合问题，其明显的特征是与数据对应的 函数未知，要找到一个函数来比较准确地表述这些数 据蕴藏的规律。显然，我们找出的函数不一定会通过 这些点，也没有必要，因为观测数据本身并不是完全准确的。现在学习的是第5页，共69页2 数据插值的基本原理数据插值的基本原理一般地，对于给定的n+1组数据 互不相等，确定一个n次多项式 使。其中 称为插值函数，为插值节点，区间，为插值称为插值条件。当n=1时为线性插值。表示过两点 的直线方程，即 定理：满足n+1个插值节点的次数不超过n次的多多项式存在且唯一。现在学习的是第6页，共69页稍加整理，即得 记 则它们满足：称

4、为基函数，那么 是两个基函数的线性组合,也称为Lagrange 线性插值函数。现在学习的是第7页，共69页当 n=2 时为抛物插值。表示过三点 的抛物线方程，使它们满足则 可表示为三个基函数的线性组合，即 仿照线性插值的情形取基函数现在学习的是第8页，共69页也称为Lagrange 抛物插值函数。一般地，满足插值条件的n次多项式为：其中基函数满足上述多项式插值又称为n次Lagrange插值。现在学习的是第9页，共69页说明：1、多项式插值的基函数仅与节点有关，而与被插值的原函数 无关；2、插值多项式仅由数对 确定，而与数对的排列次序无关。3、多项式插值除拉格朗日多项式插值法外，还有 牛顿(Ne

5、wton)插值法、埃尔米特(Hermite)插值法、三次样条插值法等，可参看有关数值分析的书籍。其中Newton插值是拉格朗日插值的一种等价变形，Hermite插值一种带导数插值条件的插值。现在学习的是第10页，共69页例例 将 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

6、)=8(x-/4)(x-/2)/2-16x(x-/2)0.7071/2 cos(/6)=P2(/6)0.8508 精确值：精确值：cos(/6)0.8660现在学习的是第11页，共69页下面来求解引例1(课堂练习)。引例1：对于下面给定的4组数据，求在x=175处 y 的值。x144169225y121315解：用一次拉格朗日插值：所以取 为插值节点，则 计算得 因为插值点 位于 和 之间，于是 现在学习的是第12页，共69页用二次拉格朗日插值：取，则 计算得，于是（的准确值为）由上例看出，二次插值的精度明显要比一次插值要高。但对于拉格朗日多项式插值，是否插值其精度就一定越高呢？现在学习的是第

7、13页，共69页答案是：对于某些函数，适当地提高插值多项式的 次数，会提高计算精度。但与此同时，多项式的次数 增大可能造成插值函数的收敛性和稳定性越来越差，逼近的效果往往不理想，一个典型的例子是函数 选取不同插值节点个数 n+1，其中 n 为插值多项式的 次数，使得它在结点的值与被插函数在对应结点的值相等。当n分别取2,4,6,8,10时，绘出的插值图形如下。现在学习的是第14页，共69页从图中可看出，图形显示出振荡现象，在5和-5附近误差很大。这种现象叫做Runge 现象。现在学习的是第15页，共69页这说明，在大范围内使用高次插值，逼近效果往往并不理想。解决此问题的思路是化整为零，采用分段

8、 插值，即在小范围内使用低次多项式插值。不是去寻求整个插值区间上的一个高次多项式，也就是说插值区间划分为若干个小区间，在每个小区间上用低而是把次多项式插值，在整个插值区间上就得到一个分段插值函数。区间的划分可以是任意的，各个区间上插值多项 式的次数的选取也可按具体问题选择。在分段插值中，较为简单的是分段线性插值。现在学习的是第16页，共69页实际数学建模中，在光滑性要求不高的条件下，分段线性或二次插值基本可以满足需要。问题中提出的插值问题，有一些插值函数曲线要求然而实际具有较高的光滑性，如飞机机翼的下轮廓线。分段线性插值虽然简单，但插值函数在结点处的 一阶导数一般不存在，光滑性不高，样条插值的

9、提出。这就导致了三次在数学上，光滑程度的定量描述是：函数(曲线)的 k 阶导数存在且连续，则称该曲线具有 k 阶光滑性。现在学习的是第17页，共69页光滑性的阶次越高，则越光滑。分段多项式达到较高阶光滑性的方法？是否存在较低次的就是一个很好的例子。三次样条插值现在学习的是第18页，共69页3 三次样条插值三次样条插值三次样条插值是一种非常有效的插值方法，它在实际工程中有着非常重要的应用。三次样条插值的理论推导是比较复杂的，但在数学软件MATLAB中有现成的调用程序，这样我们就可直接借助计算机来进行运算。下面简单介绍一下三次样条插值的基本原理。定义：设给定区间 上的一个划分 如果函数 满足条件：

10、现在学习的是第19页，共69页定义：设给定区间 上的一个划分 如果函数 满足条件：（1）在每个子区间 是三次多项式；（2）在区间 上连续，记作（3）对于在节点上给定的函数值 满足 则称 为 在区间 上的三次样条插值函数三次样条插值函数。现在学习的是第20页，共69页简单地说，已经知道函数 在节点 上的函数值 多项式函数，现要求一个三次，使满足且。由定义可知，是区间 上的分段分段三次插值多项式，即现在学习的是第21页，共69页由于，这个函数的曲线具有二阶光滑度，看起来就很光顺了，能满足一般工程上的需要。其中 是子区间 插值于两点 的三次多项式，即 下面简单介绍一下三次样条插值函数的推导。现在学习

11、的是第22页，共69页现要求 为待定系数，共4n个。已知条件：1）共 n+1个方程；2）共 3(n-1)个方程。现在学习的是第23页，共69页现要求 4n 个待定系数，但只有(n+1)+3(n-1)=4n-2个方程，故需要补充两个方程，即所谓的边界条件边界条件。通常有以下三类边界条件：3.1）给定两个端点 处的导数，即3.2）给定两个端点 处的导数 即3.3）周期性条件，即现在学习的是第24页，共69页4 用用MATLAB软件求解插值问题软件求解插值问题在MATLAB中提供了一个一维插值函数interp1，它的调用格式为cy=interp1(x,y,cx,method)其中x、y是所给数据的横

12、纵坐标，要求x的分量按升序或降序排列，cx是待求的插值点的横坐标，返回值cy是待求的插值点的纵坐标，method是插值方法，该函数提供了四种可选的插值方法：现在学习的是第25页，共69页nearest最邻近点插值。点和这两已知点间位置的远近来进行插值，取较近已知它根据已知两点间的插值插值点处的函数值作为未知插值点处的函数值。linear线性插值。它将相邻的数据点用直线相连，按所生成的直线进行插值。spline三次样条插值。它利用已知数据求出样条 函数后，按样条函数进行插值。cubic三次插值。它利用已知数据求出三次多项式函数后，按三次多项式函数进行插值。缺省时插值方法为分段线性插值。现在学习的

13、是第26页，共69页下面用该函数来求解下列插值问题。对于下面给定的4组数据，求在x=110处 y的值。x100121144169y10111213输入命令：x=100 121 144 169;y=10 11 12 13;cx=110;cy=interp1(x,y,cx,linear);现在学习的是第27页，共69页运行结果为cy=10.4762。由于线性插值只需要两个点，因而在上述命令中实际上只用了前两个点。若将最后一个命令中的method改为缺省、nearest、cubic和spline，运行结果为依次为 cy=10.4762、cy=10、cy=10.4869、cy=10.4877 通过比较

14、，显然三次样条插值的结果最好。现在学习的是第28页，共69页例：在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。程序：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,+,h,t,r:)%作图xlabel(Hour),ylabel(Degrees Celsius)现在学习的是第29页，

15、共69页现在学习的是第30页，共69页xy机翼下轮廓线例 已知飞机下轮廓线上数据如下，求X每改变0.1时的Y值。现在学习的是第31页，共69页程序：程序：lch(lagr1)lch(lagr1)function y=lagr1(x0,y0,x)function y=lagr1(x0,y0,x)n=length(x0);m=length(x);n=length(x0);m=length(x);for i=1:mfor i=1:m z=x(i);z=x(i);s=0.0;s=0.0;for k=1:n for k=1:n p=1.0;p=1.0;for j=1:n for j=1:n if j=k

16、 if j=k p=p*(z-x0(j)/(x0(k)-x0(j);p=p*(z-x0(j)/(x0(k)-x0(j);end end end end s=p*y0(k)+s;s=p*y0(k)+s;end end y(i)=s;y(i)=s;endend现在学习的是第32页，共69页x0=0 3 5 7 9 11 12 13 14 15;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.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,

17、1)plot(x0,y0,k+,x,y1,r)gridtitle(lagrange)subplot(3,1,2)plot(x0,y0,k+,x,y2,r)gridtitle(piecewise linear)subplot(3,1,3)plot(x0,y0,k+,x,y3,r)gridtitle(spline)现在学习的是第33页，共69页现在学习的是第34页，共69页曲线拟合是指：已知平面上 n 个点（xi,yi)i=1,n,寻求一个函数（曲线）y=f(x),使 f(x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。+xyy=f(x)(xi,yi)ii 为点（xi,yi)与曲线 y=f

18、(x)的距离5 曲线拟合的基本原理曲线拟合的基本原理现在学习的是第35页，共69页拟合与插值的区别拟合与插值的区别函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。问题：问题：给定一批数据点，需确定满足特定要求的曲线或曲面。解决方案：解决方案：若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是曲线数据拟合，又称曲线拟合或曲面拟合。若要求所求曲线（面）通过所给所有数据点，就是插值问题；现在学习的是第36页，共69页根据曲线拟合问题的定义，其关键在于准则的选取，选取的准则不同，其对应的拟合方法及其复杂程度也不相同。对

19、于一维曲线拟合，设n个不同的离散数据点为，要寻找的拟合曲线方程为 记拟合函数在 处的偏差为 常用的准则有：准则1：选取，使所有偏差的绝对值之和最小，即 现在学习的是第37页，共69页准则2：选取，使所有偏差的绝对值的最大值最小，即准则3：选取，使所有偏差的平方和最小，即 相对而言，准则3最便于计算，因而通常根据准则3 来选取拟合曲线。准则3又称为最小二乘准则，对应的曲线拟合方法称为最小二乘法。现在学习的是第38页，共69页线性最小二乘法的基本思路线性最小二乘法的基本思路第一步:先选定一组函数 r1(x),r2(x),rm(x),mn,令 f(x)=a1r1(x)+a2r2(x)+amrm(x)

20、其中 a1,a2,am 为待定系数。第二步:确定a1,a2,am 的最小二乘准则：使n个点（xi,yi)与曲线 y=P(x)的距离i 的平方和最小。记问题归结为，求 a1,a2,am 使 J(a1,a2,am)最小。现在学习的是第39页，共69页线性最小二乘法的求解线性最小二乘法的求解要使 J(a1,a2,am)最小，的必要条件得则由多元函数取得极值即亦即是未知量 的线性方程组，称之为正定方程组。现在学习的是第40页，共69页选定一组函数 组解出 后，就可由正规方程，于是就可得线性最小二乘拟合函数所给数据的散点图，观察数据所呈现出来的曲线的大致一般的做法是首先绘出形状，再结合该问题所在专业领域

21、内的相关规律和结论，来确定拟合函数的形式。实际操作时可在直观判断的基础上，选几种常用的曲线分别进行拟合，比较选择拟合效果最好的曲线。面对一组数据，作线性最小二乘拟合时，恰当选定函是一个难点。数现在学习的是第41页，共69页常用的曲线有直线、多项式、双曲线和指数曲线等。另外，曲线拟合又可分为线性曲线拟合和非线性曲线 拟合。一般地，如果拟合函数中的系数 以线性形式出现，全部如拟合函数 为线性拟合，也称为多项式拟合；若拟合函数中的系数 不能全部以线性形式出现，如指数拟合函数为非线性曲线拟合。实际应用中，多项式最小二乘拟合用的较多，MATLAB中也有专用函数。现在学习的是第42页，共69页线性最小二乘

22、拟合线性最小二乘拟合 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)：现在学习的是第43页，共69页6 用用MATLAB软件求解拟合问题软件求解拟合问题在MATLAB中提供了一个多项式最小二乘拟合函数 polyfit(x,y,n)，它的调用格式为P=polyfit(

23、x,y,n)拟合多项式按自变量拟合多项式按自变量降幂排列的系数向量降幂排列的系数向量 输入同长度输入同长度的数组的数组X，Y拟合多项拟合多项式次数式次数下面用该函数来求解拟合问题引例2：现在学习的是第44页，共69页例2：观测物体的直线运动，得以下数据，求物体的运动方程。t（秒）00.91.93.03.95.0s（米）010305080110输入命令：t=0 0.9 1.9 3 3.9 5;s=0 10 30 50 80 110;plot(t,s,*-)xlabel(运动时间 t（秒）)ylabel(运动位移 s（米）)gtext(物体运动的时间与位移散点图)现在学习的是第45页，共69页下面

24、显示的是物体运动的时间与位移散点图：现在学习的是第46页，共69页不难看出图形近似为一条直线，因此猜测用一次多项式来拟合，输入命令：P=polyfit(t,s,1)运行结果为：P=22.2538 -7.8550即 下面绘出的是拟合曲线和散点图对比图形，现在学习的是第47页，共69页可以看出拟合效果并不理想。现在学习的是第48页，共69页根据物理学中物体运动的方程，我们用二次曲线来拟合，输入命令：P=polyfit(t,s,2)得到拟合函数为：对比图形如下，现在学习的是第49页，共69页现在学习的是第50页，共69页可见曲线拟合本身就是一个猜测的过程，通常是不断地修正拟合函数，使拟合效果达到满意

25、的程度。可以看出拟合效果有明显改善，拟合曲线与散点图 基本上是吻合的，因此该物体运动的方程是现在学习的是第51页，共69页7 建模案例建模案例(1992年年A题：农作物施肥效果分析题：农作物施肥效果分析)某地区作物生长所需要的营养元素主要有氮（N）、钾（K）、磷（P）。某作物研究所在该地区对土豆 与生菜做了一定数量的实验，实验数据如下列表格所示，其中ha表示公顷，t表示吨，kg表示公斤。当一个营养素的施肥量变化时，总将另两个营养素的施肥量保持在第七个水平上，于N的施肥量做实验时，P与K的施肥量分别取为 如对土豆产量关196kg/ha与372kg/ha。试分析施肥量与产量之间的关系。现在学习的是

26、第52页，共69页土豆：土豆：N P K施肥量(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

27、现在学习的是第53页，共69页生菜生菜：N P K施肥量(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.111

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

29、的，互不影响。（此假设澄清了在连续进行的实验中，后期实验产量与前期施肥无关）现在学习的是第56页，共69页符号说明：符号说明：农作物产量；：施肥量；N、K、P：氮、磷、钾肥的施肥量；：农产品价格；：肥料价格；Tn、Tp、Tk：氮、磷、钾肥的价格；：常数。现在学习的是第57页，共69页问题分析：问题分析：1、普遍规律施肥量与产量满足下图所示关系，它分为三个不同的区段，第二区段，随着施肥量的增加，作物产量平缓上升，一定限度后，第三区段，当施肥量超过产量反而随施肥量的增加而减少。施肥量的增加而迅速增加，在第一区段，当施肥量较小时，作物产量随现在学习的是第58页，共69页现在学习的是第59页，共69页

30、2、数据分析通过绘制散点图，初步得到农作物产量与施肥量间的定性认识。现在学习的是第60页，共69页从散点图可以发现，氮肥施加量与农作物的产量大致呈指数关系，磷肥施加量与农作物产量大致呈分段直线关系，钾肥施加量与土豆产量大致呈指数关系，与生菜产量产量关系规律不明显。但有一点，钾肥施加量的增加时，生菜产量上升幅度不大，波动也不大，这说明钾肥对生菜产量的影响较小。3、理论支撑1）Nicklas和Miller理论：设h为达到最高产量时的施 肥量，边际产量（即产量W对施肥量 x的导数）与 成正比例关系，即现在学习的是第61页，共69页这是一个简单的微分方程，积分得即 2）米采利希学说：只增加某种养分时，

31、引起产量的增加与该种养分供应充足时达到的最高产量 A 和现在产量 W 之差成正比，即现在学习的是第62页，共69页这是一个可分离变量的微分方程，积分得即 3）博伊德理论：英国科学家博伊德发现，在某些 情况下，将施肥对象按施肥水平分成几组，则各组效应曲线就呈直线形式，若按水平分为两组，可以用下式表示：现在学习的是第63页，共69页模型建立：模型建立：本问题的性质是确定变量间的相关关系，理论和数据分析的结果，可以初步确定施肥量与产量结合前述之间的关系，即拟合曲线的形式。具体地，氮肥对土豆和生菜的效应方程为：磷肥对土豆的效应方程为：现在学习的是第64页，共69页磷肥对生菜的效应方程为：钾肥对土豆的效

32、应方程为：钾肥对生菜的效应方程为：现在学习的是第65页，共69页模型求解：模型求解：上述模型中的拟合曲线既有线性的，又有非线性的，下面仅对氮肥对土豆的效应方程进行求解。txn=0,34,67,101,135,202,259,336,404,471;tyn=15.18,21.36,25.72,32.29,34.03,39.45,43.15,43.46,40.83,30.75;a,s=polyfit(txn,tyn,2);输入命令：运行结果：a=-0.0003 0.1971 14.7416因此氮肥对土豆的效应的拟合曲线方程是：现在学习的是第66页，共69页散点图与拟合曲线对比图形如下：现在学习的是第67页，共69页上图显示出拟合效果较好。对于二维插值方法的基本原理与MATLAB求解、非线性拟合的MATLAB求解，由于时间关系没有讲解以后有机会再和大家交流。大家可以有兴趣自己上机求解。对于其余的几种情况，现在学习的是第68页，共69页练习题：练习题：(2)下面数据是某次实验所得，希望得到 x 和 f 的关系 x1247912131517f1.53.96.611.715.618.819.620.621.1(1)已知下列表值x0123y13927试用线性与二次Lagrange插值多项式分别计算当x=1.25时y的值。现在学习的是第69页，共69页