插值与数据拟合建模.pptx

插值与数据拟合就是通过一些已知数据去插值与数据拟合就是通过一些已知数据去确定某类函数的参数或寻找某个近似函数，使确定某类函数的参数或寻找某个近似函数，使所得的函数与已知数据具有较高的精度，并且所得的函数与已知数据具有较高的精度，并且能够使用数学分析的工具分析数据所反映的对能够使用数学分析的工具分析数据所反映的对象的性质象的性质 几种常用的方法：几种常用的方法：1 1、一般插值法、一般插值法 2 2、样条插值法、样条插值法 3 3、最小二乘曲线、最小二乘曲线 4 4、曲面的拟合、曲面的拟合数据拟合与插值建模数据拟合与插值建模 上大学二年级的小华正在做作业，上大学二年级的小华正在做作业，“爸爸，计爸爸，计算这道题要用到算这道题要用到sin sin ，可是我的计算器坏了，可是我的计算器坏了，怎么办。怎么办。”当工程师的老张从厚厚的一摞旧书底下当工程师的老张从厚厚的一摞旧书底下抽出一本数学用表来，抽出一本数学用表来，“给你，这是我念大学时用给你，这是我念大学时用的，那时候啊，计算器听都没听说过。的，那时候啊，计算器听都没听说过。”小华拿着表翻了一会儿，无奈的说：小华拿着表翻了一会儿，无奈的说：“表上每表上每 才有一个函数值，这里只才有一个函数值，这里只sin sin 和和sin “sin “表中没有的，都可以用插值方法计算表中没有的，都可以用插值方法计算”“”“插值！我们插值！我们的数学实验课就要学了，不过，今天我要先自己想个的数学实验课就要学了，不过，今天我要先自己想个办法，用这个算出办法，用这个算出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.5783-0.5760)0.6=0.5774 聪明的小华用的这个办法是一种插值聪明的小华用的这个办法是一种插值方法方法分段线性插值。实际上，分段线性插值。实际上，插值可插值可以理解为，要根据一个用表格表示的函数，以理解为，要根据一个用表格表示的函数，计算表中没有的函数值。计算表中没有的函数值。表中有的，如（表中有的，如（sin sin ，0.57600.5760）(sin (sin ,0.5783),0.5783)称为称为节点节点；要计算的，如；要计算的，如sin sin ，称为，称为插值点插值点，结果（，结果（0.57740.5774）即为）即为插值插值。小华作的线性函数为小华作的线性函数为插值函数插值函数，插值函数，插值函数所表示的直线当然要通过节点。所表示的直线当然要通过节点。插值最初来源于天体计算插值最初来源于天体计算由若干观由若干观测值（即节点）计算任意时刻星球的位置测值（即节点）计算任意时刻星球的位置（即插值点和插值）（即插值点和插值）的需要。现在，虽的需要。现在，虽然人们已很少需要用它从函数表计算函数值然人们已很少需要用它从函数表计算函数值了，但是插值仍然在诸如机械加工等工程技了，但是插值仍然在诸如机械加工等工程技术和数据处理等科学研究中有着许多直接的术和数据处理等科学研究中有着许多直接的应用，另一方面，插值又是数值微分、数值应用，另一方面，插值又是数值微分、数值积分、常微分方程数值等数值计算的基础。积分、常微分方程数值等数值计算的基础。几天后，小华在物理实验里又碰几天后，小华在物理实验里又碰到一个看起来非常类似的问题：有一到一个看起来非常类似的问题：有一只对温度敏感的电阻，已经测得了一只对温度敏感的电阻，已经测得了一组温度组温度T T和电阻和电阻R R数据。数据。现在想知道现在想知道60600C时的电阻多大。时的电阻多大。温度温度t(0C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032 小华征求老师的意见，老师给了他两点提示：小华征求老师的意见，老师给了他两点提示：一是在直角坐标系中把一是在直角坐标系中把5 5个点（个点（T T，R R）画一下，看）画一下，看看电阻看电阻R R和温度和温度T T之间大致有什么关系之间大致有什么关系；二是；二是测量测量数据总有相当大的误差数据总有相当大的误差，这与用函数表作插值计，这与用函数表作插值计算应该有不同之处吧（虽然函数表也存在舍入误算应该有不同之处吧（虽然函数表也存在舍入误差，但很小，可以认为表中数值是精确的）差，但很小，可以认为表中数值是精确的）通过图形小华看到，通过图形小华看到，R R与与T T大致呈直线关系，大致呈直线关系，于是用手画了一条靠近于是用手画了一条靠近5 5个点的直线，又想起中学个点的直线，又想起中学物理学过，金属材料的电阻率与温度成正比，从物理学过，金属材料的电阻率与温度成正比，从而确定而确定R R与与T T的关系应该是的关系应该是 R=at+bR=at+b 其中其中a a，b b为待定常数。为待定常数。正是由于测量误差的存在，由正是由于测量误差的存在，由R=at+bR=at+b表表示的直线不可能通过全部示的直线不可能通过全部5 5个点，所以，与插个点，所以，与插值曲线要通过全部节点不同，小华打算作一条值曲线要通过全部节点不同，小华打算作一条尽量靠近所有的点的直线，求出尽量靠近所有的点的直线，求出a a，b b待定常数，待定常数，由此计算由此计算t=t=60600C 的的R R就十分简单了。就十分简单了。根据一组（二组）数据，即平面上的根据一组（二组）数据，即平面上的若干点，确定一个一元函数，即曲线，若干点，确定一个一元函数，即曲线，使这些节点与曲线总体来说尽量接近，使这些节点与曲线总体来说尽量接近，这就是曲线拟合。这就是曲线拟合。函数值与曲线拟合都是要根据一组函数值与曲线拟合都是要根据一组数据构造一个函数据构造一个函 数作为近似，由于近似数作为近似，由于近似的要求不同，二者的数学方法是完全不的要求不同，二者的数学方法是完全不同的。同的。数数 据据 拟拟 合合1.拟合的基本原理；拟合的基本原理；2.最小二乘拟合；最小二乘拟合；3.用用Matlab作最小二乘拟合；作最小二乘拟合；4.4.如何用拟合解决实际问题。如何用拟合解决实际问题。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采集血样，测得血药浓度采集血样，测得血药浓度c如下表：如下表：求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).半对数坐标系半对数坐标系(semilogy)下的图形下的图形Log10c(t)=a t+b引例引例1 1：血药浓度的变化规律：血药浓度的变化规律曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法 已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi)i=1,n,寻求一个函数（曲线）寻求一个函数（曲线）y=f(x),使使 f(x)在某种准则下与所有数据点最为接近，即曲线在某种准则下与所有数据点最为接近，即曲线拟合得最好。拟合得最好。+xyy=f(x)(xi,yi)i i 为点为点（xi,yi)与与曲线曲线 y=f(x)的距离的距离最小二乘拟合最小二乘拟合 第一步：第一步：先选定一类函数先选定一类函数f(x，a1，a2，am)其其准则为（最小二乘准则）：使准则为（最小二乘准则）：使n个点个点（xi,yi)与与曲线曲线 y=f(x，a1，a2，am)的距离的距离 i 的平方和最小的平方和最小。其中其中a1,a2,am为待定常数。为待定常数。f f可以为一些简单的可以为一些简单的“基函数基函数”（如幂函数，三角函数等等）的线性组合：（如幂函数，三角函数等等）的线性组合：第二步：确定参数第二步：确定参数a1,a2,am，记记问题归结为，求问题归结为，求 a1,a2,am 使使 J(a1,a2,am)最小。最小。这样的拟合称为最小二乘拟合。这样的拟合称为最小二乘拟合。除了最小二乘准则（即各点误差的平方和除了最小二乘准则（即各点误差的平方和最小），你认为还可以用怎样的拟合准则？最小），你认为还可以用怎样的拟合准则？比较起来，最小二乘准则有什么优点？比较起来，最小二乘准则有什么优点？思考思考最小二乘拟合函数最小二乘拟合函数 f(xf(x，a a1 1，aam m)的选取的选取 +f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1exp(a2x)+f=a1exp(a2x)1.1.通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 f f；2.2.将数据将数据 (xi,yi)i=1,n 作图，通过直观判断确定作图，通过直观判断确定 f：2.2.作一般的作一般的最小二乘曲线拟合，可利用已有程序最小二乘曲线拟合，可利用已有程序curvefit,curvefit,其调用格式为：其调用格式为：a=curvefit(f,a0,x,y)1.1.作多项式作多项式f(x)=a1xm+amx+am+1函数函数拟合拟合,可利用已有程序可利用已有程序polyfit,polyfit,其调用格式为其调用格式为:a=polyfit(x,y,m)用用MATLABMATLAB作最小二乘拟合作最小二乘拟合数据点数据点拟合多项式次数拟合多项式次数系数系数注：注：f f为拟合函数为拟合函数y=f(a,x)y=f(a,x)的函数的函数MM文件，文件，f(a,x)f(a,x)为拟合函数。为拟合函数。数据点数据点待定常数待定常数a a的初值的初值函数函数M文件文件用用MATLABMATLAB作多项式最小二乘拟合作多项式最小二乘拟合2.2.用命令用命令polyfit(x,y,m)得到得到 a1=3.3940,a2=702.49181.选取函数选取函数 R=a1t+a2温度温度t(0C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032例例.由数据由数据拟合拟合R=f(t)用用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 0.3946 0.3903 0.5474 0.3459 0.1370 0.2211 0.1704 0.2636用命令用命令curvefit(f,a0,x,y)拟合的应用拟合的应用参数辨识参数辨识 数学建模的方法：机理分析和测试分析。数学建模的方法：机理分析和测试分析。机理分析是根据对客观事物特性的认识，找机理分析是根据对客观事物特性的认识，找出反映内部机理的数量规律，建立的模型常有明出反映内部机理的数量规律，建立的模型常有明确的物理意义。确的物理意义。测试分析将研究的对象看作一个测试分析将研究的对象看作一个“黑箱黑箱”，通过对实验数据的统计分析，找出与数据拟合得通过对实验数据的统计分析，找出与数据拟合得最好的模型。最好的模型。机理分析机理分析模型结构模型结构实验数据实验数据未知参数未知参数范例：薄膜渗透率的测定范例：薄膜渗透率的测定 一、问题：一、问题：某种医用薄膜，具有从高浓度的溶液向低某种医用薄膜，具有从高浓度的溶液向低浓度的溶液扩散的功能，在试制时需测定薄膜浓度的溶液扩散的功能，在试制时需测定薄膜被物质分子穿透的能力。被物质分子穿透的能力。测定方法：用面积为测定方法：用面积为S S的薄膜将容器分成的薄膜将容器分成体积分别为体积分别为 的两部份，在两部分中分的两部份，在两部分中分别注满该物质的两种不同浓度的溶液。此时该别注满该物质的两种不同浓度的溶液。此时该物质分子就会从高浓度溶液穿过薄膜向低浓度物质分子就会从高浓度溶液穿过薄膜向低浓度溶液中扩散。平均每单位时间通过单位面积薄溶液中扩散。平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，膜的物质分子量与膜两侧溶液的浓度差成正比，比例系数比例系数K K表征了薄膜被该物质分子穿透的能表征了薄膜被该物质分子穿透的能力，称为渗透率。定时测量容器中薄膜某一侧力，称为渗透率。定时测量容器中薄膜某一侧的溶液浓度，以此确定的溶液浓度，以此确定K K。VAVBS二、问题分析二、问题分析 考考察察时时段段tt，t+tt+t薄薄膜膜两两侧侧容容器器中中该该物物质质质量的变化。质量的变化。设设 ，对对容容器器的的B B部部分分溶溶液液浓度的测试结果如下表：（浓度单位浓度的测试结果如下表：（浓度单位 ）1）在容器的一侧，物质质量的增加是由于另一在容器的一侧，物质质量的增加是由于另一侧的物质向该侧渗透的结果，因此物质质量的增侧的物质向该侧渗透的结果，因此物质质量的增量应等于另一侧的该物质向这侧的渗透量。量应等于另一侧的该物质向这侧的渗透量。以容器以容器A A侧为例，在时段侧为例，在时段tt，t+tt+t物质质量物质质量的增量为：的增量为：分别表示在时刻分别表示在时刻t t膜两侧溶液膜两侧溶液设设的浓度，浓度单位的浓度，浓度单位:由于平均每单位时间通过单位面积薄膜的由于平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比，比物质分子量与膜两侧溶液的浓度差成正比，比例系数为例系数为K K。因此，在时段因此，在时段tt，t+tt+t，从，从B B侧渗透至侧渗透至A A侧的该物质的质量为：侧的该物质的质量为：于是有：于是有：两边除以两边除以tt，并令，并令t0t0取极限再稍加整理取极限再稍加整理即得：即得：分别表示在初始时刻两侧溶液的浓度分别表示在初始时刻两侧溶液的浓度其中其中（1）2)注意到整个容器的溶液中含有该物质的质量不注意到整个容器的溶液中含有该物质的质量不变变,与初始时刻该物质的含量相同，因此与初始时刻该物质的含量相同，因此 从而：从而：加上初值条件：加上初值条件：代入式（代入式（1）得：）得：便可得出便可得出CB(t)的变化规律，从而根据实验数据进行的变化规律，从而根据实验数据进行拟合，估计出参数拟合，估计出参数K,。三、数学模型三、数学模型假设：假设：1）薄膜两侧的溶液始终是均匀的；薄膜两侧的溶液始终是均匀的；2）平均每单位时间通过单位面积薄膜的物质分平均每单位时间通过单位面积薄膜的物质分子量与膜两侧溶液的浓度差成正比。子量与膜两侧溶液的浓度差成正比。3 3）薄膜是双向同性的即物质从膜的任何一侧向）薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的性能是相同的。另一侧渗透的性能是相同的。基于假设和前面的分析，基于假设和前面的分析，B B侧的浓度侧的浓度CB(t)应满足应满足如下微分方程和初始条件：如下微分方程和初始条件：四、求解方法：四、求解方法：1.函数拟合法函数拟合法前面得到的模型是一个带初值的一阶线性微分方前面得到的模型是一个带初值的一阶线性微分方程，解之得：程，解之得：问题归结为利用问题归结为利用C CB B在时刻在时刻t tj j的测量数据的测量数据C Cj j(j=1,2,.,N)(j=1,2,.,N)来辨识来辨识 K K 和和 。引入引入从而从而 用函数用函数CB(t)来拟合所给的实验数据，来拟合所给的实验数据，从而估计出其中的参数从而估计出其中的参数a,b,K。将将代入上式有：代入上式有：用用MATLABMATLAB软件进行计算软件进行计算.1 1）编写函数编写函数M-M-文件文件 nongdu.mnongdu.mfunction f=nongdu(x,tdata)function f=nongdu(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata);f=x(1)+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,0.05;x0=0.2,0.05,0.05;x=curvefit(nongdu,x0,tdata,cdata)x=curvefit(nongdu,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,进一步求得：进一步求得：2.非线性规划法非线性规划法 利用利用C CB B在时刻在时刻t tj j的测量数据的测量数据C Cj j(j=1,2,.,N)(j=1,2,.,N)来辨来辨识识 K K 和和 。问题可转化为求函数问题可转化为求函数即求函数即求函数的最小值点（的最小值点（K K，a a，b b）。）。3.导函数拟合法导函数拟合法前面得到的微分方程为：前面得到的微分方程为：令令上式变为：上式变为：这可以看作这可以看作随随CB的变化规律的变化规律(j=1,2,.,N)若知道一组数据若知道一组数据则可用最小二乘拟合的方法来求出函数则可用最小二乘拟合的方法来求出函数中的未知参数中的未知参数K和和h。即为求参数即为求参数K,a使下列误差函数达到最小：使下列误差函数达到最小：该问题等价于用函该问题等价于用函 数数 f(K,a,CB)=K(0.01a-0.02CB)来拟合数据来拟合数据(j=1,2,.,N)用用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.mbaomof.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)3)3)输出结果输出结果:x=0.1009 0.014:x=0.1009 0.014 即即 k=0.1009,h=0.014k=0.1009,h=0.014作函数拟合作函数拟合x0=0.2,0.1;x=curvefit(baomof,x0,cdata,dcj)4.线性化迭代法线性化迭代法前面带初始条件的一阶线性微分方程的解为前面带初始条件的一阶线性微分方程的解为其中：其中：如果得到了参数如果得到了参数K的一个较好的近似值的一个较好的近似值K*，则，则将将CB(t)关于关于K在在K*处展开，略去处展开，略去K的二次及以上的二次及以上的项得的项得CB(t)的一个近似式的一个近似式通过极小化确定a,b,d,再由K=d/0.02b得到K*的修正值K。K*K*-K,得到K的一个新的近似值，用同样的方法再求新的修正值K。这个过程可以不断重复，直到修正值足够小为止。1）当K的初值取为k=0.3时，出现奇异情况，迭代不收敛；2）当K的初值取为k=0.2时，经四次迭代，已经收敛到一个很好的解。迭代结果如下表。五、结果及误差分析五、结果及误差分析 几种方法得出的结果及相应的误差总结于几种方法得出的结果及相应的误差总结于下表，误差为计算数据与实验数据之差的平方下表，误差为计算数据与实验数据之差的平方和。和。注：导函数拟合法得出的参数值精度有限，线性化迭代法要求参数的初值比较接近精确值。因此可将导函数拟合法和线性化迭代法结合起来使用，把前者得到的参数K的值作为迭代法中K的初值，这样可使迭代法收敛或收敛更快。3）取K的初值为k=0.1009,只一次迭代就得到2）中的最后结果。函数拟合法的拟合效果函数拟合法的拟合效果求解参数辨识模型的方法：求解参数辨识模型的方法：函数拟合；函数拟合；非线性规划；非线性规划；导函数拟合；导函数拟合；线性化迭代；线性化迭代；其它方法其它方法。用用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.841.743.744.845.545.3生长观测记录数据模型表达式：模型表达式：程序！程序！关于关于polyfit命令命令命令：p=polyfit(x,y,n)（1）x与y为模拟数据（2）n为拟合多项式的次数（3）当n=1时为用最小二乘法进行直线拟合（4）得到的向量p为长度n+1向量，对应p的分量依次是次数从高到底各多项式系数用用Richard模拟模拟水稻叶伸长生长水稻叶伸长生长关于关于inline函数函数例如：y=inline(sin(x)-cos(x),x)输入y(0)，可得：-1作图：x=0:0.1:2*pi;plot(x,y(x)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、插值法：求插值函数、插值法：求插值函数P(x)的方法的方法 插插 值值二、常用插值函数1、多项式函数、多项式函数2、样条函数、样条函数1、多项式插值方法（1）n次代数插值次代数插值（2）拉格朗日插值）拉格朗日插值几点说明：几点说明：（1）拉格朗日插值基函数仅与节点有关，而）拉格朗日插值基函数仅与节点有关，而 与被插值函数与被插值函数f(x)无关。无关。（2）拉格朗日插值多项式仅由数对）拉格朗日插值多项式仅由数对(xi,yi)(i 1,2,n)确定，而与数对排列次序无关。确定，而与数对排列次序无关。（3）多项式插值除了上述插值法外还有其它）多项式插值除了上述插值法外还有其它 插值法，如插值法，如newton插值法、插值法、hermite插插 值法等。值法等。2、样条插值方法（1）样条函数）样条函数m次半截幂函数次半截幂函数（2）k次次B样条或样条或k次基本样条函数的定义次基本样条函数的定义（一）广泛使用的样条函数（一）广泛使用的样条函数（1）广泛采用：）广泛采用：二次样条、三次样条及二次样条、三次样条及B样条样条（2）力学意义：）力学意义：A:二次样条在力学上解释为集中力偶作用二次样条在力学上解释为集中力偶作用 下的弹性细梁挠度曲线。下的弹性细梁挠度曲线。B:弹性细梁受集中载荷作用形成的挠度曲弹性细梁受集中载荷作用形成的挠度曲 线，在小挠度的情况下，恰好表示为三次线，在小挠度的情况下，恰好表示为三次 样条函数，集中载荷的作用点，恰好就是三次样条样条函数，集中载荷的作用点，恰好就是三次样条函数的节点。函数的节点。（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)（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)例：某实验对一根长例：某实验对一根长10米的钢轨进行热源的温米的钢轨进行热源的温度传播测试。用度传播测试。用x表示测量点表示测量点0:2.5:10(米米)，用，用h表示测量时间表示测量时间0:30:60(秒秒)，用，用T表示测试所得表示测试所得各点的温度各点的温度()。试用线性插值求出在一分钟。试用线性插值求出在一分钟内每隔内每隔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)例：设例：设z=x2+y2，对，对z函数在函数在0,10,2区域内进行插值。区域内进行插值。为了说明这个维数的插值，再考虑一个问题。设人们对为了说明这个维数的插值，再考虑一个问题。设人们对平板上的温度分布估计感兴趣，给定的温度值取自平板平板上的温度分布估计感兴趣，给定的温度值取自平板表面均匀分布的格栅。表面均匀分布的格栅。采集了下列的数据：采集了下列的数据：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 如同在标引点上测量一样，矩阵如同在标引点上测量一样，矩阵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 interpolation 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)图图6 在深度在深度d=2处的平板温度处的平板温度 另一种方法，我们可以在两个方向插值。先另一种方法，我们可以在两个方向插值。先在三维坐标画出原始数据，看一下该数据的粗糙在三维坐标画出原始数据，看一下该数据的粗糙程度程度(见图见图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)图图7 平板温度平板温度然后在两个方向上插值，以平滑数据然后在两个方向上插值，以平滑数据(见图见图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图图8 二维插值后的平板温度二维插值后的平板温度 上面的例子清楚地证明了，二维插值上面的例子清楚地证明了，二维插值更为复杂，只是因为有更多的量要保持更为复杂，只是因为有更多的量要保持跟踪。跟踪。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个元素个元素y(i)相关，当相关，当y变化时，变化时，z的第的第j列与列与x的第的第j个元素个元素x(j)相关。相关。xi是沿是沿x-轴插值轴插值的一个数值数组；的一个数值数组；yi是沿是沿y-轴插值的一个轴插值的一个数值数组。数值数组。可选的参数可选的参数method可以是可以是 linear，cubic或或nearest。在这种情况下，。在这种情况下，cubic不意味着不意味着3次样条，而是使用次样条，而是使用3次多次多项式的另一种算法。项式的另一种算法。linear方法是线性插方法是线性插值，仅用作连接图上数据点。值，仅用作连接图上数据点。nearest方法方法只选择最接近各估计点的粗略数据点。只选择最接近各估计点的粗略数据点。插值的优点：插值的优点：利用已知点确定未知点利用已知点确定未知点 粗糙粗糙 精确精确 集合大的集合大的 简化的简化的例：某观测站测得某日例：某观测站测得某日6:00时至时至18:00时之间每时之间每隔隔2小时的室内外温度小时的室内外温度()，用，用3次样条插值分次样条插值分别求得该日室内外别求得该日室内外6:30至至17:30时之间每隔时之间每隔2小小时各点的近似温度时各点的近似温度()。设时间变量设时间变量h为一行向量，温度变量为一行向量，温度变量t为一为一个两列矩阵，其中第一列存放室内温度，第二个两列矩阵，其中第一列存放室内温度，第二列储存室外温度。命令如下：列储存室外温度。命令如下：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次样条插值计算次样条插值计算 在讨论二维插值之前，强调在讨论二维插值之前，强调interp1所所强制的二个强约束是很重要的。强制的二个强约束是很重要的。首先，人们不能要求有独立变量范围首先，人们不能要求有独立变量范围以外的结果，例如，以外的结果，例如，interp1(hours,temps,13.5)导致一个错误，因为导致一个错误，因为hours在在1到到12之间变化。之间变化。其次，独立变量必须是单调的。即独其次，独立变量必须是单调的。即独立变量在值上必须总是增加的或总是减立变量在值上必须总是增加的或总是减小的。小的。在我们的例子里，在我们的例子里，hours是单调的。然而，如是单调的。然而，如果我们已经定义独立变量为一天的实际时间，果我们已经定义独立变量为一天的实际时间，time_of_day=7:12 1:6%start at 7AM,end at 6PM time_of_day=7 8 9 10 11 12 1 2 3 4 5 6 则独立变量将不是单调的，因为则独立变量将不是单调的，因为time_of_day增加到增加到12，然后跌到，然后跌到1，再然后增，再然后增加。如果用加。如果用time_of_day代替代替interp1中的中的hours，将会返回一个错误。同样的理由，人们不能，将会返回一个错误。同样的理由，人们不能对对temps插值来找出产生某温度的时间插值来找出产生某温度的时间(小时小时)，因为因为temps不是单调的。不是单调的。案例3估计水箱的水流量模型长度单位：长度单位：E(3024cm)容积单位：容积单位：G(=3.785L(升升)某些镇的用水管理机构需估计公众的用水速度某些镇的用水管理机构需估计公众的用水速度(单位是单位是G/h)和每天总用水量的数据许多地方没有测量流入或流和每天总用水量的数据许多地方没有测量流入或流出水箱流量的设备，而只能测量水箱中的水位出水箱流量的设备，而只能测量水箱中的水位(误差不超过误差不超过5%)当水箱水位低于某最低水位当水箱水位低于某最低水位L时，水泵抽水，灌入水时，水泵抽水，灌入水箱，直至水位达到最高水位箱，直至水位达到最高水位H为止。但这也无法测量水泵为止。但这也无法测量水泵的流量，因此在水泵启动时不易建立水箱中水位和水泵工的流量，因此在水泵启动时不易建立水箱中水位和水泵工作时用水量之间关系。水泵一天灌水作时用水量之间关系。水泵一天灌水12次，每次约次，每次约2h。试估计在任意时刻试估计在任意时刻(包括水泵灌水期间包括水泵灌水期间)t流出水箱的流量流出水箱的流量f(t)，并估计一天的总用水量。，并估计一天的总用水量。案例1数据 表中给出了某镇中某一天的真实用水数据，表中测表中给出了某镇中某一天的真实用水数据，表中测量时间以量时间以s为单位，水位以为单位，水位以10-2 E为单位例如为单位例如3316s以以后，水箱中的水深降至后，水箱中的水深降至3110 E时，水泵自动启动把水时，水泵自动启动把水输入水箱；而当水位回升至输入水箱；而当水位回升至35.5E时水泵停止工作。时水泵停止工作。1、模型假设（1）1）除了问题中特别说明的数据以外，其它给）除了问题中特别说明的数据以外，其它给 定的数据其测量误差不超过定的数据其测量误差不超过52）影响水箱流量的唯一因素是该镇公众对水）影响水箱流量的唯一因素是该镇公众对水 的普通需要所给的数据反映该镇在通常的普通需要所给的数据反映该镇在通常 情况下一天的用水量，不包括任何非常情情况下一天的用水量，不包括任何非常情 况如水箱中水的短缺、水管破裂、自然况如水箱中水的短缺、水管破裂、自然 灾害等灾害等 一个物理定理一个物理定理 1、模型假设（2）3）水泵的灌水速度为常数，不随时间变化也）水泵的灌水速度为常数，不随时间变化也 不是巳灌水量的函数，因此假设水泵大约不是巳灌水量的函数，因此假设水泵大约 在水位在水位27E时开始灌水，在水位时开始灌水，在水位355E时时 停止灌水同时假设水泵不会损坏或不需停止灌水同时假设水泵不会损坏或不需 要维护要维护4）从水箱中流出的最大流速小于水泵的灌水）从水箱中流出的最大流速小于水泵的灌水 速度为了满足公众的用水需求不让水棺速度为了满足公众的用水需求不让水棺 中的水流尽中的水流尽 1、模型假设（3）5）