第二章(函数插值)

1、1计算方法电子教案函数插值姜 潮湖南大学现代车身先进设计制造国家重点实验室2 插值法插值法1 引言引言2 拉格朗日插值多项式拉格朗日插值多项式3 牛顿插值多项式牛顿插值多项式4 埃尔米特插值埃尔米特插值5 分段低次插值分段低次插值6 三次样条插值三次样条插值31 引引 言言1. 1插值问题的提法插值问题的提法 在生产和科研中出现的函数是多种多样的。常遇到这种情况：在某个实际问题中，虽然可以断定所考虑的函数 在区间 上存在且连续，但却难以但却难以找到它的解析表达式找到它的解析表达式，只能通过实验和观测得到在有限个点上的函数值（即一张函数表）。显然，要利用这张函数表来分析函数的性态、甚至直接求出其

2、 xfba ,实际问题无实际问题无法解析表达法解析表达4它一些点上它一些点上的函数值是非常困难的。在有些情况下，虽然可以写出函数 的解析表达式，但由于结构相当复杂，使用起来很不方便。面对这些情况，总希望根据所得函数表（或结构复杂的解析表达式），构造某个简单函数P(x)作为 的近似。 xf xfBLACK BOXINPUTOUTPUT有限元无网格刚体动力学.各类求解器快速计算的快速计算的需要需要 薄壁梁碰撞性能优化（薄壁梁碰撞性能优化（Jiang and Han, Numerical Methods in Engineering, 2008)基于近似模型的高效优化求解基于近似模型的高效优化求解

3、(Jiang and Han, 2008, CMES Journal)Outer layer IP-GAIndividual of design vectorInner layer IP-GAInterval of the approximate objective functionInner layer IP-GAIntervals of the approximate constraintsPossibility degree of the approximate objective functionPossibility degrees of the approximate constr

4、aintsApproximate penalty functionCurrent design spaceUncertainty spaceStopping criterionApproximation models of the uncertain objective function and constraintsCallCallOptimal design vectorYNUncertainty space and current design spaceSampling points from LHDActual simulation models作为高效近似模型滤波，降噪滤波，降噪带

5、噪声带噪声通过插值方通过插值方法降噪后法降噪后 插值法插值法是解决此类问题的一种比较古老的、然而却是目前常用的方法，它不仅直接广泛地应用于生产实际和科学研究中，而且也是进一步学习数值计算方法的基础。9定义定义 设函数 y = f(x) 在区间 a,b上连续，且在n+1个不同的点 上分别取值 ，在一个性质优良、便于计算的函数类 中，求一简单函数p(x) ，使 而在其它点 上，作为 f(x) 的近似。称区间为插值区间插值区间，点 为插插值节点值节点，称（1.1）为 f(x)的插值条件插值条件，称函数类 为插值函数类插值函数类，称 p(x)为函数在bxxxan,10nyyy,10 niyxPii,

6、1 , 0(1.1)ixx nxxx,1010节点 处的插值函数插值函数。求插值函数 p(x) 的方法称为插值法插值法。插值函数类的取法不同，所求得的插值函数p(x)逼近f(x)的效果就不同它的选择取决于使用上的需要。常用的有代数多项式、三角多项式和有理函数等。 当选用代数多项式作为插值函数时，相应的插值问题就称为多项式插值。多项式插值。 在多项式插值中，最常见、最基本的问题是：求一次数不超过n的代数多项式nxxx,1011 (1.2) nnnxaxaaxP10使其中 为实数。满足插值条件(1.3)的多项式(1.2)，称为函数f(x) 在节点 处的n次插值多项式次插值多项式。 n次插值多项式次

7、插值多项式 的几何意义的几何意义：过曲线y = f(x) 上的n+1个点 作一条n次代数曲线 ，作为曲线y = f(x) 的近似，如图图2-1。 niyxPiin, , 1 , 0naaa,10 xPn)(xPyn), 1 , 0)(,(niyxii (1.3)12 xPyn xfy 0 x1xnxXab0y1yny0Y13 1 .2 插值多项式存在唯一性插值多项式存在唯一性 由插值条件（1.3）知，插值多项式 的系数满足线性方程组 （1.4）由线性代数知，线性方程组的系数行列式（记为V）是n+1阶范德蒙（Vandermonde）行列式，且 xPnniai, 1 , 0nnnnnnnnnyxa

8、xaayxaxaayxaxaa101111000010niijjinnnnnnxxxxxxxxxxxV11021211020011114 因 是区间 上的不同点，上式右端乘积中的每一个因子 ，于是 ，方程组（1.4）的解存在且唯一。故有下面的结论：定理定理1 若节点 互不相同，则满足插值条件（1.3）的次插值多项式（1.2）存在且唯一。 0Vnxxx, 1, 0ba,0jixxnxxx, 1,0152 拉格朗日插值多项式拉格朗日插值多项式 在上一节里，我们不仅指出了插值多项式的存在唯一性，而且也提供了它的一种求法，即通过解线性方程组（1.4）来确定其系数 ，但是，这种作法的计算工作量大，不便于

9、实际应用，下面介绍几种简便的求法。 2.1 插值基函数插值基函数 先考虑一下简单的插值问题：对节点 中任一点 ,作一n次多项式 , 使它在该点上取值为1,而在其余点 上取值为零, 即 （2.1）（2.1）表明n个点 都是n次多项式 的零点，故可设)(xlkianixi, 1 , 0nkxk0nkkixi, 1, 1, 1 , 0kikixlik01)()(xlknkkixi, 1, 1, 1 , 0)()()()(1110nkkkkxxxxxxxxxxAxl16其中 为待定系数，由条件 可得 故 (2.2)对应于每一节点 ，都能求出一个满足插值条件（2.1）的n次插值多项式（2.2），这样，由

10、（2.2）式可以求出n+1个n次插插多项式 。容易看出，这组多项式仅与节点的取法有关，称它们为在n+1个节点上的n次基本插值多项式次基本插值多项式或n次插值次插值基函数。基函数。kA1)(xlk)(,),(),(10 xlxlxln)()()(1110nkkkkkkkxxxxxxxxA)()()()()()()(110110nkkkkkknkkkxxxxxxxxxxxxxxxxxlnkxk017 2.2 拉格朗日插值多项式拉格朗日插值多项式 利用插值基函数立即可以写出满足插值条件（1.3）的n次插值多项式 （2.3） 事实上，由于每个插值基函数 都是n次多项式，故其线性组合（2.3）必是不高于

11、n次的多项式，同时，根据条件（2.1）容易验证多项式（2.3）在节点 处的值为 ，因此，它就是待求的n次插值多项式 。 形如(2.3)的插值多项式称为拉格朗日插值多项式，记为 （2.4)ix xLn)()()(1100 xlyxlyxlynn), 1 , 0)(nkxlkniyi, 1 , 0 xPn)()()(1100 xlyxlyxlynnnknkkkkkknkkkxxxxxxxxxxxxxxxxy0110110)()()()()()(18 作为的特例，令n=1，由（2.4）即得两点插值公式两点插值公式即这是一个线性函数，用线性函数 近似代替函数 ，在几何上就是通过曲线 上两点 作一直线

12、近似代替曲线 (见图图2-2)，故两点插值又名线性插值线性插值。 若令n=2，由（2.4）又可得常用的三点插值公式 （2.5） （2.6）（2.7）)()()()()()()(1202102210120120102102xxxxxxxxyxxxxxxxxyxxxxxxxxyxL010110101)(xxxxyxxxxyxL)()(0010101xxxxyyyxL)(1xL),(00yx),(11yx)(1xyL)(xfy )(xfy )(xf19这是一个二次函数，用二次函数 近似代替函数 ，在几何上就是通过曲线 上的三点 ，作一抛物线 近似地代替曲线 （图图2-3），故三点插值三点插值(二次插

13、二次插值值)。例例1 已知 分别用线性插值和抛物插值 求 的值。 xLy1 xfy x0 x0 x1),00(yx),11(yxy图图2-2)(2xL)(xf)(xfy )(xfy ),(),(),(221100yxyxyx)(2xLy 12144,11121,1010011520解解 因为115在100和121之间，故取节点x0=100，x1=121相应地有 y0=10，y1=11，于是，由线性插值公式（2.5）可得 故用线性插值求得的近似值为图图2-3 xLy2 xfy ),11(yx),00(yxyxx0 x10),22(yx714.10100121100115*111211001211

14、15*10)115(1151L100121100*11121100121*10)(1xxxL 21723.10)121144)(100144()121115)(100115(*12)144121)(100121()144115)(100115(*11)144100)(121100()144115)(121115(*10)115(1152 L115仿上，用抛物插值公式（2.7）所求得的近似值为将所得结果与 的精确值10.7238相比较，可以看出抛物插值的精确度较好。 为了便于上机计算，我们常将拉格朗日插值多项式（2.4）改写成公式（2.8）的对称形式编程框图如图图2-4，可用二重循环来完成 值的

15、计算，先通过内循环，即先固定k，令j从0到 ，累乘求得 nknkjjjkjknxxxxyxL00)(（2.8）)(kjn)(xLn22图图2-4输入xi,yi,n,xy=0k=0,1,nP=1j=0,1,nk=j?P=P*(x-xj)/(xk-xj)y=y+P*yk输出x,y否是23 然后再通过外循环，即令k从0到n，累加得出插值结果 。 2.3 插值余项插值余项 在插值区间a,b上用插值多项式 近似代替 ，除了在插值节点xi上没有误差外，在其它点上一般是存在有误差的。若记则 就是用 近似代替 时所产生的截断误差，称 为插值多项式 的余项余项。 关于误差有如下定理2中的估计式。定理定理2 设

16、在区间 上有直到n+1阶导数， 为区间 上n+1个互异的节点， 为满足条件:nkjjjkjkxxxxxl0)()()()(xPxfxRnn)(xLn)(xPn)(xf)(xRn)(xPn)(xPn)(xfba,nxxx,10ba,)(xPn)(xf24 的n次插值多项式，则对于任何 ，有其中 且依赖于 。 证明证明 由插值条件 知 ，即插值节点都是 的零点，故可设 (2.10)其中 为待定函数。下面求 ，对区间 上异于 的任意一点 作辅助函数 不难看出 具有如下特点: (2.11)（2.9）b)(a, )()(01niinxxx)()!1()()(1)1(xnfxRnnn), 1 , 0)()

17、(nixfxPiinbax,x)()(iinxfxP), 1 , 0(0)(nixRin)(xRn)()()(1xxKxRnn)(xK)(xKba,ixixx )()()()()(1txKtPtftFnn)(tF), 1 , 0(0)()(nixFxFi25例例2 在例1中分别用线性插值和抛物插值计算了的 近似值，试估计它们的截断误差。解解 用线性插值求 的近似值，其截断误差由插值余项公 式（2.9）知现在x0=100，x1=121，x=115，故115xxf)(xxxxxxxfxR10102/321,)(81)()( 21)(01125. 010*6*15*81max)121115)(100

18、115(*81)115(32/3121,1001R26 当用抛物插值求 的近似值时，其截断误差为 将 代入，即得 2.4 插插值误差的事后估计法值误差的事后估计法 在许多情况下，要直接应用余项公式（2.9）来估计误差是很困难的，下面将以线性插值为例，介绍另一种估计误差的方法。 设 若将用 两点作线性插值求得 的近似值记为 ，用 两点作线性插值所求得 的近似值记为 ，则由余项公式（2.9）知.)2 , 1 , 0)(210为已知且ixfxxxxi0017. 010*)144115)(121115)(100115(161)115(52RxxxxxRxxxxfx,202102/532),)()(16

19、1)()( ! 31)(xxf)(115,144,121,100210 xxxx10,xx)(xfy 1y20,xx2y)(xfy 27（2.13） 假设 在区间 中变化不大，将上面两式相除，即得近似式 即近似式（2.13）表明，可以通过两个结果的偏差 来估计插值误差 ，这种直接利用计算结果来估计误差的方法，称为事后估事后估计法计法。例例3 在例1中，用 做节点，算得的 近似值为 ，同样，用 做节点，可算得 的另一近似值 ，（2.13）可以估计出插值结果 的误差为:,),)()( 21,),)()( 2120220221011011xxxxxxfyyxxxxxxfyy)(121211yyxxx

20、xyy115115)( xf20,xx2121xxxxyyyy12yy 1yy 121,10010 xx144,10020 xx714.101y682.102y1y28 00835. 0)714.10682.10(1211441211151151 y 3 牛顿插值多项式牛顿插值多项式 由线性代数可知，任何一个不高于n次的多项式，都可表示成函数 的线性组合，即可将满足插值条件 的n次多项式写成形式其中 为待定系数。这种形式的插值多项式称为牛顿牛顿Newton插值多项式插值多项式，我们把它记成nx,即 （3.1)()()()(110102010nnxxxxxxaxxxxaxxaa 11012010

21、nnonxxxxxxaxxxxaxxaaxNnkak, 1 , 0)()(,),)(, 1110100nxxxxxxxxxxxx), 1 , 0()(niyxPii29 因此，牛顿插值多项式 是插值多项式 的另一种表示形式,与拉格朗日插值多项式相比较，不仅克服了“增加一个节点时整个计算机工作必须重新开始”见例1的缺点，而且可以节省乘除法运算次数。同时，在牛顿插值多项式中用到的差分与差商等概念，又与数值计算的其它方面有着密切的关系. 3.1 向前差分与牛顿插值公式向前差分与牛顿插值公式 设函数x 在等距节点 处的函数值 为已知，其中h是正常数，称为步长步长，称两个相邻点 和 处函数值之差 为函数

22、x在点 处以h为步长的一阶向前差分一阶向前差分简称一阶差分，记作 ，即于是，函数x 在各节点处的一阶差分依次为 又称一阶差分的差分为二阶差分二阶差分。 xNn xPnnkkhxxk, 1 , 00kkyxf1kxkkyy1kxkykkkyyy1。11121010,nnnyyyyyyyyykkkkyyyy12kx30kmkmkmyyy111一般地，定义函数 x在点 处的m阶差分阶差分为 为了便于计算与应用，通常采用表格形式计算差分，如表表2-1所示。表表2-1xkykykyk2yk3yk4x0 x2x3x4x1y0y1y2y3y4y0y1y2y3y02y12y22y03y13y04kx31 在等

23、距节点 情况下，可以利用差分表示牛顿插值多项式3.1 的系数，并将所得公式加以简化。事实上，由插值条件 立即可得 再由插值条件 可得由插值条件 可得 一般地，由插值条件 可得 2020121202020022! 222hyhhyyyxxxxhxxyyya), 1 , 0(0nkkhxxk00yxNn00ya 11yxNnhyxxyya001011kknyxN22yxNn32 于是，满足插值条件 的插值多项式为令 ,并注意到 ，则可简化为 这个用向前差分表示的插值多项式，称为牛顿向前插值公式牛顿向前插值公式，简称前插公式前插公式。它适用于计算表头 附近的函数值。 由插值余项公式2.9，可得前插公

24、式的余项为：), 2 , 1(!nkhkyakokk iinyxN 110010202000!2nnnnxxxxxxhnyxxxxhyxxhyyxN)0(0tthxxkhxxk0002000!11!21ynntttyttytythxNnn0 x3233 （3.3)例例4 从给定的正弦函数表表表2-2左边两列出发计算 ,并估计截断误差。 ),(,!110110nnnnxxfhnntttthxR)12. 0sin(表表22xxsinyy2y30.10.20.30.40.50.60.295520.198670.099830.479430.389420.564640.098840.096850.093

25、900.090010.08521-0.00389-0.00295-0.00094-0.00096-0.00480-0.00091-0.0019934解解 因为0.12介于0.1与0.2之间，故取 ，此时 。为求 ,构造差分表表22。表中长方形框中各数依次为 在 处的函数值和各阶差分。若用线性插值求sin0.12的近似值，则由前插公式3.2立即可得用二次插值得用三次插值得：1 . 00 x2 . 01 . 01 . 012. 00hxxtyyy03020,1 . 00 x11960. 009884. 02 . 009983. 0)12. 0()12. 0sin(1 N11976. 000016.

26、 0)12. 0()00199. 0(2) 12 . 0(2 . 009884. 02 . 009983. 0)12. 0()12. 0sin(12NNxsin35 因 很接近，且由差分表表22可以看出，三阶差分接近于常数（即 接近于零），故取 作为 的近似值，此时由余项公式（3.3）可知其截断误差 3.2 向后差分与牛顿向后插值公式向后差分与牛顿向后插值公式 在等距节点 下， 除了向前差分外，还可引入向后差分和中心差分向后差分和中心差分，其定义和记号分别如下： 在点 处以h为步长的一阶向后差分一阶向后差分和m阶向后差阶向后差分分分别为)12. 0sin( 11971.000096.0622

27、.012 .02 .0)12.0()12.0()12.0sin(23NN)12. 0()12. 0(23NN与04y11971. 0)12. 0(3N xfy kx000002. 0)4 . 0sin() 1 . 0(24)32 . 0()22 . 0() 12 . 0(2 . 0)12. 0(43R), 1 , 0(0nkkhxxk36 在 点处以h为步长的一阶中心差分和m阶中心差分分别为其中 各阶向后差分与中心差分的计算，可通过构造向后差分表与中心差分表来完成?参见表表2?。 利用向后差分，可简化牛顿插值多项式（.1），导出与牛顿前插公式?3.2?类似的公式，即，若将节点的排列次序看作 ，

28、那么?.1)可写成 xfy kx), 3 , 2(2112112121myyyyyykmkmkmkkk.2,22121hxfyhxfykkkk01,xxxnn), 3 , 2(1111myyyyyykmkmkmkkk37根据插值条件 , 可得到一个用向后差分表示的插值多项式其中t0，插值多项式(3.4)称为牛顿向后插值公式牛顿向后插值公式，简称后插公式。它适用于计算表尾 附近的函数值。由插值余项公式（.9），可写出后插公式的余项(3.4) 111210 xxxxxxbxxxxbxxbbxNnnnnnnn )0 , 1 , 1,(nniyxNiinnnnnnnnynntttyttytythxN!

29、11! 212nx38 (3.5)例例已知函数表同例，计算 ，并估算截断误差。解解因为.58位于表尾 附近，故用后插公式（3.4）计算sin(0.58)的近似值。 一般地为了计算函数在 处的各阶向后差分，应构造向后差分表。但由向前差分与向后差分的定义可以看出，对同一函数表来说，构造出来的向后差分表与向前差分表在数据上完全相同。因此，表表-用“”线标出的各数依次给出了 在 处的函数值和向后差分值。因三阶向后差分接近于常数，故用三次插值进行计算，且 ，于是由后插公式（3.4）得 ),(!11011nnnnnxxfhnntttthxR)58. 0sin(6 . 05x5xxsin6 . 05x2 .

30、 01 . 0/ )6 . 058. 0(/ )(5hxxt39 因为在整个计算中，只用到四个点 上的函数值，故由余项公式（.5）知其截断误差 54802. 000091. 0622 . 012 . 02 . 000480. 0212 . 02 . 03 . 0 , 4 . 0 , 5 . 06 . 0 ，x58. 058. 0sin3N08521. 02 . 056464. 0000002. 06 . 0sin1 . 04 2432 . 022 . 012 . 02 . 058. 03R40 3.3 差商与牛顿基本插值多项式差商与牛顿基本插值多项式 当插值节点非等距分布时，就不能引入差分来简

31、化牛顿插值多项式，此时可用差商这个新概念来解决。 设函数 在一串互异的点 上的值依次为 。我们称函数值之差 与自变量之差 的比值为函数 关于 点的一阶差商一阶差商，记作例如210iiixxx、)()()(210iiixfxfxf、0101)()(iiiixxxfxf,10iixxf,)()(,)()(,121221010110 xxxfxffxxxfxffxxxx01iixx )(xf)()(01iixfxf)(xf01,iixx41称一阶差商的差商 为函数 关于点 的二阶差商二阶差商（简称二阶差商二阶差商），记作 ，例如 一般地，可通过函数 的m-1阶差商定义的m阶差商如下：021021,i

32、iiiiixxxxfxxf xf210iiixxx、021021210,xxxxfxxfxxxf,210ixxxfii010110,iiiiiiiiixxxxfxxfxxxfmmmm)(xf42 差商计算也可采用表格形式（称为差商表），如表表23所示， 表表23 1xf 0 xf32,xxf 一阶差商 二阶差商 三阶差商 kxf10,xxf210,xxxf21,xxf3210,xxxxf2xf321,xxxf3xfkx3x0 x2x1x43差商具有下列重要性质（证明略）：（1） 函数 的m阶差商 可由函数值 的线性组合表示，且（2） 差商具有对称性，即任意调换节点的次序，不影响差商的值。 例如

33、（）当 在包含节点 的某个区间上存在时， 在 之间必有一点 使,10mxxxf mxfxfxf、10 mimiiiiiiimxxxxxxxxxfxxxf011010)()()(,.,021201210 xxxfxxxfxxxf xfmmjxji, 1 , 0miiixxx,10, !,10mfxxxfmiimi，)(xf44（）在等距节点 情况下，可同时引入 阶差分与差商，且有下面关系： 引入差商的概念后，可利用差商表示牛顿插值多项式（.1）的系数。事实上，从插值条件出发，可以象确定前插公式中的系数那样，逐步地确定（.1）中的系数故满足插值条件 的n次插值多项式为nkkhxxk, 1 , 00

34、nmmmnmmnnnmmmhmyxxxfhmyxxxf!,!,1010)(, 2 , 1,0010 xfankxxxfakk niyxNiin, 1 , 045 （3.6）（3.6）称为牛顿基本插值多项式牛顿基本插值多项式，常用来计算非等距节点上的函数常用来计算非等距节点上的函数值。值。例例 试用牛顿基本插值多项式按例要求重新计算 的近似值。解解 先构造差商表。 由上表可以看出牛顿基本插值多项式(3.6)中各系数依次为11010,nnxxxxxxxxxf 102100100,xxxxxxxfxxxxfxfxNn115一阶商差 二阶商差xx10012111100.0434780.047619-0

35、.0000941441246 故用线性插值所得的近似值为 用抛物插值所求得的近似值为所得结果与例1相一致。比较例1和例6的计算过程可以看出，与拉格朗日插值多项式相比较，牛顿插值多项式的优点是明显的。 由插值多项式的存在唯一性定理知，满足同一组插值条件的拉由插值多项式的存在唯一性定理知，满足同一组插值条件的拉格朗日插值多项式格朗日插值多项式（2.42.4）与牛顿基本插值多项式与牛顿基本插值多项式（3.63.6）是同一多是同一多项式。因此，余项公式项式。因此，余项公式（2.92.9）也适用于牛顿插值。但是在实际计也适用于牛顿插值。但是在实际计算中，有时也用差商表示的余项公式算中，有时也用差商表示的

36、余项公式000094. 0,047619. 0,10)(210100 xxxxxxfff7143.10)100115(047619. 010)115(1151 N7228.10)121115()100115()000094. 0()115()115(11512NN47 （3.73.7）来估计截断误差（证明略）。注意注意： 上式中的n+1阶商差 与 的值有关，故不 能准确地计算出 的精确值，只能对它作一 种估计。例，当四阶差商变化不大时，可用 近似代替 。)(xf)(,)(110 xxxxxfxRnnn,10mxxxf,10mxxxf,43210 xxxxxf,3210 xxxxxf 4 埃尔米

37、特插值埃尔米特插值.Hermite ,)插值问题)插值问题埃尔米特(埃尔米特( . .种插值问题称为这相等导数值甚至高阶导数值在某些节点上对应的的而且要求在节点上函数值相等不少插值问题不仅要求(5.1) ., 1 , 0 ,)( ,)( )(12 , 1 , 0,)()(112121210nifxHfxHxHnnifxffxfbxxxaniiniinniiiin满足次的多项式超过要求一个次数不和上已知个节点着重讨论一种情况：在次多项式，且满足是都和函数采用基函数法，插值基12), 1 , 0)()(nnjxxjj(5.2) ), 1 , 0,( .)( , 0)( ; 0)( ,)( nkjx

38、xxxjkkjkjkjjkkj (5.3) . )()()(012njjjjjnxfxfxH ),()()( 2xlbaxxjj令(5.4) ).(1)(21)( 20 xlxxxxxjnjkkkjjj ),()()( 2xlxxAxjjj令(5.5) ).()()( 2xlxxxjjj .唯一性，. ),(5.6) ),()!22()()()()( ),(22),()(2)22(12)22(1xbaxnfxHxfxRxfnbaxfnnnn且依赖于其中插值余项为则阶导数内有在插值区间若. )()()(1)( 21 )()()(00220012njnjjjjjjnjkkkjjnjjjjjnfxl

39、xxfxlxxxxxfxfxH尔米特插值多项式：时，应用广泛的三次埃当1n .)()( 2121)(12010102101012010101021010103fxxxxxxfxxxxxxfxxxxxxxxfxxxxxxxxxH ).,( ,)()(! 4)( )()()( 102120)4(33xxxxxxfxHxfxR余项为例例 给定 f ( 1)=0, f (1)=4, f ( 1)=2, f (1)=0, 求H3(x), 并计算 f (0.5).解解 x0 = 1, x1 = 1, f(0.5)H3(0.5) = 3.5625.3010110( )( ) 0( ) 4( ) 2( ) 0

40、4( )2( ).Hxh xh xgxg xh xgx 221220( 1)11( )12(2)(1) ,1 11 ( 1)411( )( 1)(1)(1) ,1 14xxh xx xxgxxxx 2231( )(2)(1)(1)(1) ,2Hxx xxx54 4 分段低次插值分段低次插值 例2、例4表明，适当地提高插值多项式的次数，有可能提高计算结果的准确程度。但是决不可由此提出结论，认为插值多项式的次数越高越好。例如，对函数 先以 为节点作五次插值多项式 ，再以 为节点作十次插值多项式 ，并将曲线 描绘在同一坐标系中，如图图2-52-5所示。)10, 1 ,0(511iixi)11(251

41、1)(2xxxf)5, 1 ,0(521iixi)(5xP)(10 xP)1 , 1)(),(,2511)(1052xxPyxPyxxf55 -1 0 1 x y 1y=1/(1+25x2)y=P5(x)图图2-52-5y=P10(x)56 由上图可看出，虽然在局部范围中，例如在区间-0.2 ，0.2中， 比 较好地逼近 ,但从整体上看， 并非处处都比 较好地逼近 ，尤其是在区间-1 ，1的端点附近。进一步的分析表明，当n增大时，该函数在等距接点下的高次插值多项式 ，在-1 ，1两端会发生激烈的振荡。这种现象称为龙格龙格（Runge)现象现象。这表明，在大范围内使用高次插值，逼近的效果可在大范

42、围内使用高次插值，逼近的效果可能不理想的。能不理想的。 另一方面，插值误差除来自截断误差外，还来自初始数据的误差和计算过程中的舍入误差。插值次数越高，计算工作越大，积累误差也可能越大。 因此，在实际计算中，常用分段低次插值进行计算 ，即把整个插值区间分成若干小区间，在每个小区间上进行低次插值。 例如例如，当给定n+1个点 上的函数值 后，若要计算点 处函数值 的近似值，可先选取两个节点 使 然后在小区间 上作线性插值，即得 )(5xP)(10 xP)(xf)(10 xP)(5xP)(xPnnxxx10nyyy10ixx )(xf)(xfiixx,1,1iixxx,1iixx57 （4.1）这种

43、分段低次插值叫分段线性插值分段线性插值。在几何上就是用折线代替曲线，如图图2-6所示。故分段线性插值又称折线插值折线插值.11111)()( iiiiiiiixxxxyxxxxyxPxf xy=f(x)yox0 x1xnx2图2-658 类似地，为求 的近似值，也可选取距点 最近的三个节点 进行二次插值，即取这种分段低次插值叫分段二次插值分段二次插值。在几何上就是用分段抛物线代替曲线，故分段二次插值又称分段抛物插值分段抛物插值。为了保证 是距点 最近的三个节点，（4.2）中的 可通过下面方法确定:)(xfyx11,iiixxx)()()(11112ikjijjkjiikkxxxxyxPxf（4

44、.2）11,iiixxxxixxxxnxxxxxjxxxxinnnjjjj1211210211212121159I n-1输入xi ,yi (I=0 ,1, ,n)及n,xj=1,2,n-2P 0.5*(xj +xj+1)xP?I j按公式（4.2）计算P2(x)输出x, P2(x)图图2-7605 三次样条插值三次样条插值 分段低次插值虽然具有计算简单、稳定性好、收敛性有保证且易在电子计算机上实现等优点，但它只能保证各小段曲线在连接点上的连续性，却不能保证整条曲线的光滑性(如图图2-6中的折线)，这就不能满足某些工程技术上的要求。从六十年代开始，首先由于航空、造船等工程设计的需要而发展起来的

45、61所谓样条（Spline)的插值方法，既保留了分既保留了分段低次插值多项式的各种优点，又提高了插段低次插值多项式的各种优点，又提高了插值函数的光滑性。值函数的光滑性。今天，样条插值方法已成为数值逼近的一个极其重要的分支，在许多领域里得到越来越广泛的应用。 本节介绍应用最广泛且具有二阶连续导数的三次样条插值函数。 5.1 三次样条插值函数的定义三次样条插值函数的定义 对于给定的函数表62x)(xfy0 x0 x1xny1yn其中 ，若 函数 满足: （1） 在每个子区间 上是不高于三次的多项式； （2） 在a,b上连续； （3）满足插值条件 bxxxan10)(xS), 2 , 1(,1nix

46、xii)(xS)( ),( ),(xSxSxS), 1 , 0()(niyxSii63i)(xS)(xfnxxx,10)(xSi)(xS,1iixx)(xSi则称 为函数 关于节点 的三三次样条插值次样条插值。5.2 边界条件问题的提出与类型边界条件问题的提出与类型 注：注：单靠一张函数表是不能完全确定一个三次样条插值函数的。 事实上，由条件（1）知，三次样条插值函数 是一个分段三次多项式， 若用 表示它在第 个子区间 上的表达式，则 形如：,)(1332210iiiiiiixxxxaxaxaaxS64这里有四个待定系数 。子区间共有n n个，确定 需要确定4n n个待定系数。 另一方面，要求

47、分段三次多项式 及其导数 在整个插值区间a,b上连续，只要在各子区间的端点 连续即可。故由条件（2），（3）可得待定系数应满足的4n n-2个方程为： ) 3 , 2 , 1 , 0( jaij)(xS)(xS)( ),( xSxS) 1, 2 , 1(nixi65),.,1 , 0()() 1,.,2 , 1()0( )0( ) 1,.,2 , 1()0( )0( ) 1,.,2 , 1()0()0(niyxSnixSxSnixSxSnixSxSiiiiiiii(5.1) 由此可以看出，要确定4n个待定系数还缺少两个条件，这两个条件通常在插值区间a,b的边界点a,b处给出，称为边界条件边界条

48、件。边界条件的类型很多，常见的有：（1）给定一阶导数值（2）给定二阶导数值00)( ,)( nnyxSyxS00)( ,)( nnyxSyxS66 特别地， 称为自然边界自然边界条件条件，满足自然边界条件自然边界条件的三次样条插值函数称为自然样条插值函数自然样条插值函数。（3）当 是周期为 的函数时，则要求 及其导数都是以 为周期的函数，相应的边界条件为0)( )( 0nxSxS)(xfab)(xSab)0( )0( ),0( )0( 00nnxSxSxSxS675.3 三次样条插值函数的求法三次样条插值函数的求法 虽然可以利用方程组（5.1）和边界条件求出所有待定系数 aij 从而得到三次样

49、条插值函数 在各个子区间 的表达式 。但是，这种做法的计算工作量大，不便于实际应用。下面介绍一种简便的方法。 设在节点 处 的二阶导数为ix)(xS)(xS,1iixx)(xSi), 1 , 0()( niMxSii68 因为在子区间 上 是不高于三次的多项式，其二阶导数必是线性函数（或常数）。于是，有记 则有,)(11111 iiiiiiiiiiixxxxxxxMxxxxMxS,1iixx)()(xSxSi1iiixxhiiiiiiihxxMhxxMxS11 )(69连续积分两次得：BAhxxMhxxMxSiiiiiiixxiii1331616)(5.2)其中 为积分常数。利用插值条件 ii

50、BA,iiiiiiyxSyxS)(,)(11易得70 将它们代入（5.2） ，整理得112111()616iiiiiiiiiiyyAhByhMMM ), 2 , 1,()6()6(616)(112211331nixxxhxxiyhxxhyhihixiiiiiiiiiiiiiiiihMMxxMxxMS(5.3)71 综合以上讨论可知，只要确定这n+1个值，就可定出三次样条插值函数。 为了求出 ，利用一阶函数在子区间连接点上连续的条件 ), 1 , 0(niMi), 1 , 0(niMi)0()0( )0( )0( 1iiiiiixSxSxSxS即（5.4）由（5.3）可得72 116iiiiii

51、MMnhyy(5.5) iiiiiihxxMhxxMxS222121故iiiiiiiiiMhMhhyyxS36011(5.6) 将（5.5）中的 改为 ，即得 在子区间 上的表达式 ，并由此得： i1i)( xS,1iixx)(1xSi73将(5.6),(5.7)代入(5.4)整理后得111111630iiiiiiiiiMhMhhyyxS(5.7)iiiiiiiiiiiiihyyhyyMhMhhMh1111111636两边同乘以 ，即得方程组16iihh1, 2 , 162111111111nihyyhyyhhMhhhMMhhhiiiiiiiiiiiiiiiii74若记1, 2 , 1,61,

52、111111nixxfxxfhhghhhhhhiiiiiiiiiiiiiiii(5.8)则所得方程组可简写成1, 2 , 1211nigMMMiiiiii75即1112 1232212121101222nnnnnngMMMgMMMgMMM(5.9)这是一个含有 n+1个未知数、n-1 个方程的线性方程组。要确定 的值，还需用到边界条件。在第 (1) 种边界条件下，由于iMoyxS0nnyxS和76已知，可以得到包含 另外两个线性方程。由(5.5)知， 在子区间 上的导数为 xS10,xxiM 011101120112101622MMhhyyhxxMhxxMxS故由条件 立即可得oyxS0011

53、10110062MMhhyyhMy77即010111062yhyyhMM(5.10)同理，由条件 可得nnyxSnnnnnnnhyyyhMM1162（5.11） 将(5.9)、(5.10)、(5.11)合在一起，即得确定 的线性方程组：nMMM,1078nnnnnnggggMMMM1101101111212212(5.12)其中nnnnnxxfyhgyxxfhg,6,6101010(5.13)79已知，在方程组(5.13)中实际上只包含有n-1个未知数 ，并且可以改写成 000yxSM nnnyxSM121,nMMM在第(2) 种边界条件下，由 nnnnnnnnnygggygMMMM11220

54、1112211222212222(5.14)80在第第(3)种边界种边界条件下，000nxSxS(5.15)nMM0由条件 可得000 nxSxS1101 1101106262nnnnnnnnMMhhyyhMMMhhyyhM由81注意到 和 ，上式整理后得 nyy 0nMM0nnnnnnnnnhyyhyyhhMMhhhMhhh110111111162若记nnnnnnnnnnnxxfxxfhhghhhhhh,61,10111182则所得方程可简写成nnnnngMMM211(5.16)将 (5.9) 、(5.15) 、(5.16) 合在一起，即得确定 的线形方程组：nMMM,21nnnnnnnng

55、gggMMMM121121112212222(5.17)83利用线性代数知识，可以证明方程组 (5.12) 、(5.14) 及(5.17)的系数矩阵都是非奇异的，从而都有唯一确定的解。 针对不同的边界条件，解相应的方程组(5.12) 、 (5.14)或(5.17) ，求出的值，将它们代入（5.3），就可以得到 在各子区间上的表达式。综上分析，有 nMMM,20 xS84xnx xfy ny1y0y1x0 xbxxxan10定理定理3 对于给定的函数表满足第（1）或第（2）或第（3）种边界条件的三次样条插值函数是存在且唯一的存在且唯一的。 三次样条插值函数 的具体求解过程，在下面例子中给出了详细

56、说明。 xS85例例 7 已知函数 的函数值如下 xfy yx-1.5 0 1 20.125 -1 1 9在区间 上求三次样条插值函数 ，使它满足边界条件：2, 5 . 1 142,75. 05 . 1SS xS解解 先根据给定数据和边界条件算出 giii与,iM，写出确定 的线性方程组。86在本例中，给出的是第（1）种边界条件，确定 的线性方程组形如（5.12） 。由所给函数表知 3 , 2 , 1 , 0iMi5 . 11h12h13h75. 0,10 xxf2,21xxf8,32xxf于是由 的算式(5.8) 知1, 2 , 1,nigiii与186 . 65 . 04 . 05 . 0

57、6 . 0212121gg8736,66,63233301010 xxfyhgyxxfhg 由第（1）边界条件下 与 的计算公式（5.13）知0gng故确定 的方程组为3210,MMMM与36186 . 6621005 . 025 . 0004 . 026 . 000123210MMMM（5.18）88然后解所得方程组，得到 在各节点 上的值 。在本例中，解（5.18）得 xS ixiM16, 4, 4, 53210MMMM最后将所得 代入（5.3），即得 在各子区间上的表达式 。由（5.3）知， 在 上的表达式为 xSiM10,xx nixSi, 2 , 1 xS 1021111121001

58、301131016666hxxhMyhxxhMyhxxMhxxMxS89在本例中，将4, 5, 1,125. 0, 0, 5 . 1101010MMyyxx代入，整理后得 12231xxxS0, 5 . 1x同理可得 3642233xxxxS1 , 0 x 1222 xxS2, 1x90故所求三次样条插值函数为 3642121223223xxxxxxxS211005 . 1xxx第（1）边界条件下计算三次样条插值函数S（x）在 x 处函数值 的程序框图如图2-891xnyyniyxnii, , ), 1 , 0(,0及输入), 2 , 1(,1nixxfhiii与计算) 1., 2 , 1(,

59、)8 . 5(nigiii计算按公式), 1 , 0()12. 5(niMi得解方程组值式计算并按式中的确定根据)() 3 . 5() 3 . 5(xsix)(,xsx输出图图2-8ggn与计算按公式0)13. 5(92 上述求三次样条插值函数的方法，其基本思路和特点是: 先利用一阶导数 在内节点 上的连续性以及边界条件，列出确定二阶导数 的线性方程组（在力学上称为三弯矩方程组三弯矩方程组），并由此解出 ，然后用 来表达 。)(xS) 1., 2 , 1(nixi nixSMii, 1 , 0 xSiMiM93 通过别的途径也可求三次样条插值函数。例如，可以先利用二阶导数在内节点上的连续性以及

60、边界条件，列出确定一阶导数 的线性方程组（在力学上称为三转角方程组三转角方程组），并由此解出 ，然后用 来表达 。在有些情况下，这种表达方法与前者相比较，使用起来更方便1。 nixSmii, 1 , 0imim ixS946 数值微分（敏感性分析）数值微分（敏感性分析） 作为多项式插值的应用，本节介绍两种求函数导数的近似值的方法。 6.1 利用插值多项式求导数利用插值多项式求导数 若函数 在节点 处的函数值已知，就可作 的n次插值多项式 ，并用 近似代替 ，即 xfnixi, 1 , 0 xf)(xPn xf)(xPn)()(xPxfn(6.1)95由于 是多项式，容易求导数，故对应于 的每一