数据拟合与插值1课件.ppt

1、数据拟合与插值数据拟合与插值山东建筑大学山东建筑大学贺长伟贺长伟一、引言一、引言 在解决实际问题的生产（或工程）实践和科学在解决实际问题的生产（或工程）实践和科学实验过程中，通常需要通过研究某些变量之间的函实验过程中，通常需要通过研究某些变量之间的函数关系来帮助我们认识事物的内在规律和本质属性，数关系来帮助我们认识事物的内在规律和本质属性，而这些变量之间的未知函数关系又常常隐含在从试而这些变量之间的未知函数关系又常常隐含在从试验、观测得到的一组数据之中。因此，能否根据一验、观测得到的一组数据之中。因此，能否根据一组试验观测数据找到变量之间相对准确的函数关系组试验观测数据找到变量之间相对准确的函

2、数关系就成为解决实际问题的关键。就成为解决实际问题的关键。例如在工程实践和科学实验中，常常需要从一例如在工程实践和科学实验中，常常需要从一组试验观测数据组试验观测数据 之中找到自之中找到自变量变量x x与因变量与因变量y y 之间的函数关系，一般可用一个之间的函数关系，一般可用一个近似函数近似函数y y=f f(x x)来表示。函数来表示。函数y y=f f(x x)的产生的产生办法因观测数据和要求不同而异，通常可采用数据办法因观测数据和要求不同而异，通常可采用数据拟合与函数插值两种办法来实现。拟合与函数插值两种办法来实现。二、拟合和插值二、拟合和插值 数据拟合数据拟合主要是考虑到观测数据受随

3、机观测误差的主要是考虑到观测数据受随机观测误差的影响，进而寻求整体误差最小、能较好反映观测数影响，进而寻求整体误差最小、能较好反映观测数据的近似函数据的近似函数y y=f f(x x)，此时，此时并不要求并不要求所得到的所得到的近似函数近似函数y y=f f(x x)满足满足 函数插值函数插值则要求近似函数则要求近似函数y y=f f(x x)在每一个观测在每一个观测点点 处处一定要满足一定要满足 在这种情况下，通常要求观测数据相对比较准确，在这种情况下，通常要求观测数据相对比较准确，即不考虑观测误差的影响。即不考虑观测误差的影响。ix 在实际问题中，通过观测数据能否正确揭示在实际问题中，通过

4、观测数据能否正确揭示某些变量之间的关系，进而正确认识事物的内在某些变量之间的关系，进而正确认识事物的内在规律与本质属性，往往取决于规律与本质属性，往往取决于两方面因素两方面因素。其一是观测数据的准确性或准确程度，这是其一是观测数据的准确性或准确程度，这是因为在获取观测数据的过程中一般存在随机测量因为在获取观测数据的过程中一般存在随机测量误差，导致所讨论的变量成为随机变量。误差，导致所讨论的变量成为随机变量。其二是对观测数据处理方法的选择，即到底其二是对观测数据处理方法的选择，即到底是采用插值方法还是用拟合方法，插值方法之中、是采用插值方法还是用拟合方法，插值方法之中、拟合方法之中又选用哪一拟合

5、方法之中又选用哪一种插值或拟合技巧来处种插值或拟合技巧来处理观测数据。理观测数据。三、数据拟合三、数据拟合 在科学计算中经常要建立实验数据的数学模型。在科学计算中经常要建立实验数据的数学模型。给定函数的实验数据，需要用比较简单和合适的给定函数的实验数据，需要用比较简单和合适的函数来函数来逼近逼近（或拟合）实验数据。这种逼近的特（或拟合）实验数据。这种逼近的特点是：点是：(a)(a)适度的精度是需要的；适度的精度是需要的；(b)(b)实验数据有小的误差；实验数据有小的误差；(c)(c)对于某些问题，可能有某些特殊的信息能够对于某些问题，可能有某些特殊的信息能够用来选择实验数据的数学模型。用来选择

6、实验数据的数学模型。数据拟合基本方法数据拟合基本方法 逼近离散数据的基本方法就是曲线拟合，常采用逼近离散数据的基本方法就是曲线拟合，常采用最小二乘拟合最小二乘拟合。曲线拟合问题的数学描述是，已知一组（二维）曲线拟合问题的数学描述是，已知一组（二维）数据数据 即平面上的即平面上的n n个点个点 互不相同。寻求一个函数（曲线）互不相同。寻求一个函数（曲线）y y=f f(x x)，使，使f f(x x)在某种准则下与所有数据点在某种准则下与所有数据点最为接近，即曲线拟合得最好。最为接近，即曲线拟合得最好。最小二乘拟合分为最小二乘拟合分为线性最小二乘拟合和非线性线性最小二乘拟合和非线性最小二乘拟合。

7、最小二乘拟合。ix3.1 线性最小二乘拟合线性最小二乘拟合 线性最小二乘法是解决曲线拟合问题最常用的方线性最小二乘法是解决曲线拟合问题最常用的方法，基本思路是，令法，基本思路是，令 （1.11.1）其中其中 是一组事先选定的线性无关是一组事先选定的线性无关 的函数，的函数，是一组待定系数。寻求系数是一组待定系数。寻求系数 使得使得 与与 的距离的距离 的的平方平方和最小和最小。这种准则称为最小二乘准则，其求系数。这种准则称为最小二乘准则，其求系数 的方法称为线性最小二乘拟合方法的方法称为线性最小二乘拟合方法。iy()if x1 1、系数系数 的求法的求法若记若记 （1.21.2）则则J J 为

8、为 的二次函数。由数学分析（或高的二次函数。由数学分析（或高等数学）的极值理论，等数学）的极值理论，J J 达到最小的充分必要条件达到最小的充分必要条件是是 满足满足于是得到求使于是得到求使J J 达到最小的达到最小的 方法方法1,.maa1,.maa1,.maa 是求解线性方程组（称为法方程组）是求解线性方程组（称为法方程组）（1.31.3）即求解线性方程组即求解线性方程组 （1.41.4）若记若记 （1.51.5）则方程组（则方程组（1.41.4）可表示为）可表示为 （1.61.6）由于当由于当 线性无关时，线性无关时，R R列满秩，列满秩，可逆，所以方程组（可逆，所以方程组（1.61.6

9、）有唯一解）有唯一解 （1.71.7）2 2、函数组函数组 的选取的选取 面对一组数据面对一组数据 ，用线，用线性最小二乘法作曲线拟合时，首要的、也是关键性最小二乘法作曲线拟合时，首要的、也是关键的一步是恰当地选取的一步是恰当地选取 如果通过机理分析，能够知道如果通过机理分析，能够知道y y 与与x x之间应该之间应该有什么样的函数关系，则有什么样的函数关系，则 容易确定。容易确定。若无法知道若无法知道y y 与与x x之间的关系，通常可以将之间的关系，通常可以将数据数据 作图作图，直观地判断应该，直观地判断应该用什么样的曲线去作拟合。用什么样的曲线去作拟合。人们常用的曲线有人们常用的曲线有（

10、1 1）直线直线（2 2）多项式多项式（3 3）双曲线双曲线（一支）（一支），拟合前作变量，拟合前作变量 替换替换t=1/x,t=1/x,求解求解a1,a2a1,a2较简单。较简单。（4 4）指数曲线指数曲线 ，拟合前作变量代换，拟合前作变量代换 z=lny,tz=lny,t=1/x,=1/x,则指数曲线转化为线性函数则指数曲线转化为线性函数 ,求解求解a1,a2a1,a2较简单。较简单。小结小结 在实际计算过程中，面对一组已知数据，到在实际计算过程中，面对一组已知数据，到底用什么样的曲线拟合最好，可以在直观判断的底用什么样的曲线拟合最好，可以在直观判断的基础上，选几种曲线分别拟合，然后比较，

11、看哪基础上，选几种曲线分别拟合，然后比较，看哪条曲线的最小二乘指标条曲线的最小二乘指标J J 最小。最小。3.2 非线性最小二乘拟合非线性最小二乘拟合 非线性最小二乘法是假设非线性最小二乘法是假设f f(x x)是待定系数是待定系数的任意非线性函数，在最小二乘准则下求其系数的任意非线性函数，在最小二乘准则下求其系数 例如上述人们常用的双曲线和指数曲线就是非线例如上述人们常用的双曲线和指数曲线就是非线性最小二乘拟合中最常用的非线性函数，只不过在性最小二乘拟合中最常用的非线性函数，只不过在上面使用中将它们转变成线性最小二乘拟合方法。上面使用中将它们转变成线性最小二乘拟合方法。对于给定的实验数据，通

12、常应根据实验数据的走对于给定的实验数据，通常应根据实验数据的走向、趋势选择合适的数学模型，即拟合函数。向、趋势选择合适的数学模型，即拟合函数。例如当实验数据具有单调性和凸性时，可选择例如当实验数据具有单调性和凸性时，可选择下述适当的数学模型下述适当的数学模型y y=f f(x x)来拟合实验数据。来拟合实验数据。等，其中等，其中a a、b b 为参数。为参数。在有可能的情况下，一般将非线性拟合函数转在有可能的情况下，一般将非线性拟合函数转化为线性拟合函数求解，这一方面是如此求解简单，化为线性拟合函数求解，这一方面是如此求解简单，另一方面也是因为一般情况下求解法方程组另一方面也是因为一般情况下求

13、解法方程组 得到的通常仅是得到的通常仅是J J 的驻点，不一的驻点，不一定是极值点。也可以直接解定是极值点。也可以直接解J J 极小化问题。极小化问题。3.3 最小二乘拟合法的最小二乘拟合法的Matlab实现实现1.解方程组方法解方程组方法例例1 1 用最小二乘法求一个形如用最小二乘法求一个形如 的经验公式，使它与下表所示的数据拟合的经验公式，使它与下表所示的数据拟合。x 19 25 31 38 44 y 19.0 32.3 49.0 73.3 97.8解解 编写程序如下编写程序如下 x=19 25 31 38 44;y=19.0 32.3 49.0 73.3 97.8;r=ones(5,1)

14、,x.2;ab=ry x0=19:0.1:44;y0=ab(1)+ab(2)*x0.2;plot(x,y,o,x0,y0,r)2yabx2.线性最小二乘拟合（多项式拟合）方法线性最小二乘拟合（多项式拟合）方法 在线性最小二乘拟合中，用的较多的是多项在线性最小二乘拟合中，用的较多的是多项式拟合。如果取式拟合。如果取 即用即用m m 次多项式拟合给定数据，则次多项式拟合给定数据，则MatlabMatlab中有现中有现成的函数成的函数 a=polyfit(x0,y0,m)a=polyfit(x0,y0,m)其中输入参数其中输入参数x0,y0 x0,y0为要拟合的数据，为要拟合的数据，m m为拟合多为

15、拟合多项式的次数，输出参数项式的次数，输出参数a a为拟合多项式为拟合多项式 系数系数 多项式在多项式在x x处的值处的值y y可用下面的函数计算可用下面的函数计算 y=polyval(a,xy=polyval(a,x)10.mmya xa xa10,.,maaa a将例将例1 1中的数据拟合成一个二次方程中的数据拟合成一个二次方程 x=19 25 31 38 44;x=19 25 31 38 44;y=19.0 32.3 49.0 73.3 97.8;y=19.0 32.3 49.0 73.3 97.8;ab=polyfit(x,y,2)ab=polyfit(x,y,2)x0=19:0.1:

16、44;x0=19:0.1:44;y0=ab(3)+ab(2)y0=ab(3)+ab(2)*x0+ab(1)x0+ab(1)*x0.2;x0.2;plot(x,y,o,x0,y0,r)plot(x,y,o,x0,y0,r)结果：结果：abab=0.0497 0.0193 0.6882 0.0497 0.0193 0.6882例例2 2 某乡镇企业某乡镇企业1990-19961990-1996年的生产利润如下表：年的生产利润如下表：年份年份 1990 1991 1992 1993 1994 1995 19961990 1991 1992 1993 1994 1995 1996利润（万元）利润（万元

17、）70 122 144 152 174 196 20270 122 144 152 174 196 202试预测试预测19971997年和年和19981998年的利润。年的利润。解解 作已知数据的的散点图，作已知数据的的散点图，x0=1990 1991 1992 1993 1994 1995 1996;x0=1990 1991 1992 1993 1994 1995 1996;y0=70 122 144 152 174 196 202;y0=70 122 144 152 174 196 202;plot(x0,y0,plot(x0,y0,*)发现该乡镇企业的年生产利润几乎直线上升。发现该乡镇企

18、业的年生产利润几乎直线上升。因此，我们可以用因此，我们可以用 作为拟合函数来预测作为拟合函数来预测该乡镇企业未来的年利润。编写程序如下：该乡镇企业未来的年利润。编写程序如下：x0=1990 1991 1992 1993 1994 1995 1996;x0=1990 1991 1992 1993 1994 1995 1996;y0=70 122 144 152 174 196 202;y0=70 122 144 152 174 196 202;a=polyfit(x0,y0,1)a=polyfit(x0,y0,1)y97=polyval(a,1997)y97=polyval(a,1997)y98

19、=polyval(a,1998)y98=polyval(a,1998)求得求得19971997年的生产利润年的生产利润y97=233.4286y97=233.4286，19981998年的生产利年的生产利润润y98=253.9286y98=253.9286。10ya xa3 非线性最小二乘拟合非线性最小二乘拟合 MatlabMatlab的优化工具箱中提供了两个求非线性最的优化工具箱中提供了两个求非线性最小二乘拟合的函数：小二乘拟合的函数：curvefitcurvefit和和leastsqleastsq。使用这。使用这两个命令时，都要先建立两个命令时，都要先建立M M文件文件fun.mfun.m

20、，但它们定义，但它们定义f f(x x)的方式是不同的。的方式是不同的。1.curvefit1.curvefit 设已知设已知 xdataxdata=(xdata1,xdata2,=(xdata1,xdata2,xdata,xdatan n ),),ydata ydata=(ydata1,ydata2,=(ydata1,ydata2,ydata,ydata n n),),curvefit curvefit用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数 中的参变量中的参变量x x（向量），使得（向量），使得 最小。最小。输入格式为：输入格式为：1.x=curvefit(f

21、un,x0,xdata,ydata)1.x=curvefit(fun,x0,xdata,ydata)；2.x=curvefit(fun,x0,xdata,ydata,options)2.x=curvefit(fun,x0,xdata,ydata,options)；3.x=curvefit(fun,x0,xdata,ydata,options,grad)3.x=curvefit(fun,x0,xdata,ydata,options,grad)；4.x4.x，options=curvefit(fun,x0,xdata,ydata,options=curvefit(fun,x0,xdata,ydat

22、a,)5.x,options,funval=curvefit(fun,x0,xdata,ydata,5.x,options,funval=curvefit(fun,x0,xdata,ydata,)6.x,options,funval,Jacob=curvefit(fun,x0,xdata,ydat6.x,options,funval,Jacob=curvefit(fun,x0,xdata,ydata,a,).).输出目标函数值格式：输出目标函数值格式：f=fun(x,xdataf=fun(x,xdata).).其中其中x0 x0为迭代初值，为迭代初值，optionsoptions为控制参数。为

23、控制参数。2.leastsq2.leastsq设已知设已知 xdataxdata=(xdata1,xdata2,=(xdata1,xdata2,xdata,xdatan n ),),ydata ydata=(ydata1,ydata2,=(ydata1,ydata2,ydata,ydata n n),),leastsq leastsq 用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数 中的参变量中的参变量x x（向（向 量），使得量），使得 最小。其中最小。其中输入格式为：输入格式为：（1 1）x=leastsqx=leastsq(fun,x0)(fun,x0)；（2 2

24、）x=leastsqx=leastsq(fun,x0,options)(fun,x0,options)；（3 3）x=leastsqx=leastsq(fun,x0,options,grad)(fun,x0,options,grad)；（4 4）xx，options=leastsqoptions=leastsq(fun,x0,(fun,x0,)；（5 5）xx，options,funvaloptions,funval=leastsq(fun,x0,=leastsq(fun,x0,)例例3 3 用下面一组数据拟合函数用下面一组数据拟合函数 中的参数中的参数a a,b b,k k。_ 100 20

25、0 300 400 500 600 700 800 900 1000100 200 300 400 500 600 700 800 900 1000 4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59 4.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59_ 解解 该问题即解最优化问题该问题即解最优化问题1.1.用命令用命令curvefitcurvefit。此时。此时（1 1）编写编写M M文件文件curvefun1.mcurvefun1.mfunction f=curvefun1(x,tdata)fu

26、nction f=curvefun1(x,tdata)f=x(1)+x(2)f=x(1)+x(2)*exp(-0.02exp(-0.02*x(3)x(3)*tdatatdata)%)%其中其中x(1)=ax(1)=a；x(2)=bx(2)=b；x(3)=kx(3)=k；（2 2）输入命令输入命令tdatatdata=100:100:1000=100:100:1000cdatacdata=1e-=1e-0303*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.5

27、0,6.59;x0=0.2,0.05,0.005;x0=0.2,0.05,0.005;x=curvefit(curvefun1,x0,tdata,cdata)x=curvefit(curvefun1,x0,tdata,cdata)f=curvefun1(x,tdata)f=curvefun1(x,tdata)（3 3）运算结果为：）运算结果为：x=0.0070-0.0030 0.1012x=0.0070-0.0030 0.1012f=Columns 1 through 7f=Columns 1 through 70.0045 0.0050 0.0054 0.0057 0.0059 0.0061

28、 0.00630.0045 0.0050 0.0054 0.0057 0.0059 0.0061 0.0063Columns 8 through 10Columns 8 through 100.0064 0.0065 0.00660.0064 0.0065 0.0066即拟合得即拟合得a=0.0070a=0.0070，b=-0.0030b=-0.0030，k=0.0066k=0.0066。2.2.用命令用命令leastsqleastsq。此时此时x=(ax=(a，b b，k)k)（1 1）编写编写M M文件文件curvefun2.mcurvefun2.mfunction f=curvefun2

29、(x)function f=curvefun2(x)tdatatdata=100:100:1000;=100:100:1000;cdatacdata=1e-03=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;f=cdata-x(1)-x(2)f=cdata-x(1)-x(2)*exp(-0.02exp(-0.02*x(3)x(3)*tdatatdata)%)%其中其中x(1)=ax(1)=a；x(2)=b;x(3)=kx(2)=b

30、;x(3)=k；（2 2）输入命令输入命令x0=0.2,0.05,0.005;x0=0.2,0.05,0.005;x=leastsq(curvefun2,x0)x=leastsq(curvefun2,x0)f=curvefun2(x)f=curvefun2(x)（3 3）运算结果为：）运算结果为：x=0.0070-0.0030 0.1012x=0.0070-0.0030 0.1012f=1.0e-005f=1.0e-005*Columns 1 through 7Columns 1 through 70.0221 0.2081-0.3933-0.2872 0.2973 0.3561 0.0693

31、0.0221 0.2081-0.3933-0.2872 0.2973 0.3561 0.0693Columns 8 through 10Columns 8 through 10-0.2327-0.0970 0.0296-0.2327-0.0970 0.0296可以看出，两个命令的计算结果是相同的。可以看出，两个命令的计算结果是相同的。四、插值方法四、插值方法 在实际问题中所遇到的插值问题一般分为在实际问题中所遇到的插值问题一般分为 一维插值问题一维插值问题和二维插值问题。和二维插值问题。一维插值问题的一维插值问题的数学描述数学描述为：已知某一函数为：已知某一函数 y y=g g(x x)（g

32、g(x x)的解析表达式可能十分复杂也的解析表达式可能十分复杂也 可以是未知的）在区间可以是未知的）在区间 a a,b b 上上n n+1+1个互异点个互异点 处的函数值处的函数值 ,j j=0,1,.,=0,1,.,n n，还知道，还知道g g(x x)在在 a a,b b 上有若干阶导数，如何求出上有若干阶导数，如何求出g g(x x)在在 a a,b b 上任一点上任一点x x的近似值。的近似值。一维插值方法的一维插值方法的基本思想基本思想是：根据是：根据g g(x x)在区间在区间 a a,b b 上上n n+1+1个互异点个互异点 （称为节点）的函数（称为节点）的函数 值值 ，j j

33、=0,1,.,=0,1,.,n n，求一个足够光滑、简单，求一个足够光滑、简单 便于计算的函数便于计算的函数f f(x x)（称为插值函数）作为（称为插值函数）作为 g g(x x)的近似表达式，使得的近似表达式，使得 （4.14.1）然后计算然后计算f f(x x)在区间在区间 a a,b b(称为插值区间称为插值区间)上点上点x x（称为插值点）的值作为原函数（称为插值点）的值作为原函数g g(x x)（称为被插（称为被插函数）在此点的近似值。求插值函数函数）在此点的近似值。求插值函数f f(x x)的方法的方法称为插值方法，（称为插值方法，（4.14.1）称为插值条件。）称为插值条件。4

34、.1 常用插值方法常用插值方法1 1、代数多项式代数多项式比较简单，常用多项式作为插值函数。比较简单，常用多项式作为插值函数。假设假设f f(x x)是一个满足插值条件（是一个满足插值条件（4.14.1）的次数不）的次数不 超过超过n n的代数多项式，即的代数多项式，即 （4.24.2）为满足（为满足（4.14.1）的插值函数，则）的插值函数，则f f(x x)的的n n+1+1个待个待 定系数定系数 满足满足。（4.34.3）记此方程组的系数矩阵为记此方程组的系数矩阵为A A，则，则 是范德蒙是范德蒙(Vandermonde(Vandermonde)行列式。当行列式。当 互不相同时，此行列式

35、值不为零。因此方程组互不相同时，此行列式值不为零。因此方程组 （4.34.3）有唯一解。这表明，只要）有唯一解。这表明，只要n n+1+1个插值节个插值节 点点 互异，满足插值条件（互异，满足插值条件（4.14.1）的插）的插值多项式（值多项式（4.24.2）存在唯一。）存在唯一。从几何上看，上述多项式插值就是过从几何上看，上述多项式插值就是过n n+1+1个数据个数据 点点 作一条多项式曲线作一条多项式曲线 y y=f f(x x)近似原有曲线近似原有曲线y y=g g(x x)。二二 拉格朗日插值法拉格朗日插值法1 1 拉格朗日插值的基本结果拉格朗日插值的基本结果：用多项式函数（：用多项式

36、函数（4.24.2）作为插值函数时，希望通过解方程组（作为插值函数时，希望通过解方程组（4.34.3）而得）而得到待定系数到待定系数 的做法当的做法当n n比较大时是比较大时是不现实的。方便的方法是先构造一组不现实的。方便的方法是先构造一组基函数基函数 显然显然 是是n n次多项式且满足次多项式且满足 称称n n次多项式次多项式 为节点为节点 上的上的n n次次 插值基函数。令插值基函数。令 则多项式则多项式 显然满足插值条件显然满足插值条件 ，从而，从而 为插值多项式。插值多项式为插值多项式。插值多项式 又称又称 为为n n次拉格朗日（次拉格朗日（LagrangeLagrange）插值多项式

37、，由）插值多项式，由 方程组（方程组（4.34.3）解的存在唯一性，）解的存在唯一性，n n+1+1个互异节个互异节点点 的的n n次拉格朗日插值多项式次拉格朗日插值多项式 存在唯一。存在唯一。2 2 拉格朗日插值的拉格朗日插值的MatlabMatlab实践实践 MatlabMatlab中没有现成的拉格朗日插值函数，必须编写一个中没有现成的拉格朗日插值函数，必须编写一个M M文件实文件实现拉格朗日插值。设现拉格朗日插值。设n n个节点数据以数组个节点数据以数组x x0,0,y y0 0输入（注意输入（注意MatlabMatlab的数组下标从的数组下标从1 1开始），开始），m m 个插值点以数

38、组个插值点以数组x x输入，输出数输入，输出数组组y y 为为m m 个插值。编写一个名为个插值。编写一个名为lagrange.mlagrange.m的的M M文件：文件：function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:m z=x(i);s=0.0;for k=1:n p=1.0;for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j);end end s=p*y0(k)+s;end y(i)=s;end三三 分段线性插值法分段线性插值法 假设区间假设区间 a a,b b 上的连续函数上的连

39、续函数g g(x x)在在n n+1+1个节点个节点 上的函数值上的函数值 则得到则得到xyxy平面上的平面上的n n+1+1个数据点个数据点 连接相邻数据点连接相邻数据点 得到得到n n条线条线 段，它们组成一条折线。把区间段，它们组成一条折线。把区间 a a,b b 上这上这n n条条 折线段表示的函数称为被插函数折线段表示的函数称为被插函数g g(x x)关于这关于这 n n+1+1个数据点的分段线性插值函数，记作个数据点的分段线性插值函数，记作I I(x x)，则则I I(x x)具有如下具有如下性质性质：1 1）I I(x x)可以分段表示，在每个小区间可以分段表示，在每个小区间上，

40、它是线性函数，即上，它是线性函数，即 2 2）3 3）I I(x x)在在 a a,b b 上连续上连续。若构造插值基函数若构造插值基函数 则则 用用I I(x x)计算计算x x点的插值时，只用到点的插值时，只用到x x左右的两个左右的两个点，计算量与节点个数点，计算量与节点个数n n+1+1无关。无关。但但n n越大，分越大，分段越多，插值误差越小段越多，插值误差越小。分段线性插值的分段线性插值的MatlabMatlab实现实现用用MatlabMatlab实现分段线性插值不需要编制函数程序，实现分段线性插值不需要编制函数程序，MatlabMatlab中有现成的一维插值函数中有现成的一维插值

41、函数interp1interp1。y=interp1(x0,y0,x,method)y=interp1(x0,y0,x,method)methodmethod指定插值的方法，默认为线性插值。其值可为：指定插值的方法，默认为线性插值。其值可为：nearest nearest 最近项插值最近项插值linear linear 线性插值线性插值splinespline 立方样条插值立方样条插值cubic cubic 立方插值。立方插值。所有的插值方法要求所有的插值方法要求x0 x0是单调的。是单调的。当当x0 x0为等距时可以用快速插值法，使用快速插值法为等距时可以用快速插值法，使用快速插值法的格式为

42、的格式为*nearestnearest、*linear linear、*splinespline 、*cubiccubic。二维插值的二维插值的MatlabMatlab实现实现MatlabMatlab中有一些计算二维插值的程序。如中有一些计算二维插值的程序。如z=interp2(x0,y0,z0,x,y,method)z=interp2(x0,y0,z0,x,y,method)其中其中x0,y0 x0,y0分别为分别为m m 维和维和n n维向量，表示节点，维向量，表示节点，z0z0为为N N*m m 维矩阵，表示节点值，维矩阵，表示节点值，x,yx,y为一维数组，表示插为一维数组，表示插值点

43、，值点，x x与与y y应是方向不同的向量，即一个是行向量，应是方向不同的向量，即一个是行向量，另一个是列向量，另一个是列向量，z z为矩阵，表示得到的插值，为矩阵，表示得到的插值，methodmethod的用法同上面的一维插值。的用法同上面的一维插值。如果是如果是三次样条插值三次样条插值，可以使用命令，可以使用命令pp=csape(x0,y0,z0,conds,valconds)pp=csape(x0,y0,z0,conds,valconds)y=ppval(pp,x,yy=ppval(pp,x,y)其中其中x0,y0 x0,y0分别为分别为m m 维和维和n n维向量，维向量，z0z0为为

44、m m*n n 维矩维矩阵，具体使用方法同一维插值。阵，具体使用方法同一维插值。例例2 2 在一丘陵地带测量高程，在一丘陵地带测量高程，x x和和y y 方向每隔方向每隔100100米米测一个点，得高程如下表，试拟合一曲面，确定测一个点，得高程如下表，试拟合一曲面，确定合适的模型，并由此找出最高点和该点的高程。合适的模型，并由此找出最高点和该点的高程。x/y x/y 100 200 300 400 500100 200 300 400 500100 636 697 624 478 450100 636 697 624 478 450200 698 712 630 478 420200 698

45、712 630 478 420300 680 674 598 412 400300 680 674 598 412 400400 662 626 552 334 310400 662 626 552 334 310解解 编写程序如下：编写程序如下：clear,clcclear,clcx=100:100:500;x=100:100:500;y=100:100:400;y=100:100:400;z=636 697 624 478 450z=636 697 624 478 450698 712 630 478 420698 712 630 478 420680 674 598 412 400680

46、 674 598 412 400662 626 552 334 310;662 626 552 334 310;pp=csape(x,y,zpp=csape(x,y,z)cz1=fnval(pp,100:10:500,100:10:400)cz1=fnval(pp,100:10:500,100:10:400)cz2=interp2(x,y,z,100:10:500,100:10:400,spline)cz2=interp2(x,y,z,100:10:500,100:10:400,spline)四 三次样条插值法 1 1 三次样条的基本结果三次样条的基本结果 分段线性插值函数在节点处的一阶导数一

47、般不存分段线性插值函数在节点处的一阶导数一般不存在，在，光滑性不高光滑性不高，这影响了它在诸如机械加工等领，这影响了它在诸如机械加工等领域（希望插值曲线光滑）中的应用域（希望插值曲线光滑）中的应用。在船舶、飞机等设计中，面对已知的一些数据在船舶、飞机等设计中，面对已知的一些数据点点 ，绘图员的做法是：，绘图员的做法是：首先将这些数据点描绘在平面图纸上，再把一根富首先将这些数据点描绘在平面图纸上，再把一根富有弹性的细直条（称为样条）弯曲，使其一边通过有弹性的细直条（称为样条）弯曲，使其一边通过这些数据点，用压铁固定细直条的形状，沿样条边这些数据点，用压铁固定细直条的形状，沿样条边沿绘出一条光滑曲

48、线。往往要用几根样条分段完成沿绘出一条光滑曲线。往往要用几根样条分段完成上述工作，这时应当让样条在连接处也保持光滑。上述工作，这时应当让样条在连接处也保持光滑。根据力学理论进行分析，这样画出的曲线，在相邻根据力学理论进行分析，这样画出的曲线，在相邻两数据点之间实际上是次数不高于两数据点之间实际上是次数不高于3 3的多项式，在整的多项式，在整个区间上有连续的曲率。对绘图员用样条画出的曲个区间上有连续的曲率。对绘图员用样条画出的曲线进行数学模拟，导出了样条函数的概念线进行数学模拟，导出了样条函数的概念。定义于区间定义于区间 a a,b b 上的分段函数上的分段函数S S(x x)，若满足，若满足(

49、1)(1)S S(x x)在每一小区间在每一小区间 上是一个三次多上是一个三次多 项式函数；项式函数；(2)(2)S S(x x)在整个区间在整个区间 a a,b b 上有连续的二阶导数。上有连续的二阶导数。则称则称S S(x x)为为 a a,b b 上关于剖分上关于剖分 的一个三次样条函数。的一个三次样条函数。2 三次样条插值的Matlab实现MatlabMatlab中三次样条插值也有现成的函数：中三次样条插值也有现成的函数：y=interp1(x0,y0,x,spline)y=interp1(x0,y0,x,spline)；y=spline(x0,y0,x)y=spline(x0,y0,

50、x)；pp=csape(x0,y0,conds),pp=csape(x0,y0,conds),pp=csape(x0,y0,conds,valconds)pp=csape(x0,y0,conds,valconds)y=ppval(pp,xy=ppval(pp,x)。其中其中x0,y0 x0,y0是已知数据点，是已知数据点，x x是插值点，是插值点，y y是插值是插值点点的函数值。的函数值。对于三次样条插值，我们提倡使用函数对于三次样条插值，我们提倡使用函数csape,csapecsape,csape的返回值是的返回值是pppp形式，要求插值点的函形式，要求插值点的函数值，必须调用函数数值，必须