哈尔滨工程大学工程算法课件04插值方法

1、1第四章第四章 插值方法插值方法4.0 4.0 引言引言4.1 4.1 多项式插值问题的一般提法多项式插值问题的一般提法4.2 4.2 拉格朗日拉格朗日(Lagrange)(Lagrange)插值插值4.3 4.3 差商与差分及其性质差商与差分及其性质4.4 4.4 牛顿插值公式牛顿插值公式4.5 4.5 分段插值法分段插值法4.6 4.6 三次样条插值三次样条插值4.7 4.7 曲线拟合的最小二乘法曲线拟合的最小二乘法2引引 言言 1 插值法插值法是广泛应用于理论研究和生产实践的重要数值方法，是广泛应用于理论研究和生产实践的重要数值方法，它是用简单函数它是用简单函数(特别是多项式或分段多项式

2、特别是多项式或分段多项式)为各种离散数为各种离散数组建立连续模型；为各种非有理函数提供好的逼近方法。众所组建立连续模型；为各种非有理函数提供好的逼近方法。众所周知，反映自然规律的数量关系的函数有三种表示方法：周知，反映自然规律的数量关系的函数有三种表示方法：x= ysiny，(开普勒开普勒(Kepler)方程方程)。 325f xxxA 解析表达式解析表达式 。(1865年，瓦里斯年，瓦里斯Walis；1690年，年，Raphson拉夫逊；拉夫逊；1669年，牛顿年，牛顿Newton；历史悠久的方程；历史悠久的方程)。悬链线方程；悬链线方程； y = cos(x/)。3B 图像法图像法C 表格

3、法表格法2 2 事实上，许多数据都是用表格法给出的事实上，许多数据都是用表格法给出的( (如观测和实验而得到的函数如观测和实验而得到的函数数据表格数据表格) )，可是，从一个只提供离散的函数值去进行理论分析和进行，可是，从一个只提供离散的函数值去进行理论分析和进行设计，是极不方便的甚至是不可能的。因此需要设法去寻找与已知函数设计，是极不方便的甚至是不可能的。因此需要设法去寻找与已知函数值相符，并且形式简单的插值函数值相符，并且形式简单的插值函数( (或近似函数或近似函数) )。3 3 另外一种情况是，函数表达式完全给定，但其形式不适宜计算机使用，另外一种情况是，函数表达式完全给定，但其形式不适

4、宜计算机使用，因为计算机只能执行算术和逻辑操作，因此涉及连续变量问题的计算都因为计算机只能执行算术和逻辑操作，因此涉及连续变量问题的计算都需要经过离散化以后才能进行。如数值积分方法、数值微分方法、差分需要经过离散化以后才能进行。如数值积分方法、数值微分方法、差分方程以及有限元法等，都必须直接或间接地应用到插值理论和方法。方程以及有限元法等，都必须直接或间接地应用到插值理论和方法。44.1 4.1 多项式插值问题的一般提法多项式插值问题的一般提法 01n01nii01ny=f xabxxxabn+1f xyyyP xP x=y i=0,1,2n1xxx1P xy=f xab1假设函数是，上的实值

5、函数， ， ， ，是，上个互异的点，在这些点上的取值分别为， ， ， 。求已确定的函数，使之满足：， ，称 ， ， ，为，关系式称为，函数称为插值节点插值原则插值函数的，区间 ，函数插值区为。间称插值法的概念5 插值函数插值函数p(x)作为作为f(x)的近似，可以选择不同类型的函数，的近似，可以选择不同类型的函数，如如p(x)为代数多项式、三角多项式、有理分式；其函数性态为代数多项式、三角多项式、有理分式；其函数性态可以是光滑的、亦可以是分段光滑的。其中，可以是光滑的、亦可以是分段光滑的。其中，代数多项式代数多项式类类的插值函数占有重要地位的插值函数占有重要地位: (a) 结构简单、计算机容易

6、处理、任何多项式的导数和积分结构简单、计算机容易处理、任何多项式的导数和积分也易确定，并且仍是多项式。也易确定，并且仍是多项式。 (b) 著名的著名的Weierstrass逼近定理逼近定理(定义在闭区间上的任何定义在闭区间上的任何连续函数连续函数f(x)，存在代数多项式，存在代数多项式p(x)一致逼近一致逼近f(x)，并达，并达到所要求的精度到所要求的精度.)。 因此，我们主要考虑代数多项式的插值问题。因此，我们主要考虑代数多项式的插值问题。2 例题分析：例题分析： 已知函数已知函数f(x)有如下数据：有如下数据： 求求f(x)的插值多项式的插值多项式p(x)，并求，并求f(x)在在x=0.5

7、处的近似值。处的近似值。6 2340123401400123401234112340123423,000111224816120001234119310,0,424931424P xaa xa xa xa xaaaPafPaaaaafPaaaaafPafPaaaafaaaaaP xxx 解：设其中， ， ， ， 为待定系数。由给定的条件：联立上式得：于是40.50.50.3906xfP74.2 4.2 拉格朗日拉格朗日(Lagrange)(Lagrange)插值插值 01001 000101 1110101,0,0,11,11ijnnniinnnnnnnnnnnxxijnnP xaa xa x

8、P xyinP xaa xa xyP xaa xa xyP xaa xa xya aan 从如下数据表着手，并假定求 次多项式使得， ，。根据插值条件，有显然，这是一个关于的元线性方程组，其系数矩阵行多项式插值的存在唯一性：列式为：8定理：由定理：由n+1个不同插值节点个不同插值节点 可以惟一确定一个可以惟一确定一个 n 次多项式次多项式 满足插值条满足插值条件件 。从理论上说，由方程组。从理论上说，由方程组(1)可以求出可以求出 的惟一解的惟一解,从而确定从而确定 。但从数值计算上看，当。但从数值计算上看，当n 较大时较大时求解线性方程组的工作量较大且不便应用。为解决此问题，求解线性方程组的

9、工作量较大且不便应用。为解决此问题，现已提出了不少构造现已提出了不少构造 的巧妙办法。的巧妙办法。 0011011010ji n0111V,11,2,V,=01,nnnnnninnijnxxxxx xxxxx inx xxxxa aa注意到插值节点两两相异，而故方程组有唯一解于是满足插值条件的多项式存在且唯一。01,nxxx01,na aa 01nnnPxaa xa x iiP xy nPx nPx9 1011001100110110011110011101001010000101,;-,0,11;2iiilxlxinx x yyL xaa xL xy L xyL xxyxyyyL xxxxx

10、xxxxyx xxxyyl x yl xilxlagrange首先讨论时的情形已知求使得显然是过，和，两点的一条直线。由点斜式容易求得：其中具有如下特点：插值的基函数构造法01101100;1lxlxlx称其为线性插值基函数。10 102012210011122220011221200102021101201220212,;.,0,1inx x xyy yLxL xy L xy LxyLxxyxyxyxxxxlxxxxxxxxxlxxxxxxxxxlxxxxxl xi再讨论时的情形已知求使得显然是过，、 ，、 ，三点的一条抛物线。仿照线性插值基函数的构造方法，令其中具有如 0001021011

11、122021221;0;00;1;00;0;1lxlxlxlxlxlxlxlxlx下特点：称其为抛物线插值基函数。如图所示11 12200102021101220120202101,2,.,1,2,iiinniinjij iijjixxxxLxyxxxxxxxxyxxxxxxxxyl x yxxxxLxLxy innxxl xxxl xin于是，最后讨论一般情形求使得令 次多项式插值基函数为：具有如下特点：12 00021,0,1,0,nLagrangeLagrange1,0,1,2,3,4,5ijijnnnjniiiiij iijjniiijlxi jnijxxLxylx yxxLxlxxi

12、ilx 于是，满足插值条件的 次多项式可以直接写为：我们称为多项式，为插值基函数。思考 给定。下面哪个是的图像？13a) 已知连续函数已知连续函数f(x)的函数表如右：的函数表如右：求方程求方程f(x)=0在在(-1,2)内的近似根。内的近似根。解：利用解：利用Lagrange插值法有，解方程插值法有，解方程取初值取初值x=0.5，利用牛顿法求解可得，利用牛顿法求解可得f(x) 在在(-1,2) 内的近似根为内的近似根为0.67433 。3 3 例题例题 33201211222101 11201 01 021212 1111221 202115914126xxxxxxLxxx xxx xxxx

13、 325914120 xxx14b) 利用利用100，121的开方计算的开方计算 。解：由于解：由于利用利用Lagrange 插值法有插值法有于是，于是， 的准确值为的准确值为10.72380529，因此近似值，因此近似值10.71428 有有3位有效数字。位有效数字。 11211001011100121121100 xxLx1115121115100115115101110.71428100121121 100L115115154 4 插值余项插值余项如图所示，其截断误差如图所示，其截断误差 ，称为称为Lagrange插值多项式的余项。插值多项式的余项。 nnRxf xLx定理：假设定理：假

14、设f(x)在在a, b上有上有n+1阶导数，且在不同插值节阶导数，且在不同插值节点点 取值为取值为 是经过插值样点是经过插值样点 的的Lagrange插值多项式，若引进记号插值多项式，若引进记号:则有如下的误差估计：则有如下的误差估计：01,nx xx ,iinf xy Lx ,0,1,iix yin 1010()nnniiwxxxxxxxxx16 1011010010011 !,1 !00,1,2nnnniinnnnnnnniinnnniinfRxf xLxxxnfwxa bnRxf xLxinRxRxk xxxxxxxk xxxtf tLtk xtxtxtxf tLtk xtxtx x x

15、xn证明：因为于是可假设具有如下形式作辅助函数容易看出，有共个相异零点，且 1111110,1Pr1,|0nnnnnnnitnia bnRolleinciplettabnta bdfLk xtxdt在上存在阶导数。根据，在的两个零点之间至少有一个零点，故在，上至少有个零点。如此类推，在上至少有一个零点 ，使得17注意到Ln是n次多项式， 的首项为tn+1，故 100;nnniiLttx1101 !nninidtxndt由上述方程解得 1,1 !nfk xa bn于是 101 !nnniifRxxn注注1：如果 在(a,b)上有界，即 1nfx 10,( , )nMfxMxa b 则有余项估计：

16、 01 !nniiMRxxxn注注2：当f(x)为任一个次数n的多项式时，由 f(n+1)(x) 0，可知Rn(x)0,因此，插值多项式Ln(x)对于次数n的多项式的估计是精确的。18思考思考: :(1) (1) 线性插值的余项估计：线性插值的余项估计：(2) (2) 抛物线插值的余项估计：抛物线插值的余项估计： 0121018maxxx xxxR xfx 201201,6fRxxxxxxxx x n01f xLxf -1 .454520242584245245020405450 xx xxx xxlxxxxxxxlx 例1：已知函数y=f的观察数据为试构造的拉格朗日多项式并估计解：先构造基函

17、数例题19 2333032325254240452424245250543545245584402524515531243542142151552411421421721/ ,k kkxx xx xxlxxx xx xxlxx xxxxxLxy lxx xxx xxxxxLf xx 所求三次多项式为例 ：已知连续函数其函数表如下：求方 3fLagrange程的值并估计误差。解：利用插值法有20 22242422222.540.522.524240.42.522.5422.50.254242.50.050.4251.15330.32563max2,81 33333232.5346 80.0312

18、51333-0.3250.00833xxxLxxxxxxxfLfxxfxfRfLRfL 故因为故实际误差214.3 4.3 差商与差分及其性质差商与差分及其性质 1001101201012201010100001,)112nnnnf xf xf x xf xxxf x xf x xf x x xf xxxf xxf xxf x xxxxf xf xf xf xxan定义 ：称为函数的称为函数的二一般地，称为函数的；特别的，定义为函数关于的零阶差商。由此可知，高阶差商总是由比它低一阶的两个差商组合而成。一阶差商；性阶差商阶差商质；：差商的概念 差商性质010100111,nnknkkkkkkkn

19、nnyyyyf x xxxxxxxxxx阶差商可以表示成个函数值的线性组合，即220101010101100112012020112020110021221012010202101220210101021012.,1111f xf xyyf x xxxxxxxf x xf x xf x x xxxyyyyxxxxxxxxxxxxyyyxxxxxxxxxxxxxxyyxxxxxxxx例： 22021011001210202100101-1)2,1)3:-1,-2,-kyxxxxbf xxf xxf x x xf x x xf x x xcf xxnf xxxnf xx xxnf xkf xx x

20、xxn kknknk性质对称性 ：差商与节点的顺序无关。即，这一点可以从性质 看出。性质若是 的 次多项式，则一阶差商，是 的次多项式；二阶差商，是 的次多项式；一般地，函数的 阶差商，是 的次多项式，而时， 阶差商为零。233 3 利用差商表计算差商利用差商表计算差商利用差商的递推定义利用差商的递推定义,可以用递推来计算差商。可以用递推来计算差商。如下表：如下表：如要计算四阶差商，应再增加一个节点，表中还要增加一行。如要计算四阶差商，应再增加一个节点，表中还要增加一行。例例1：已知下表，计算三阶差商：已知下表，计算三阶差商 f 1 , 3 , 4 , 73244 差分的概念差分的概念2：设函

21、数：设函数y=f(x)在等距节点在等距节点 上的函数值上的函数值 ，其中，其中，h为常数称作步长。称为常数称作步长。称分别为分别为f(x)在在 处以处以h为步长的一阶向前差分，一阶为步长的一阶向前差分，一阶向后差分和一阶中心差分。称符号、向后差分和一阶中心差分。称符号、分别分别为向前差分算子，向后差分算子和中心差分算子。为向前差分算子，向后差分算子和中心差分算子。解：列表计算解：列表计算00,1,ixxih in iif xf111122/ 2/ 2iiiiiiiiiiiffffffff xhf xhffix25例已知函数例已知函数y=f(x)y=f(x)的数据如表中第的数据如表中第2 2、3

22、 3列。计算它的各阶差商。列。计算它的各阶差商。解：依据差商计算公式，结果列表中。解：依据差商计算公式，结果列表中。5 算例：算例：264.4 4.4 牛顿插值公式牛顿插值公式1 1 牛顿插值公式的构造牛顿插值公式的构造 Lagrange Lagrange插值虽然易算，但若要增加一个节点时，全部基函插值虽然易算，但若要增加一个节点时，全部基函数数l li i(x) (x) 都需重新算过。本节介绍另外一种方法都需重新算过。本节介绍另外一种方法牛顿插值法，牛顿插值法，并用它解决上面所述问题。并用它解决上面所述问题。 将将LagrangeLagrange插值公式插值公式 00001020101010

23、1000000001011, , , , ,0, ,nnnjniiiiij iijjnnnnnnx xL xyl x yxxN xaa x xa x xx xa x xx xa aaa aaf xf xf x xf xf xf x xx xx xf x xf x xf x x xx xL x其中，为待定系数。如何求？解：因为， 所以又中的改写成27 00101101012012201012012201000100,1,2,/-12110nnnnnnnnf x xf x xf x x xxxf x x xf x x xf x x x xxxf x x xf x x xf x x x xxxf x

24、 xxf xxf x xxxxf x xxf xxf x xxxxnnn有又一般地，将式代入， ，式代入式，式代入式， 00100120100100,11nnnnnnnnf xf xf x xxxf x x xxxxxf xxxxxxf x xxxxxxxRxnxnxnNxf xNxRx如此可得：尤为注意的是：最后一项中，差商部分含有 乃是余项部分，记作；而前面项中，差商部分都不含 ，因而前面项是关于 的 次多项式，记牛顿插值公作这就是。于是，上式式成为。282 2 算例算例 001001010110010000100100120101201220021011,n,Nnf xf xf x xx

25、xf x x xxxxxyyxf xf x xxxyxxxxf xf xf x xxxf x x xxxxxf x x x xxxxxxxxf xNxf x例 ：当时其中N这就是牛顿一次插值多项式，也就是点斜式直线方程。当 =2时这就是牛顿二次插值多项式。显然， 011010101220200101122021202011221Nxf xf xxxf xxxf xf xNxf xxxxxf xf xf xf xxxxxf xxxxxxx即满足二次插值条件29例例2：已知：已知求满足以上插值条件的牛顿型插值多项式。求满足以上插值条件的牛顿型插值多项式。解：由于解：由于则牛顿三次插值多项式为则牛顿

26、三次插值多项式为例例3：已知：已知f(x)在六个点的函数值如下表，运用牛在六个点的函数值如下表，运用牛顿型插值多项式求顿型插值多项式求f(0.596)的近似值。的近似值。00101201230,1,4,1.25;f xf x xf x x xf x x x x 3N014131.25134xxxxxxx3230 2001001201232012301234301,0.5960.410751.11600.1960.280.1960.0460.632010,0.5960.6320100.19700.1960.0460.0540.6319145,Nxf xf x xxxf x x xxxxxNNxN

27、xf x x x xxxxxxxNNxNxf x x 欲求，只需在之后再加一项： 2340123640,0.03440.1960.0460.0540.2043.4 100.63191450.00000340.63191791)n,0,1,3nnnknkknnx xxxxxxxxxxNxP xNxP xNxf xknP xNxf x x ，和均为 次多项式，且满足插值条件：由多项式的唯一性：因而，两个公式的余项是相等的，即拉格朗日插值与牛顿插值的比较 11,1 !nnnnfxxwxwxn2) 2) 当插值多项式从当插值多项式从n-1n-1次增加到次增加到n n次时，拉格朗日型插值必须重新次时，拉

28、格朗日型插值必须重新计算所有的基本插值多项式；而对于牛顿型插值，只需用表格再计计算所有的基本插值多项式；而对于牛顿型插值，只需用表格再计算一个算一个n n阶差商，然后加上一项即可。阶差商，然后加上一项即可。31插值节点为等距节点插值节点为等距节点 如如下下图：图： 牛顿插值公式牛顿插值公式设等距节点设等距节点 ，记记令令 例如例如(右右图图)x 在在 的中点时，的中点时， 。将牛顿插值公式中的差商用差分代替，而将牛顿插值公式中的差商用差分代替，而从而，牛顿插值公式在等距插值节点下的形式为：从而，牛顿插值公式在等距插值节点下的形式为：余项为余项为0,0,1,kxxkh kn4 4 等距牛顿插值公

29、式等距牛顿插值公式0kxxkh0,0,1,kknyf xknxx x 。当0,0 xxthtn 23,x x02.5xxh00kxxxthxkhtk h 20000111112!nnNxyt yt tyt ttnyn 11111111nnnnnR xfw xfht ttnnn！等距牛顿向前插值公式等距牛顿向前插值公式32上式为等距牛顿向前插值公式。上式为等距牛顿向前插值公式。下面来推导下面来推导等距牛顿向后插值公式等距牛顿向后插值公式： 2110 ,;111112!11141,1.5,2,2.5,3,2.2nn knn knnnnnnnnnxxxthntxxkh xxtk hNxyt yt t

30、yt ttnynRxfht ttnnyf xexf 令这时余项为！例 ：设插值节点为相应的函数值如下表，求33 2.2122000323032432.29.025011,0.5,2.21,2.4,12.218.872322!2.22.211122.42.4 12.420.742100.16623362.2 =2.2 +0.16623=9.038552.22.2kkfehxxxthtNyt yt tyNNt ttyNNNN 解：精确值此时故于是求时，在后加一项！所以求在 4043234111232.42.4 12.422.430.481464240.016182.2 =2.2 -0.01618=

31、9.02237R =0.15269R =0.01354R =0.00264t tttyNN 后加一项！所以，344.5 4.5 分段插值法分段插值法1 1 问题的提出问题的提出 在区间在区间a, ba, b上用插值多项式上用插值多项式P P逼近函数逼近函数f f 时，时，f f 和和P P在在每个节点上的差异每个节点上的差异( (理论上理论上) )应该为零。自然，我们期应该为零。自然，我们期望在一切中间点上也能很好地逼近望在一切中间点上也能很好地逼近f f，并且当插值点，并且当插值点增加时这种逼近效果应该越来越好。但上述的期望不增加时这种逼近效果应该越来越好。但上述的期望不可能实现的。当认识到

32、这一点时，在数学界曾引起强可能实现的。当认识到这一点时，在数学界曾引起强烈的震动。烈的震动。2 2 反例反例(A) (A) 狄里克莱狄里克莱(Dirichlet)(Dirichlet)函数函数如果取插值节点为有理数时：如果取插值节点为有理数时：如果取插值节点为无理数时：如果取插值节点为无理数时：插值多项式插值多项式P P的图象如左下：的图象如左下：( (逼近函数效果极差逼近函数效果极差) ) 01xQfxxQ 0P xxR 1P xxR 35(B) 龙格龙格(Runge)函数函数将区间将区间-5,5分成分成n 等分，以等分，以 表示取表示取n+1个等分的插值个等分的插值多项式，右上图给出了多项

33、式，右上图给出了 的图象：的图象：可以看出，随着插值节点数增加，插值多项式的次数也相应可以看出，随着插值节点数增加，插值多项式的次数也相应增加，而对于高次插值容易带来剧烈振荡，带来数值不稳定；增加，而对于高次插值容易带来剧烈振荡，带来数值不稳定；越靠近端点逼近的效果越差越靠近端点逼近的效果越差(Runge现象现象)。因此，既要增加插值结点，减小插值区间，又要不增加插值因此，既要增加插值结点，减小插值区间，又要不增加插值多项式的次数以减少误差，我们可以采用分段插值的办法。多项式的次数以减少误差，我们可以采用分段插值的办法。 3P x 10Px 2P x nP x363 3 分段线性插值分段线性插

34、值分段线性插值问题的提出分段线性插值问题的提出给定区间给定区间a, b，将其分割成，将其分割成 已知已知函数函数 在这些插值结点的函数值为在这些插值结点的函数值为 ，求一个分，求一个分段函数段函数 ，使其满足：，使其满足：(1)(2) 在每个区间在每个区间 上，上， 是个一次函数；是个一次函数；(3) 在区间在区间a, b上，上， 为连续函数。为连续函数。易知，易知，P(x)是个折线函数，在是个折线函数，在每个区间每个区间 上上 于是，于是， 在在a, b上是连续的，但其一阶导数是不上是连续的，但其一阶导数是不连续的连续的(即不光滑的即不光滑的)。01,naxxxb yf x,0,1,kkyf

35、 xkn P x,0,1, ;kkP xy kn1,kkx x 1111kkkkkkkkxxxxP xyyxxxx P x P x1,0,1,kkx xkn P x374 分段线性函数的基函数分段线性函数的基函数我们从整体上来构造分段线性函数的基函数。每个插值节点上所对应的插值函数：我们从整体上来构造分段线性函数的基函数。每个插值节点上所对应的插值函数： 101010011111111111111120,0,1,1,0,0,iikiiiiiiiiiiiiinnnnnnnnlxkilxkixxxxxxxlxxxxinxxxxxxxxxlxxxxxxxxxxxxxxxxlxxxx是分段线性函数；3

36、8 01111211,10,514.5,111115nk kkkkkkk kkkkkkkP xy lxx xlxlxP xy lxylxxxfk kxkxkP xyykkkkyxkyx 于是，注意：此表达式在区间上，只有是非零的，其它基函数均为零，即。已知函数y=f在区间上取等矩插值节点 如下表求区间上分段线性插值函数，并利用它求出近似值。解：在每个分段区间上，例题k于是39 10.50,10.520.211,20.230.122,30.140.0588233,40.0588250.0384644,54.50.058824.550.038464.540.048644.50.0470588235

37、294174.50.04762270 xxxxxxP xxxxxxxxxxPfn 实际值：当时， 321996104.50.04705882352941n 当时，,a bx1,kkxxx由此可见，对于光滑性要求不高的插值问题，分由此可见，对于光滑性要求不高的插值问题，分段线性段线性插值的效果非常好！计算也简单！插值的效果非常好！计算也简单！根据拉格朗日一根据拉格朗日一次插值函数的余项，可以得到分次插值函数的余项，可以得到分段线性插值函数的插值误差估计：段线性插值函数的插值误差估计：对，当对，当 时，时，40 21101,128maxmaxkkkkknxa bhR xfxxxxR xMhxxMf

38、xh 则其中于是可以加密插值结点，缩小插值区间，使 减小，从而减小插值误差。414.6 4.6 三次样条插值三次样条插值1 样条函数的概念：样条函数的概念：1：在：在a, b上取上取n+1个插值结点个插值结点 已知函数已知函数 在这在这n+1个点的函数值为个点的函数值为 则则在在a, b上函数上函数 的的m次样条插值函数次样条插值函数 满足：满足：(1) 在在(a, b)上直到上直到m-1阶导数连续；阶导数连续；(2) (3) 在区间在区间 上，上， 是是m次多项式。次多项式。 yfx,kkyf x1,0,1,1kkx xkn 高次插值函数可以保证曲线的光滑性，但计算量大，有剧高次插值函数可以

39、保证曲线的光滑性，但计算量大，有剧烈振荡，数值稳定性差；而分段线性插值在分段点上仅连续而烈振荡，数值稳定性差；而分段线性插值在分段点上仅连续而不光滑不光滑( (导数不连续导数不连续) )，这往往不能满足工程设计的需要。样条函，这往往不能满足工程设计的需要。样条函数可以同时解决这两个问题，使插值函数既是低阶分段函数，数可以同时解决这两个问题，使插值函数既是低阶分段函数，又是光滑的函数。又是光滑的函数。01,naxxxb yf x S x S x,0,1, ;kkS xy kn S x422 2 三次样条函数三次样条函数 1,1,01 200 ,00 ,001,2,12,0,1, ;3,0,1,1

40、kkkkkkkkkka byf xS xS xa bS xS xS xS xSxSxknS xy knx xknS x在上函数的三次样条插值函数满足：在上 、 、阶导数连续；即在区间上，是三次多项式。43 1111111111,0,1, ,3kkkkkkkkkkkkkkkkkkkkkkkkkkkSxm kn mS xSxxxxxx xSxmmxxxxxxxxhxxSxmmhhx xS xyS xy由二阶导数连续，设是未知、待定的数。因为是分段三次多项式，则是分段一次多项式在每个区间内，记则 将上式在区间上积分两次，并且由三次样条函数的计算 11,kkkxx xS x来确定两个积分常数。当时，利

41、用一阶导数连续的性质，对上式求导，得： 332211111S6666kkkkkkkkkkkkkkkkxxxxhxxhxxxmmymymhhhh44 22111111111111k1S226,S063kk 1,S x,S063S0 =S0mkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkxxxxhxmmmmyyhhhhhyyxmmhxxhhyyxmmhxx 在上式中，令x=x 得：将上式中的 换成 - 得在上表达式，用x=x 代入，而联立上述两式，得到关于的方程11111111111111111+1,2,1636662kkkkkkkkkkkkkkkkkkkkkkkkkkkkk

42、kkkhhhhyyyymmmknhhhhhhyyyymmmhhhhhhhh：两边乘以，得：45111111111111101,6,1,21,2,111kkkkkkkkkkkkkkkkkkkkkkkkkkkkkknmm myyyhhyyyyChhhhhhhhmmmCknnmmmn 上式中，等式左边含未知量等式右边是已知的，令则得：这是含有个未知量， ，共有个方程组成的线性方程组。欲确定方程的解，尚 00000221,;2,;30nnnnnSxy SxySxy SxySxSx缺 个方程。因此，求三次样条函数还要 个附加条件。通常有如下三种给法：给出边界端点的一阶导数值：三转角方程。给出边界端点的二

43、阶导数值：三弯矩方程。给出导数值：；自然样条。46 00122111100010001011111,12260361,634nnkkkkkkkkkkkkkknnnnnnnnnS xyS xyxx xxxxxhS xmmmmyyhhhkxxhhyyymmhknxxhhyyymmh 给出边界端点的一阶导数值： 利用前面已推导的公式，当时，取，得：取得：移三转角方程的计算项得47100100001111010101121121111112216262222221222nnnnnnnnnnnnnnnnnnyymmyChhyymmyChhmmCmmmCmmmCmmC于是，我们可以建立如下方程组：其系数矩

44、阵是严格对角占优的三角矩阵;112n48 011332111211,S6666Hnkkkkkkkkkkkkkkkkkkm mmxxxxxxhxxxmmymhhhhxxymhermite从而可以解出。解出后可以得到三次样条函数的分段表达式，即当x时，注：三次样条插值与插值的区别。4900112110212232323343232212121111223322-12-2222222225nnnnnnnnnnnnnnnnnSxmSxmmmCmmmmCmmmCmmmCmmCm附加条件为 ，。则方程组为：其系数矩阵为三弯矩方程的计算这是一个三对角矩阵，由于这是一个三对角矩阵，由于 ，因而它是严格对角占优

45、的。原方程组是个因而它是严格对角占优的。原方程组是个三对角方程组，可以用追赶法求解。三对角方程组，可以用追赶法求解。12kk506 6 算例算例例：已知例：已知y=f(x)的函数值如下所示，求函数的三次样的函数值如下所示，求函数的三次样条插值。条插值。012111222102110110322122121121,2,1,112221,1,1;1233123366433 13,12216624435,121222/32/32hhhyyyyChhhhyyyyChhhhmm 解：建立方程组123539,44mm 解得 1332111211,S6666kkkkkkkkkkkkkkkkkkxxxxxxh

46、xxxmmymhhhhxxymh当x时，1,kkkhxx11111116,1,kkkkkkkkkkkkkkkkkkhhyyyyChhhhhhhh 1121,2,1kkkkkkmmmCkn51 013332123322,1,2213120106646113113731641884,2,4423923436246246422946xx xxxxS xxxxxxx xxxxS x 从而得到函数的三次样条插值；当时，当时，322137142884xxxx 52 233332323232,4,55491950464664114345103203361888137112884137124884345103

47、33 45888xx xxxxS xxxxxxxxxS xxxxxxxxx 当时，所以534.7 4.7 曲线拟合的最小二乘法曲线拟合的最小二乘法1 拟合问题的数学提法拟合问题的数学提法通过观测、测量或试验得到某一函数在通过观测、测量或试验得到某一函数在 的的函数值函数值 。我们可以用插值的方法对这一函数进行近。我们可以用插值的方法对这一函数进行近似，而插值方法要求所得到的插值多项式经过已知的这似，而插值方法要求所得到的插值多项式经过已知的这n个插值个插值结点；在结点；在n比较大的情况下，插比较大的情况下，插值值多项式往往是高次多项式，这多项式往往是高次多项式，这也就容易出现振荡也就容易出现振

48、荡现象；现象；虽然在插值结点上没有误差，但在插值虽然在插值结点上没有误差，但在插值结点结点之外插值之外插值误差变得很大，从整体上看，插值逼近效误差变得很大，从整体上看，插值逼近效果将变得果将变得很很差。于是，我们采用曲线拟合的方法。差。于是，我们采用曲线拟合的方法。所谓曲线拟合所谓曲线拟合是求一个是求一个简单的函数简单的函数 ，例如，例如 是是一个低次多一个低次多项式，这儿不要求项式，这儿不要求 通过已知的这通过已知的这n个点，而是个点，而是要求在整体上要求在整体上“尽量好尽量好”的逼近原函的逼近原函数数.这时，在每个已这时，在每个已知点上就会有误差知点上就会有误差 数数据拟合就据拟合就是从整

49、体上使误差是从整体上使误差 ，尽量的尽量的小小一些。一些。12,nx xx12,ny yy x yx,1,2, ,kkyxkn,1,2, ,kkyxkn54如果要求如果要求 达到最小，达到最小，因误差因误差 可正可负，本可正可负，本来很大的误差可能会正负抵消，来很大的误差可能会正负抵消，这样的提法不合理。为防止正这样的提法不合理。为防止正负抵消，可以要求负抵消，可以要求 达达到最小，但是由于绝对值函数到最小，但是由于绝对值函数不可以求导，分析起来不方不可以求导，分析起来不方便，求解也很难。为了既能防止正负抵消，又能便，求解也很难。为了既能防止正负抵消，又能便于我们分析、求解，提出如下问题：便于

50、我们分析、求解，提出如下问题：求一个低次多项式求一个低次多项式 ，使得，使得 达到最小，此问题便是一个曲线拟合的最小二乘达到最小，此问题便是一个曲线拟合的最小二乘问题问题。1nkkkyxkkyx1nkkkyx x21nkkkQyx552 2 直线拟合（一次函数）直线拟合（一次函数）a) a) 问题的提法问题的提法通过观测、测量或试验得到某一函数在通过观测、测量或试验得到某一函数在 的的 函数值函数值 即得到即得到n n组数据组数据 ，如果这些数据在直角坐标系中近似地分布在一条直线上，我们可以如果这些数据在直角坐标系中近似地分布在一条直线上，我们可以用直线拟合的方法。用直线拟合的方法。问题：问题

51、：已知数据已知数据 ，求一个一次多，求一个一次多项式项式 ( (实际上，就是求实际上，就是求a, b)a, b)，使得，使得 达到最小。注意到达到最小。注意到 中，中， 均是已知的，而均是已知的，而a, ba, b是未知量是未知量, ,是未知量是未知量a, ba, b的二元的二元函数，利用高等数学中求二元函数极小值函数，利用高等数学中求二元函数极小值( (最小值最小值) )的方法，因此，的方法，因此，上述问题转化为求解下列方程组上述问题转化为求解下列方程组12,ny yy12,nx xx 1122,nnx yx yxy 1122,nnx yx yxy xabx2211,nnkkkkkkQ a

52、byxyabx,Q a b,kkxy,0,0Q a baQ a bb562111111112111),20,20nkkknkkknkkkknnnkkkkknnkkkknnnkkkkkkkbQ a byabxQ a byabxaQ a byabxxbanabxbxnaxbyxaxbx y 正则方程组由得：因为得到如下的正则方程组 这是个关于这是个关于a, ba, b的二元一次方程组，称其为最小二乘问题的正则方的二元一次方程组，称其为最小二乘问题的正则方程组。解得程组。解得a, ba, b，便得到最小二乘问题的拟合函数，便得到最小二乘问题的拟合函数 y = a + bx y = a + bx。57

53、c) c) 算例算例例例1 1：已知：已知1010对数据如下表，利用最小二乘法求拟合对数据如下表，利用最小二乘法求拟合曲线曲线y = a + bxy = a + bx。58 10502550269.12109.946.4383,0.78776.43830.78772,1,2,?11 iiiiiiabababyxx yimxP xaxbxyP xaxbyxyxyabxx yx ya 形成正则方程组解得 于是，最小二乘拟合一次函数为例 ：已知如图所示如何构造拟合曲线分析：由于。于是可令：由此。这是个线性问题。将化为后易解 和b。59牛顿牛顿(Newton)插值多项式插值多项式 ).()(.)()(

54、)(110102010nnnxxxxxxaxxxxaxxaaxN 的系数的系数 可根据插值条件推出可根据插值条件推出, 即由即由 有有 naaa,.,10nixfxNiin,.,1 , 0)()()()(000 xfaxNn)()()(101101xfxxaaxNn)()()()(21202201102xfxxxxaxxaaxNn )().()(.)()(1100110nnnnnnnnxfxxxxxxaxxaaxN这是关于这是关于 的下三角方程组的下三角方程组, ,可以求得可以求得 naaa,.,1060)(00 xfa 10010101011,)()()()()(xxfxxxfxfxxaxf

55、a2101210201202021022,)(,)()()()(xxxfxxxxfxxfxxxxxxaxfxfa一般，用数学归纳法可证明一般，用数学归纳法可证明 ),.,1 , 0(,.,10nkxxxfakk所以所以n n次牛顿次牛顿(Newton)(Newton)插值公式为插值公式为 ).()(.,.)(,)()(110100100nnnxxxxxxxxxfxxxxfxfxN其余项其余项 0101( ), .,()().()nnnR xf x xxxxxxxxx61为牛顿插值多项式的误差。由插值多项式的为牛顿插值多项式的误差。由插值多项式的存在惟一性定理知，满足同一组插值条件的存在惟一性定

56、理知，满足同一组插值条件的拉格朗日插值多项式拉格朗日插值多项式P(x)P(x)与牛顿插值多项式与牛顿插值多项式N Nn n(x)(x)实际上是同一个多项式，仅是同一插值实际上是同一个多项式，仅是同一插值多项式的不同表达形式而已，因此得到牛顿多项式的不同表达形式而已，因此得到牛顿插值多项式的误差与拉格朗日插值多项式的插值多项式的误差与拉格朗日插值多项式的误差也完全相等。故有误差也完全相等。故有(1)0100( )( ), .,()()(1)!nnnnniiiifR xf x xx xx xx xn(1)01( ), .,(1)!nnff x x xxn 可以看出，牛顿插值公式计算方便，可以看出，

57、牛顿插值公式计算方便，增加一个插值点，只要多计算一项，而增加一个插值点，只要多计算一项，而Nn(x)的各项系数恰好是各阶差商值，的各项系数恰好是各阶差商值，很有规律很有规律62000,xxxfxfxxf fx0,x(x- x0)= f(x) - f(x0)f(x)+ fx0,x(x- x0)=f(x0)101001,xxxxfxxfxxxf fx1,x0,x(x-x1) =fx0,x-fx1,x0fx0,x+ fx1,x0,x(x-x1)= fx1,x0f(x)+ (x- x0) fx1,x0=f(x0)+ (x- x0) (x-x1) fx1,x0,x牛顿插值公式牛顿插值公式(另一种推导方法

58、）另一种推导方法）63f(x)=f(x0)+(x- x0)fx1,x0+(x- x0)(x-x1)fx1,x0,x201201012,xxxxxfxxxfxxxxf fx1,x0,x= (x-x2) fx2,x1,x0,x+fx2,x1,x0f(x)=f(x0)+(x- x0)fx1,x0 + (x- x0)(x-x1)fx2,x1,x0 + (x- x0)(x-x1)(x-x2) fx2,x1,x0,x64,.,).()(,.,).()(.,)(,)()()(011001110012100100 xxxxfxxxxxxxxxfxxxxxxxxxfxxxxxxfxxxfxfnnnnnn Nn(

59、x)Rn(x)其中其中Nn(x)称为称为牛顿插值多项式牛顿插值多项式 Rn(x)称为称为牛顿插值余项牛顿插值余项65xifxifxi,xi+1fxi,xi+1,xi+2 fxi,xi+1,xi+2x0f(x0)x1f(x1)fx0,x1x2f(x2)fx1,x2fx0,x1,x2x3f(x3)fx2,x3fx1,x2,x3fx0,x1,x2 ,x3,.,).()(.,)(,)()()(01110012100100 xxxfxxxxxxxxxfxxxxxxfxxxfxNnnnn66定义 函数y= f(x)在区间xi ,xi+1上的平均变化率iiiiiixxxfxfxxf111)()(,自变量之差

60、和因变量之差之比叫自变量之差和因变量之差之比叫差商差商 称为称为f(x)关于关于xi , xi+1 的一阶差商的一阶差商,并记为并记为fxi ,xi+1 二阶差商二阶差商iiiiiiiiixxxxfxxfxxxf212121,01102110,.,.,.,xxxxxfxxxfxxxfmmmmm阶差商阶差商67fxi,xj,xk是指fxi , xj , xk=fxj , xk- fxi , xj xk- xi一般的一般的,可定义区间可定义区间xi, xi+1 , xi+n上的上的n阶差商为阶差商为ininiiiniiiniiixxxxxfxxxfxxxf ,.,.,.,1121102102121