插值与数据拟合优秀PPT.ppt

《插值与数据拟合优秀PPT.ppt》由会员分享，可在线阅读，更多相关《插值与数据拟合优秀PPT.ppt（78页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、插值与数据拟合你现在浏览的是第一页，共78页 插值与数据拟合就是通过一些已知数据插值与数据拟合就是通过一些已知数据去确定某类函数的参数或寻找某个近似函数，使所去确定某类函数的参数或寻找某个近似函数，使所得的函数与已知数据具有较高的精度，并且能够使得的函数与已知数据具有较高的精度，并且能够使用数学分析的工具分析数据所反映的对象的性质用数学分析的工具分析数据所反映的对象的性质 几种常用的方法：几种常用的方法：1 1、一般插值法、一般插值法 2 2、样条插值法、样条插值法 3 3、最小二乘曲线、最小二乘曲线 4 4、曲面的拟合、曲面的拟合数据拟合与插值建模数据拟合与插值建模你现在浏览的是第二页，共7

2、8页 上大学二年级的小华正在做作业，上大学二年级的小华正在做作业，“爸爸，计算这爸爸，计算这道题要用到道题要用到sin sin ，可是我的计算器坏了，怎么办。，可是我的计算器坏了，怎么办。”当工程师的老张从厚厚的一摞旧书底下抽出一本数当工程师的老张从厚厚的一摞旧书底下抽出一本数学用表来，学用表来，“给你，这是我念大学时用的，那时候啊，给你，这是我念大学时用的，那时候啊，计算器听都没听说过。计算器听都没听说过。”小华拿着表翻了一会儿，无奈的说：小华拿着表翻了一会儿，无奈的说：“表上每表上每 才才有一个函数值，这里只有一个函数值，这里只sin sin 和和sin sin “表中没表中没有的，都可以

3、用插值方法计算有的，都可以用插值方法计算”“插值！我们的数学实插值！我们的数学实验课就要学了，不过，今天我要先自己想个办法，用验课就要学了，不过，今天我要先自己想个办法，用这个算出这个算出sin sin ”这本四位数学用表给出这本四位数学用表给出sin sin 0.5760.576，sin =0.5783sin =0.5783。小华认为在。小华认为在sin sin 到到sin sin 这样小的范围内，正弦可以近似为线性函数，于是很容这样小的范围内，正弦可以近似为线性函数，于是很容易地得到易地得到Sin =0.576+(0.5783-0.5760)0.6=0.5774Sin =0.576+(0.

4、5783-0.5760)0.6=0.5774你现在浏览的是第三页，共78页 聪明的小华用的这个办法是一种插值聪明的小华用的这个办法是一种插值方法方法分段线性插值。实际上，分段线性插值。实际上，插值可插值可以理解为，要根据一个用表格表示的函数，以理解为，要根据一个用表格表示的函数，计算表中没有的函数值计算表中没有的函数值，表中有的。表中有，表中有的。表中有的，如（的，如（sin sin ，0.57600.5760）(sin (sin ,0.5783),0.5783)称为称为节点节点；要计算的，如；要计算的，如sin sin ，称为，称为插值点插值点，结果（，结果（0.57740.5774）即为）

5、即为插值插值。小华作的线性函数为小华作的线性函数为插值函数插值函数，插值函数所，插值函数所表示的直线当然要通过节点。表示的直线当然要通过节点。你现在浏览的是第四页，共78页 插值最初来源于天体计算插值最初来源于天体计算由若干观测由若干观测值（即节点）计算任意时刻星球的位置（即插值值（即节点）计算任意时刻星球的位置（即插值点和插值）点和插值）的需要。现在，虽然人们已很的需要。现在，虽然人们已很少需要用它从函数表计算函数值了，但是插少需要用它从函数表计算函数值了，但是插值仍然在诸如机械加工等工程技术和数据处值仍然在诸如机械加工等工程技术和数据处理等科学研究中有着许多直接的应用，另一理等科学研究中有

6、着许多直接的应用，另一方面，插值又是数值微分、数值积分、常微方面，插值又是数值微分、数值积分、常微分方程数值等数值计算的基础。分方程数值等数值计算的基础。你现在浏览的是第五页，共78页 几天后，小华在物理实验里又碰到一几天后，小华在物理实验里又碰到一个看起来非常类似的问题：有一只对温度个看起来非常类似的问题：有一只对温度敏感的电阻，已经测得了一组温度敏感的电阻，已经测得了一组温度T T和电和电阻阻R R数据。数据。现在想知道现在想知道 时的电阻多大。时的电阻多大。温度温度t(0C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032你现在浏览的

7、是第六页，共78页 小华征求老师的意见，老师给了他两点提示：小华征求老师的意见，老师给了他两点提示：一是在直角坐标系中把一是在直角坐标系中把5 5个点（个点（T T，R R）画一下，看）画一下，看看电阻看电阻R R和温度和温度T T之间大致有什么关系之间大致有什么关系；二是；二是测量数测量数据总有相当大的误差据总有相当大的误差，这与用函数表作插值计算应，这与用函数表作插值计算应该有不同之处吧（虽然函数表也存在舍入误差，但该有不同之处吧（虽然函数表也存在舍入误差，但很小，可以认为表中数值是精确的）很小，可以认为表中数值是精确的）通过图形小华看到，通过图形小华看到，R R与与T T大致呈直线关系，

8、大致呈直线关系，于是用手画了一条靠近于是用手画了一条靠近5 5个点的直线，又想起中学个点的直线，又想起中学物理学过，金属材料的电阻率与温度成正比，从物理学过，金属材料的电阻率与温度成正比，从而确定而确定R R与与T T的关系应该是的关系应该是 R=at+b R=at+b 其中其中a a，b b为待定常数。为待定常数。你现在浏览的是第七页，共78页 正是由于测量误差的存在，由正是由于测量误差的存在，由R=at+bR=at+b表表示的直线不可能通过全部示的直线不可能通过全部5 5个点，所以，与插个点，所以，与插值曲线要通过全部节点不同，小华打算作一条值曲线要通过全部节点不同，小华打算作一条尽量靠近

9、所有的点的直线，求出尽量靠近所有的点的直线，求出a a，b b待定常数，待定常数，由此计算由此计算t=t=的的R R就十分简单了。就十分简单了。你现在浏览的是第八页，共78页 根据一组（二组）数据，即平面上的根据一组（二组）数据，即平面上的若干点，确定一个一元函数，即曲线，若干点，确定一个一元函数，即曲线，使这些节点与曲线总体来说尽量接近，使这些节点与曲线总体来说尽量接近，这就是曲线拟合。这就是曲线拟合。函数值与曲线拟合都是要根据一组数函数值与曲线拟合都是要根据一组数据构造一个函据构造一个函 数作为近似，由于近似的要数作为近似，由于近似的要求不同，二者的数学方法是完全不同的。求不同，二者的数学

10、方法是完全不同的。你现在浏览的是第九页，共78页 数数 据据 拟拟 合合1.拟合的基本原理；拟合的基本原理；2.最小二乘拟合；最小二乘拟合；3.用用Matlab作最小二乘拟合；作最小二乘拟合；4.4.如何用拟合解决实际问题。如何用拟合解决实际问题。你现在浏览的是第十页，共78页t(h)0.25 0.5 1 1.5 2 3 4 6 8c(g/ml)19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01 对某人用快速静脉注射方式一次性注射某种药物对某人用快速静脉注射方式一次性注射某种药物300mg后，经过时间后，经过时间t采集血样，测得血药浓度采集血样，

11、测得血药浓度c如下表：如下表：求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).半对数坐标系半对数坐标系(semilogy)下的图形下的图形Log10c(t)=a t+b引例引例1 1：血药浓度的变化规律：血药浓度的变化规律你现在浏览的是第十一页，共78页曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法 已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi)i=1,n,寻求一个函数（曲线）寻求一个函数（曲线）y=f(x),使使 f(x)在某种准在某种准则下与所有数据点最为接近，即曲线拟合得最好。则下与所有数据点最为接近，即曲线拟合得最好。+xyy

12、=f(x)(xi,yi)i i 为点为点（xi,yi)与与曲线曲线 y=f(x)的距离的距离你现在浏览的是第十二页，共78页最小二乘拟合最小二乘拟合 第一步：第一步：先选定一类函数先选定一类函数f(x，a1，a2，am)其其准则为（最小二乘准则）：使准则为（最小二乘准则）：使n个点个点（xi,yi)与与曲线曲线 y=f(x，a1，a2，am)的距离的距离 i 的平方和最小的平方和最小。其中其中a1,a2,am为待定常数。为待定常数。f f可以为一些简单的可以为一些简单的“基函数基函数”（如幂函数，三角函数等等）的线性组合：（如幂函数，三角函数等等）的线性组合：第二步：确定参数第二步：确定参数a

13、1,a2,am，你现在浏览的是第十三页，共78页记记问题归结为，求问题归结为，求 a1,a2,am 使使 J(a1,a2,am)最小。最小。这样的拟合称为最小二乘拟合。这样的拟合称为最小二乘拟合。除了最小二乘准则（即各点误差的平方和除了最小二乘准则（即各点误差的平方和最小），你认为还可以用怎样的拟合准则？最小），你认为还可以用怎样的拟合准则？比较起来，最小二乘准则有什么优点？比较起来，最小二乘准则有什么优点？思考思考你现在浏览的是第十四页，共78页最小二乘拟合函数最小二乘拟合函数 f(x f(x，a a1 1，a am m)的选取的选取 +f=a1+a2xf=a1+a2x+a3x2f=a1+a

14、2x+a3x2f=a1exp(a2x)+f=a1exp(a2x)1.1.通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 f f；2.2.将数据将数据 (xi,yi)i=1,n 作图，通过直观判断确定作图，通过直观判断确定 f：你现在浏览的是第十五页，共78页2.2.作一般的作一般的最小二乘曲线拟合，可利用已有程序最小二乘曲线拟合，可利用已有程序curvefit,curvefit,其调用格式为：其调用格式为：a=curvefit(f,a0,x,y)1.1.作多项式作多项式f(x)=a1xm+amx+am+1函数函数拟合拟合,可利可利用已有程序用已有程序polyfit,polyfit

15、,其调用格式为其调用格式为:a=polyfit(x,y,m)用用MATLABMATLAB作最小二乘拟合作最小二乘拟合数据点数据点拟合多项式次数拟合多项式次数系数系数注：注：f f为拟合函数为拟合函数y=f(a,x)y=f(a,x)的函数的函数M M文件，文件，f(a,x)f(a,x)为拟合函数。为拟合函数。数据点数据点待定常数待定常数a a的初值的初值函数函数M文件文件你现在浏览的是第十六页，共78页用用MATLABMATLAB作多项式最小二乘拟合作多项式最小二乘拟合2.2.用命令用命令polyfit(x,y,m)得到得到 a1=3.3940,a2=702.49181.选取函数选取函数 R=a

16、1t+a2温度温度t(0C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032例例.由数据由数据拟合拟合R=f(t)你现在浏览的是第十七页，共78页用用MATLABMATLAB作最小二乘曲线拟合作最小二乘曲线拟合例：用函数例：用函数f(x)=a1*exp(-a2*x)+a3*exp(-a4*x)拟合拟合下列数据点：下列数据点：xdata=0:.1:2 ydata=5.8955 3.5639 2.5173 1.9790 1.8990 1.3938 1.1359 1.0096 1.0343 0.8435 0.6856 0.6100 0.5392

17、0.3946 0.3903 0.5474 0.3459 0.1370 0.2211 0.1704 0.2636用命令用命令curvefit(f,a0,x,y)你现在浏览的是第十八页，共78页拟合的应用拟合的应用参数辨识参数辨识 数学建模的方法：机理分析和测试分析。数学建模的方法：机理分析和测试分析。机理分析是根据对客观事物特性的认识，找出机理分析是根据对客观事物特性的认识，找出反映内部机理的数量规律，建立的模型常有明确的反映内部机理的数量规律，建立的模型常有明确的物理意义。物理意义。测试分析将研究的对象看作一个测试分析将研究的对象看作一个“黑箱黑箱”，通，通过对实验数据的统计分析，找出与数据拟

18、合得最好过对实验数据的统计分析，找出与数据拟合得最好的模型。的模型。机理分析机理分析模型结构模型结构实验数据实验数据未知参数未知参数你现在浏览的是第十九页，共78页范例：薄膜渗透率的测定范例：薄膜渗透率的测定 一、问题：一、问题：某种医用薄膜，具有从高浓度的溶液向低浓某种医用薄膜，具有从高浓度的溶液向低浓度的溶液扩散的功能，在试制时需测定薄膜被物度的溶液扩散的功能，在试制时需测定薄膜被物质分子穿透的能力。质分子穿透的能力。测定方法：用面积为测定方法：用面积为S S的薄膜将容器分成体的薄膜将容器分成体积分别为积分别为 的两部份，在两部分中分别注的两部份，在两部分中分别注满该物质的两种不同浓度的溶

19、液。此时该物质分满该物质的两种不同浓度的溶液。此时该物质分子就会从高浓度溶液穿过薄膜向低浓度溶液中扩子就会从高浓度溶液穿过薄膜向低浓度溶液中扩散。平均每单位时间通过单位面积薄膜的物质分散。平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，比例系数子量与膜两侧溶液的浓度差成正比，比例系数K K表征了薄膜被该物质分子穿透的能力，称为渗透表征了薄膜被该物质分子穿透的能力，称为渗透率。定时测量容器中薄膜某一侧的溶液浓度，以率。定时测量容器中薄膜某一侧的溶液浓度，以此确定此确定K K。VAVBS你现在浏览的是第二十页，共78页二、问题分析二、问题分析 考考察察时时段段tt，t+tt+

20、t薄薄膜膜两两侧侧容容器器中中该该物物质质质质量量的变化。的变化。设设 ，对对容容器器的的B B部部分分溶溶液液浓度的测试结果如下表：（浓度单位浓度的测试结果如下表：（浓度单位 ）1）在容器的一侧，物质质量的增加是由于另一在容器的一侧，物质质量的增加是由于另一侧的物质向该侧渗透的结果，因此物质质量的增侧的物质向该侧渗透的结果，因此物质质量的增量应等于另一侧的该物质向这侧的渗透量。量应等于另一侧的该物质向这侧的渗透量。你现在浏览的是第二十一页，共78页 以容器以容器A A侧为例，在时段侧为例，在时段tt，t+tt+t物质质量物质质量的增量为：的增量为：分别表示在时刻分别表示在时刻t t膜两侧溶液

21、膜两侧溶液设设的浓度，浓度单位的浓度，浓度单位:由于平均每单位时间通过单位面积薄膜的物由于平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，比例系质分子量与膜两侧溶液的浓度差成正比，比例系数为数为K K。因此，在时段因此，在时段tt，t+tt+t，从，从B B侧渗透至侧渗透至A A侧侧的该物质的质量为：的该物质的质量为：你现在浏览的是第二十二页，共78页于是有：于是有：两边除以两边除以tt，并令，并令t0t0取极限再稍加整理即得：取极限再稍加整理即得：分别表示在初始时刻两侧溶液的浓度分别表示在初始时刻两侧溶液的浓度其中其中（1）2)注意到整个容器的溶液中含有该物质的质量不变

22、注意到整个容器的溶液中含有该物质的质量不变,与与初始时刻该物质的含量相同，因此初始时刻该物质的含量相同，因此 你现在浏览的是第二十三页，共78页从而：从而：加上初值条件：加上初值条件：代入式（代入式（1）得：）得：便可得出便可得出CB(t)的变化规律，从而根据实验数据进行的变化规律，从而根据实验数据进行拟合，估计出参数拟合，估计出参数K,。你现在浏览的是第二十四页，共78页三、数学模型三、数学模型假设：假设：1）薄膜两侧的溶液始终是均匀的；薄膜两侧的溶液始终是均匀的；2）平均每单位时间通过单位面积薄膜的物质分子量平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比。与膜两侧溶液

23、的浓度差成正比。3 3）薄膜是双向同性的即物质从膜的任何一侧向另一侧）薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的性能是相同的。渗透的性能是相同的。基于假设和前面的分析，基于假设和前面的分析，B B侧的浓度侧的浓度CB(t)应满足应满足如下微分方程和初始条件：如下微分方程和初始条件：你现在浏览的是第二十五页，共78页四、求解方法：四、求解方法：1.函数拟合法函数拟合法前面得到的模型是一个带初值的一阶线性微分方前面得到的模型是一个带初值的一阶线性微分方程，解之得：程，解之得：问题归结为利用问题归结为利用C CB B在时刻在时刻t tj j的测量数据的测量数据C Cj j(j=1,2,.,N

24、)(j=1,2,.,N)来辨识来辨识 K K 和和 。你现在浏览的是第二十六页，共78页引入引入从而从而 用函数用函数CB(t)来拟合所给的实验数据，从来拟合所给的实验数据，从而估计出其中的参数而估计出其中的参数a,b,K。将将代入上式有：代入上式有：你现在浏览的是第二十七页，共78页用用MATLABMATLAB软件进行计算软件进行计算.1 1）编写函数编写函数M-M-文件文件 nongdu.m nongdu.mfunction f=nongdu(x,tdata)function f=nongdu(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata);f=x(1)

25、+x(2)*exp(-0.02*x(3)*tdata);其中其中 x(1)=a;x(2)=b;x(3)=k;x(1)=a;x(2)=b;x(3)=k;2)2)在工作空间中执行以下命令在工作空间中执行以下命令(test1.m)(test1.m)tdata=linspace(100,1000,10);tdata=linspace(100,1000,10);cdata=4.54 4.99 5.35 5.65 5.90 6.10.cdata=4.54 4.99 5.35 5.65 5.90 6.10.6.26 6.39 6.50 6.59;6.26 6.39 6.50 6.59;x0=0.2,0.05

26、,0.05;x0=0.2,0.05,0.05;x=curvefit(x=curvefit(nongdunongdu,x0,tdata,cdata),x0,tdata,cdata)3)3)输出结果输出结果:x=0.007 -0.003 0.1012:x=0.007 -0.003 0.1012 即即 k=0.1012,a=0.007,b=-0.003,k=0.1012,a=0.007,b=-0.003,你现在浏览的是第二十八页，共78页进一步求得：进一步求得：2.非线性规划法非线性规划法 利用利用C CB B在时刻在时刻t tj j的测量数据的测量数据C Cj j(j=1,2,.,N)(j=1,2

27、,.,N)来辨识来辨识 K K 和和 。问题可转化为求函数问题可转化为求函数即求函数即求函数的最小值点（的最小值点（K K，a a，b b）。）。你现在浏览的是第二十九页，共78页3.导函数拟合法导函数拟合法前面得到的微分方程为：前面得到的微分方程为：令令上式变为：上式变为：这可以看作这可以看作随随CB的变化规律的变化规律(j=1,2,.,N)若知道一组数据若知道一组数据则可用最小二乘拟合的方法来求出函数则可用最小二乘拟合的方法来求出函数中的未知参数中的未知参数K和和h。你现在浏览的是第三十页，共78页即为求参数即为求参数K,a使下列误差函数达到最小：使下列误差函数达到最小：该问题等价于用函该

28、问题等价于用函 数数 f(K,a,CB)=K(0.01a-0.02CB)来拟合数据来拟合数据(j=1,2,.,N)你现在浏览的是第三十一页，共78页用用MATLABMATLAB软件进行计算软件进行计算.求数据点求数据点(j=1,2,.,N)tdata=linspace(100,1000,10);cdata=1e-05.*454 499 535 565 590.610 626 639 650 659;d,ifail=e01bef(tdata,cdata);cj,dcj=e01bgf(tdata,cdata,d,tdata);1 1）编写函数编写函数M-M-文件文件 baomof.m baomof

29、.mfunction f=baomof(x,cdata)f=x(1)*(0.01*x(2)-0.02*cdata)其中其中 x(1)=K;x(2)=h2)2)编写命令编写命令M M文件文件(baomo21.m)(baomo21.m)你现在浏览的是第三十二页，共78页3)3)输出结果输出结果:x=0.1009 0.014:x=0.1009 0.014 即即 k=0.1009,h=0.014 k=0.1009,h=0.014作函数拟合作函数拟合x0=0.2,0.1;x=curvefit(baomof,x0,cdata,dcj)你现在浏览的是第三十三页，共78页4.线性化迭代法线性化迭代法前面带初始

30、条件的一阶线性微分方程的解为前面带初始条件的一阶线性微分方程的解为其中：其中：如果得到了参数如果得到了参数K的一个较好的近似值的一个较好的近似值K*，则，则将将CB(t)关于关于K在在K*处展开，略去处展开，略去K的二次及以上的的二次及以上的项得项得CB(t)的一个近似式的一个近似式通过极小化你现在浏览的是第三十四页，共78页确定a,b,d,再由K=d/0.02b得到K*的修正值K。K*K*-K,得到K的一个新的近似值，用同样的方法再求新的修正值K。这个过程可以不断重复，直到修正值足够小为止。1）当K的初值取为k=0.3时，出现奇异情况，迭代不收敛；2）当K的初值取为k=0.2时，经四次迭代，

31、已经收敛到一个很好的解。迭代结果如下表。你现在浏览的是第三十五页，共78页五、结果及误差分析五、结果及误差分析 几种方法得出的结果及相应的误差总结于几种方法得出的结果及相应的误差总结于下表，误差为计算数据与实验数据之差的平方下表，误差为计算数据与实验数据之差的平方和。和。注：导函数拟合法得出的参数值精度有限，线性化迭代法要求参数的初值比较接近精确值。因此可将导函数拟合法和线性化迭代法结合起来使用，把前者得到的参数K的值作为迭代法中K的初值，这样可使迭代法收敛或收敛更快。3）取K的初值为k=0.1009,只一次迭代就得到2）中的最后结果。你现在浏览的是第三十六页，共78页你现在浏览的是第三十七页

32、，共78页函数拟合法的拟合效果函数拟合法的拟合效果你现在浏览的是第三十八页，共78页求解参数辨识模型的方法：求解参数辨识模型的方法：函数拟合；函数拟合；非线性规划；非线性规划；导函数拟合；导函数拟合；线性化迭代；线性化迭代；其它方法其它方法。你现在浏览的是第三十九页，共78页用用Logistic模拟水稻叶伸长生长模拟水稻叶伸长生长时间11.82.63.44.14.85.46.16.87.48.1重量0.30.50.91.42.53.24.37.610.114.418.5时间8.89.410.110.811.712.413.114.415.115.7重量23.025.230.433.738.84

33、1.743.744.845.545.3生长观测记录数据你现在浏览的是第四十页，共78页模型表达式：模型表达式：你现在浏览的是第四十一页，共78页程序！程序！你现在浏览的是第四十二页，共78页关于关于polyfit命令命令命令：p=polyfit(x,y,n)（1）x与y为模拟数据（2）n为拟合多项式的次数（3）当n=1时为用最小二乘法进行直线拟合（4）得到的向量p为长度n+1向量，对应p的分量依次是次数从高到底各多项式系数你现在浏览的是第四十三页，共78页用用Richard模拟模拟水稻叶伸长生长水稻叶伸长生长你现在浏览的是第四十四页，共78页关于关于inline函数函数例如：y=inline(

34、sin(x)-cos(x),x)输入y(0)，可得：-1作图：x=0:0.1:2*pi;plot(x,y(x)你现在浏览的是第四十五页，共78页1、插值问题：、插值问题：不知道某一函数不知道某一函数f(x)在待定范围在待定范围a,b上上 的具体表达式，而只能通过实验测量得到该的具体表达式，而只能通过实验测量得到该 函数在一系列点函数在一系列点ax1,x2,.,xn b上的值上的值 y0,y1,y2,.,yn，需要找一个简单的函数，需要找一个简单的函数P(x)来近似地代替来近似地代替f(x)，要求满足：，要求满足：P(xi)=yi (i=1,2,.,n)2、插值函数：、插值函数：P(x)，3、插

35、值法：求插值函数、插值法：求插值函数P(x)的方法的方法 插插 值值你现在浏览的是第四十六页，共78页二、常用插值函数1、多项式函数、多项式函数2、样条函数、样条函数你现在浏览的是第四十七页，共78页1、多项式插值方法（1）n次代数插值次代数插值（2）拉格朗日插值）拉格朗日插值几点说明：几点说明：（1）拉格朗日插值基函数仅与节点有关，而）拉格朗日插值基函数仅与节点有关，而 与被插值函数与被插值函数f(x)无关。无关。（2）拉格朗日插值多项式仅由数对）拉格朗日插值多项式仅由数对(xi,yi)(i 1,2,n)确定，而与数对排列次序无关。确定，而与数对排列次序无关。（3）多项式插值除了上述插值法外

36、还有其它）多项式插值除了上述插值法外还有其它 插值法，如插值法，如newton插值法、插值法、hermite插插 值法等。值法等。你现在浏览的是第四十八页，共78页2、样条插值方法（1）样条函数）样条函数m次半截幂函数次半截幂函数（2）k次次B样条或样条或k次基本样条函数的定义次基本样条函数的定义你现在浏览的是第四十九页，共78页（一）广泛使用的样条函数（一）广泛使用的样条函数（1）广泛采用：）广泛采用：二次样条、三次样条及二次样条、三次样条及B样条样条。（2）力学意义：）力学意义：A:二次样条在力学上解释为集中力偶作用二次样条在力学上解释为集中力偶作用 下的弹性细梁挠度曲线。下的弹性细梁挠度

37、曲线。B:弹性细梁受集中载荷作用形成的挠度曲弹性细梁受集中载荷作用形成的挠度曲 线，在小挠度的情况下，恰好表示为三线，在小挠度的情况下，恰好表示为三 次样条函数，集中载荷的作用点，恰好次样条函数，集中载荷的作用点，恰好 就是三次样条函数的节点。就是三次样条函数的节点。你现在浏览的是第五十页，共78页（1）二次样条的定义 设设a,b 的一个划分：的一个划分：a=x0 x1,x2,.,xn=b，函数，函数f(x)各节点的值分别为：各节点的值分别为：f(xi)=yi (i=1,2,.,n)如果二次样条函数：如果二次样条函数：满足：S(xi)=yi(i=1,2,.,n)你现在浏览的是第五十一页，共78

38、页（2）三次样条函数的定义 设设a,b 的一个划分：的一个划分：a=x0 x1,x2,.,xn=b，函数函数f(x)各节点的值分别为：各节点的值分别为：f(xi)=yi (i=1,2,.,n)如果三次样条函数：如果三次样条函数：3满足：S(xi)=yi(i=1,2,.,n)你现在浏览的是第五十二页，共78页例：某实验对一根长例：某实验对一根长10米的钢轨进行热源的温米的钢轨进行热源的温度传播测试。用度传播测试。用x表示测量点表示测量点0:2.5:10(米米)，用，用h表示测量时间表示测量时间0:30:60(秒秒)，用，用T表示测试所得各点表示测试所得各点的温度的温度()。试用线性插值求出在一分

39、钟内每隔。试用线性插值求出在一分钟内每隔20秒、钢轨每隔秒、钢轨每隔1米处的温度米处的温度TI。命令如下：命令如下：x=0:2.5:10;h=0:30:60;T=95,14,0,0,0;88,48,32,12,6;67,64,54,48,41;xi=0:10;hi=0:20:60;TI=interp2(x,h,T,xi,hi)你现在浏览的是第五十三页，共78页 例：设例：设z=x2+y2，对，对z函数在函数在0,10,2区域内进行插值。区域内进行插值。为了说明这个维数的插值，再考虑一个问题。设人们对为了说明这个维数的插值，再考虑一个问题。设人们对平板上的温度分布估计感兴趣，给定的温度值取自平板

40、表平板上的温度分布估计感兴趣，给定的温度值取自平板表面均匀分布的格栅。面均匀分布的格栅。采集了下列的数据：采集了下列的数据：width=1:5;%index for width of plate(i.e.,the x-dimension)depth=1:3;%index for depth of plate(i,e,the y-dimension)temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86%temperature data temps=82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 你现在浏览的

41、是第五十四页，共78页 如同在标引点上测量一样，矩阵如同在标引点上测量一样，矩阵temps表示整个平板表示整个平板的温度分布。的温度分布。temps的列与下标的列与下标depth或或y-维相联系，行与维相联系，行与下标下标width或或x-维相联系维相联系(见图见图6)。为了估计在中间点的。为了估计在中间点的温度，我们必须对它们进行辨识。温度，我们必须对它们进行辨识。wi=1:0.2:5;%estimate across width of plate d=2;%at a depth of 2 zlinear=interp2(width,depth,temps,wi,d);%linear int

42、erpolation zcubic=interp2(width,depth,temps,wi,d,cubic);%cubic interpolation plot(wi,zlinear,-,wi,zcubic)%plot results xlabel(Width of Plate),ylabel(Degrees Celsius)title(Temperature at Depth=num2str(d)你现在浏览的是第五十五页，共78页图图6 在深度在深度d=2处的平板温度处的平板温度你现在浏览的是第五十六页，共78页 另一种方法，我们可以在两个方向插值。先在三另一种方法，我们可以在两个方向插值

43、。先在三维坐标画出原始数据，看一下该数据的粗糙程度维坐标画出原始数据，看一下该数据的粗糙程度(见图见图7)。mesh(width,depth,temps)%use mesh plot xlabel(Width of Plate),ylabel(Depth of Plate)zlabel(Degrees Celsius),axis(ij),grid xi,yi=meshgrid(width,depth);zi=interp2(width,depth,temps,xi,yi,cubic);mesh(xi,yi,zi)你现在浏览的是第五十七页，共78页图图7 平板温度平板温度你现在浏览的是第五十八页

44、，共78页然后在两个方向上插值，以平滑数据然后在两个方向上插值，以平滑数据(见图见图8)。di=1:0.2:3;%choose higher resolution for depth wi=1:0.2:5;%choose higher resolution for width zcubic=interp2(width,depth,temps,wi,di,cubic);%cubic mesh(wi,di,zcubic)xlabel(Width of Plate),ylabel(Depth of Plate)zlabel(Degrees Celsius),axis(ij),grid你现在浏览的是第

45、五十九页，共78页图图8 二维插值后的平板温度二维插值后的平板温度你现在浏览的是第六十页，共78页 上面的例子清楚地证明了，二维插值上面的例子清楚地证明了，二维插值更为复杂，只是因为有更多的量要保持更为复杂，只是因为有更多的量要保持跟踪。跟踪。interp2的基本形式是的基本形式是interp2(x,y,z,xi,yi,method)。这里。这里x和和y是两个独立是两个独立变量，变量，z是一个应变量矩阵。是一个应变量矩阵。x和和y对对z的的关系是关系是 z(i,:)=f(x,y(i)和和 z(:,j)=f(x(j),y).也就是，当也就是，当x变化时，变化时，z的第的第i行与行与y的第的第i个

46、元素个元素y(i)相关，当相关，当y变化时，变化时，z的第的第j列与列与x的第的第j个元素个元素x(j)相关，。相关，。xi是沿是沿x-轴插轴插值的一个数值数组；值的一个数值数组；yi是沿是沿y-轴插值的一轴插值的一个数值数组。个数值数组。你现在浏览的是第六十一页，共78页 可选的参数可选的参数method可以是可以是 linear，cubic或或nearest。在这种情况下，。在这种情况下，cubic不意味着不意味着3次样条，而是使用次样条，而是使用3次多次多项式的另一种算法。项式的另一种算法。linear方法是线性插方法是线性插值，仅用作连接图上数据点。值，仅用作连接图上数据点。neare

47、st方法方法只选择最接近各估计点的粗略数据点。只选择最接近各估计点的粗略数据点。插值的优点：插值的优点：利用已知点确定未知点利用已知点确定未知点 粗糙粗糙 精确精确 集合大的集合大的 简化的简化的你现在浏览的是第六十二页，共78页例：某观测站测得某日例：某观测站测得某日6:00时至时至18:00时之间每隔时之间每隔2小时的室内外温度小时的室内外温度()，用，用3次样条插值分别求得次样条插值分别求得该日室内外该日室内外6:30至至17:30时之间每隔时之间每隔2小时各点的小时各点的近似温度近似温度()。设时间变量设时间变量h为一行向量，温度变量为一行向量，温度变量t为一为一个两列矩阵，其中第一列

48、存放室内温度，第二列个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。命令如下：储存室外温度。命令如下：h=6:2:18;t=18,20,22,25,30,28,24;15,19,24,28,34,32,30;XI=6.5:2:17.5YI=interp1(h,t,XI,spline)%用用3次样条插值计算次样条插值计算你现在浏览的是第六十三页，共78页 在讨论二维插值之前，强调在讨论二维插值之前，强调interp1所所强制的二个强约束是很重要的。强制的二个强约束是很重要的。首先，人们不能要求有独立变量范围首先，人们不能要求有独立变量范围以外的结果，例如，以外的结果，例如，interp1

49、(hours,temps,13.5)导致一个错误，因为导致一个错误，因为hours在在1到到12之间变化。之间变化。其次，独立变量必须是单调的。即独其次，独立变量必须是单调的。即独立变量在值上必须总是增加的或总是减立变量在值上必须总是增加的或总是减小的。小的。你现在浏览的是第六十四页，共78页 在我们的例子里，在我们的例子里，hours是单调的。然而，如果是单调的。然而，如果我们已经定义独立变量为一天的实际时间，我们已经定义独立变量为一天的实际时间，time_of_day=7:12 1:6%start at 7AM,end at 6PM time_of_day=7 8 9 10 11 12 1

50、 2 3 4 5 6 则独立变量将不是单调的，因为则独立变量将不是单调的，因为time_of_day增加到增加到12，然后跌到，然后跌到1，再然后增加。，再然后增加。如果用如果用time_of_day代替代替interp1中的中的hours，将会，将会返回一个错误。同样的理由，人们不能对返回一个错误。同样的理由，人们不能对temps插插值来找出产生某温度的时间值来找出产生某温度的时间(小时小时)，因为，因为temps不不是单调的。是单调的。你现在浏览的是第六十五页，共78页案例3估计水箱的水流量模型长度单位：长度单位：E(3024cm)容积单位：容积单位：G(=3.785L(升升)某些镇的用水