1、12.5 2.5 埃尔米特插值埃尔米特插值 有些实际的插值问题不但要求在节点上函数值相等, 下面只讨论函数值与导数值个数相等的情况. 满足这种要求的插值多项式就是埃尔米特插值多项式埃尔米特插值多项式. 而且还要求对应的导数值也相等,甚至要求高阶导数也相等.2), 1 , 0(,)(,)(njmxHyxHjjjj(5.1)这里共有 个插值条件,可唯一确定一个次数不超过22n的多项式 ,12n)()(12xHxHn), 1 ,0()(njxfmjj问题是求插值多项式 ,)(xH 设在节点 上,),(jjxfybxxxan10.)(12121012nnnxaxaaxH 现在仍采用求拉格朗日插值多项式
2、的基函数方法. 满足条件 其形式为3, 1,0)(kjkjxjkkj将满足条件(5.1)的插值多项式 写成用插值基函数表示的形式 )()(12xHxHn.)()()(012njjjjjnxmxyxH(5.3) 先求出 个插值基函数 及 ,)(xj), 1 , 0()(njxj22n每一个基函数都是 次多项式,12n且满足条件;0)(kjx,0)(kjxjkkjx)(5.2)., 1 ,0,(nkj4令 ),()()(2xlbaxxjj由条件(5.2),有 , 1)()()(2jjjjjxlbaxx, 0)()(2)()()(jjjjjjjjjxlbaxxalxlx 由插值基函数所满足的条件(5
3、.2),有 ), 1 , 0(,)(,)(1212nkmxHyxHkknkkn下面的问题就是如何求出这些基函数 及 )(xj),(xj利用拉格朗日插值基函数).( xlj5解出 ).(21),(2jjjjjxlxbxla由于 ,)()()()()()()(110110njjjjjjnjjjxxxxxxxxxxxxxxxxxl整理得 .0)(2;1jjjxlabax6于是 ).()1)(21 ()(20 xlxxxxxjnjkkkjjj(5.4)两端取对数再求导,得 ,1)(0njkkkjjjxxxl同理,可得 ).()()(2xlxxxjjj(5.5)7 可以证明满足条件(5.1)的插值多项式
4、是惟一的. 用反证法,假设 及 均满足条件(5.1),)(12xHn)(12xHn)()()(1212xHxHxnn这样, 有 重根,但 是不高于 次的多 项式,)(x22n)(x12n于是在每个节点 上的值及导数值均为零,即 为二重根.kxkx.0)(x故惟一性成立.8其中 且与 有关. ),(bax 若 在 内的 阶导数存在,则其插值余项)(xf),(ba22n)()!22()()()()(21)22(12xnfxHxfxRnnn(5.6)仿照拉格朗日插值余项的证明方法,可以证明:9 插值多项式(5.3)的重要特例是 的情形. 1n这时可取节点为 及 ,kx1kx插值多项式为 ,)(3xH
5、,)(3kkyxH,)(3kkmxH.)(;)(113113kkkkmxHyxH(5.7)相应的插值基函数为),(),(),(),(11xxxxkkkk它们满足条件 满足10,0)(,1)(1kkkkxx;0)()(111kkkkxx,1)(,0)(111kkkkxx,0)()(1kkkkxx,0)()(111kkkkxx,0)()(1kkkkxx,0)(,1)(1kkkkxx.1)(,0)(111kkkkxx11.21)(,21)(211112111kkkkkkkkkkkkkkxxxxxxxxxxxxxxxxxx(5.8).)(,)(2111211kkkkkkkkkkxxxxxxxxxxxx
6、xx(5.9)根据 及 的一般表达式(5.4)及(5.5),)(xj)(xj可得到12),()()()()(11113xmxmxyxyxHkkkkkkkk(5.10)其余项 ,)()()(33xHxfxR).,(,)()(! 41)(1212)4(3kkkkxxxxxxfxR于是满足条件(5.7)的插值多项式是 由(5.6)得13 求满足 及 )2, 1 ,0()()(jxfxPjj)()(11xfxP 由给定的4个条件,可确定次数不超过3的插值多项式. 由于此多项式通过点),(,(),(,(),(,(221100 xfxxfxxfx)(,)()(0100 xxxxfxfxP)(,10210
7、xxxxxxxf的插值多项式及其余项表达式.例例4 4故其形式为),)()(210 xxxxxxA14.)(,)(,)(210121001101xxxxxxxfxxxxfxfA待定常数 ,可由条件 确定,)()(11xfxPA其中 为待定函数. )(xk为了求出余项 的表达式,)()()(xPxfxR),()()()(2210 xxxxxxxkxR通过计算可得 可设15).()()()()()(2210 xtxtxtxktPtft显然),2, 1 ,0(0)(jxj故 在 内有5个零点(二重根算两个). )(t),(ba, 0)(! 4)()()4()4(xkf 反复应用罗尔定理,得 在 内至
8、少有一个零点,)()4(t),(ba构造,0)(,0)(1xx且故有16),()()(! 41)(2210)4(xxxxxxfxR(5.11)式中 位于 和 所界定的范围内.210,xxxx余项表达式为 于是 ),(!41)()4(fxk172.6 2.6 分段低次插值分段低次插值 2.6.1 2.6.1 高次插值的病态性质高次插值的病态性质 这是因为对任意的插值节点,当 时, 不一定收敛到 .n)(xLn)(xf在次数 增加时逼近 的精度不一定也增加.n)(xf 根据区间 上给出的节点做出的插值多项式),(xLn,ba18), 1 , 0(,105nknkxk所构造的拉格朗日插值多项式为 以
9、 上的 个等距节点5 , 51n 考虑函数 ,它在 上的各阶导数均存在.5 , 5)1/(1)(2xxf.)()()(11)(1102jnjnnjjnxxxxxxL令),(2112/1nnnxxx,552/1nxn则19 表2-5列出了 时的 的计算结果及 20, 4 , 2n)(2/1nnxL在 上的误差2/1nx).(2/1nxR994889.39952449.39042440.020080751.20123671.20042920.018217397.10173867.10043530.016288409.5332743.5044334.014800440.2755000.2045440
10、.012531662.1578721.1047059.010880668.0831017.0049651.08553416.0607879.0054463.06423216.0356826.0066390.04621684.0759615.0137931.02)()()(2/12/12/1nnnnxRxLxfn5表220 可见,随 的增加, 的绝对值几乎成倍增加. n)(2/1nxR 这说明当 时 在 上是不收敛的. nnL5 , 5 Runge证明了,存在一个常数 ,使得当 63. 3ccx时, 而当 时 发散.cx )(xLn),()(limxfxLnn21 取 根据计算画出 及 ,10n
11、)1/(12xy)(10 xLy 在 上的图形,见图2-5.5 , 5图2-522 从图上看到,在 附近, 与 5x)(10 xL)1/(1)(2xxf偏离很远, 这说明用高次插值多项式 近似 效果并不好.)(xLn)(xf通常不用高次插值,而用分段低次插值. 23下图是用Matlab完成的Lagrange插值(附程序):24附:Lagrange插值程序n=11; m=61;x= -5:10/(m-1):5;y=1./(1+x.2);z=0*x;x0=-5:10/(n-1):5;y0=1./(1+x0.2);y1=lagr1(x0, y0, x);plot(x, z, r, x, y, k:
12、,x, y1, r)gtext(Lagr.), gtext(y=1/(1+x2)title(Lagrange)25附:Lagrange插值子程序 lagr1:function y=lagr1(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;end26 2.6.2 2.6.2 分段线性插值分段线性插值 所谓分段线性插值就是通过插值点用折线段连接
13、起来逼近 ).(xf 由于升高插值多项式的阶数有时并不能达到提高精度的效果, 所以实际中往往采用分段插值的思想. 分段插值的基本思想是将插值区间划分为若干个小区间, 然后在每个小区间上做满足一定条件的低阶插值.27 设已知节点 上的函数值 bxxxan10,10nfff记,max,1kkkkkhhxxh,)(.1baCxIh)(xIh求一折线函数 , 满足:kkhfxI)(. 2), 1 ,0(nk 在每个小区间 上是线性函数.)(. 3xIh,1kkxx则称 为分段线性插值函数分段线性插值函数.)(xIh28 由定义可知 在每个小区间 上可表示为 )(xIh,1kkxx1111)(kkkkk
14、kkkhfxxxxfxxxxxI).(1kkxxx(6.1) 若用插值基函数表示,则在整个区间 上 为 ,ba)(xIh, )()(0njjjhxlfxI(6.2)其中基函数 满足条件 )(xlj), 1 ,0,()(nkjxljkkj其形式是29);0(1略去jxxxjj);(1略去njxxxjj.,11jjxxxbax,0,)(1111jjjjjjjxxxxxxxxxl(6.3) 利用插值余项(2.17)得到分段线性插值的误差估计)(max2)()(max1211kkxxxhxxxxxxxMxIxfkkkk30或写成 ,8)()(max22hMxIxfhbxa(6.4)其中.)(max2x
15、fMbxa 31),()()(110 xlxlxlkknjj当 时, ,1kkxxx故 ).()()()(1xfxlxlxfkk另一方面,这时 ).()()(11xlfxlfxIkkkkh这种性质称为局部非零性质局部非零性质. 分段线性插值基函数 只在 附近不为零,在其他地方均为零,)( xljjx利用 的局部非零性质及 知,)( xljnjjxl0)(132 现在证明 , )()(limxfxIhn)()()()()()()(111xlfxlfxfxlxlxIxfkkkkkkh)()()(1kkkhxlxl这里 是函数 在区间 上的连续模,即对任意)(h)(xf,ba两点 ,只要 就有,ba
16、xx ,hxx 考虑11)()()()(kkkkfxfxlfxfxl)(kh).(h33),()()(hxfxf 称 为 在 上的连续模连续模, )(h)(xf,ba当 时,,)(baCxf.0)(lim0hh就有 由前式可知,当 时有 ,bax ).()()(maxhxIxfhbxa因此,只要 ,就有,)(baCxf)()(lim0 xfxIhh在 上一致成立,,ba故 在 上一致收敛到 .,ba)(xf)(xIh34下图是用Matlab完成的分段线性插值(附程序):35附:分段线性插值程序n=11; m=61;x=-5:10/(m-1):5;y=1./(1+x.2);z=0*x;x0=-5
17、:10/(n-1):5;y0=1./(1+x0.2);y1=interp1(x0, y0, x);plot(x, z, r, x, y, k:, x, y1, r)gtext(Piece. linear.), gtext(y=1/(1+x2)title(Piecewise Linear)注:interp1(x0,y0,x)为Matlab中现成的分段线性插值程序.36 2.6.3 2.6.3 分段三次埃尔米特插值分段三次埃尔米特插值 分段线性插值函数 的导数是间断的,若在节点)(xIh上除已知函数值 外还给出导数值 ), 1 ,0(nkxkkf), 1 , 0(nkmfkk插值函数 ,)(xIh
18、,)(.11baCxIh), 1 ,0()(,)(.2nkfxIfxIkkhkkh 在每个小区间 上是三次多项式.)(.3xIh,1kkxx这样就可构造一个导数连续的分段满足条件37 根据两点三次埃尔米特插值插值多项式(5.10), 在区间 上的表达式为)(xIh,1kkxx21121121)(kkkkkkkkkkhxxxxfxxxxxxxxxI.)(211121211111kkkkkkkkkkkkkkfxxxxxxxfxxxxxxfxxxx(6.5)38 若在整个区间 上定义一组分段三次插值基函数,ba及 ,)(xj), 1 ,0()(njxj则 可表示为)(xIh),()()(0 xfxf
19、xIjnjjjjh(6.6)其中 , 分别表示为 )(xj)(xj(6.7);,0),(,21),0(,21)(1121111211其他略去略去njxxxxxxxxxxxjxxxxxxxxxxxxjjjjjjjjjjjjjjjjj39;其他略去略去,0),(,),0(,)(12111211njxxxxxxxxxjxxxxxxxxxxjjjjjjjjjjjjj(6.8) 由于 , 的局部非零性质,)(xj)(xj当 时,,1kkxxx只有 不为零,)(),(),(),(11xxxxkkkk于是 可表示为)(xIh).()()()()(1111xfxfxfxfxIkkkkkkkkh).(1kkxx
20、x(6.9)40, 1)(0 xj(6.10) 再研究 的收敛性. )(xIh,274)(kkhx,274)(1kkhx(6.11). )()()(100njjjjnjjjxxfxf 此外,当 是分段三次多项式时, 的插值多项式 就是它本身. )(xf)(xf)(xIh 例如,当 时,1)(xf 由 及 的表达式,直接得估计式 )(xj)(xj所以就有由于,0f41当 时就得 ,1kkxxx, 1)()(1xxkk(6.12)由(6.9)-(6.12),当 时还有 ,1kkxxx11)()()()()()(kkkkhfxfxfxfxxIxf)27/4(1kkkffh11)()()()(kkkk
21、hfxhfx,)27/4(1kkkffh这里 且依赖于 .),(,1kkxxx42 因此对 成立,bax . )(max2735)()(maxxfhxIxfbxahbxa(6.13)这表明用 逼近 时,它的界只依赖 ,而与 无关. )(xIh)(xfhx因此,当 时,bax)()(lim0 xfxIhh一致成立. 设 则当 时, 在 上一致收敛于 .,1baCf 0h)(xIh,ba)(xf从而得到:定理定理3 3432.7 2.7 三次样条插值三次样条插值44 2.7.1 2.7.1 三次样条函数三次样条函数上是三次多项式,其中 是给定节点,bxxxan10 若函数 且在每个小区间 ,)(2
22、baCxS,1jjxx则称 是节点 上的三次样条函数.)(xSnxxx,10 若在节点 上给定函数值 jx), 1 ,0)(njxfyjjjjyxS)(), 1 ,0(nj(7.1)则称 为三次样条插值函数三次样条插值函数.)(xS定义定义4 4并成立45 由于 在每个小区间 上有4个待定系数,)(xS,1jjxx共有 个小区间,所以共有 个待定参数. nn4),0()0(jjxSxS 由于 在 上二阶导数连续,所以在节点 )( xS,ba)1,2, 1(njxj处应满足连续性条件),0()0(jjxSxS这些共有 个条件,再加上 本身还要满足的 个插值条件,共有 个条件,还需要2个才能确定
23、. )(xS)(xS33n24n1n).0()0( jjxSxS(7.2)46 通常可在区间 端点 上各加一个条件,banxbxa,0 1. 已知两端的一阶导数值,即.)(,)(00nnfxSfxS(7.3)(7.4)称为自然边界条件自然边界条件. 2. 已知两端的二阶导数,即其特殊情况为,)(,)(00nnfxSfxS (7.4).0)()(0 nxSxS(7.4)常见的边界条件有以下3种:(称为边界条件边界条件),47 3. 当 是以 为周期的周期函数时,则要求)(xf0 xxn),0()0(0nxSxS此时插值条件(7.1)中 . nyy0 这样确定的样条函数 称为周期样条函数周期样条函
24、数. .)(xS 这时边界条件应满足 ),0()0(0nxSxS(7.5)).0()0(0 nxSxS)(xS也是周期函数.48 2.7.2 2.7.2 样条插值函数的建立样条插值函数的建立 下面利用 的二阶导数值 表示 . )(xSjjMxS )(), 1 , 0(nj)(xS 由于 在区间 上是三次多项式,故 )(xS,1jjxx)(xS 在 上是线性函数,,1jjxxjjjjjjhxxMhxxMxS 11)((7.7)对 积分两次并利用 及 , )(xS jjyxS)(11)(jjyxS可表示为可定出积分常数,于是得三次样条表达式49jjjjjjhxxMhxxMxS6)(6)()(313
25、1jjjjjjjjjjhxxhMyhxxhMy6621112这里 ,是未知的. ), 1 ,0(njMj).1, 1 ,0(nj(7.8)50 为了确定 ,对 求导得 jM)(xSjjjjjjhxxMhxxMxS2)(2)()(2112;611jjjjjjhMMhyy(7.9)由此可求得 .63)0(11jjjjjjjjhyyMhMhxS51类似地可求出 在区间 上的表达式,从而得 )(xS,1jjxx,36)0(11111jjjjjjjjhyyMhMhxS利用 可得 )0()0(jjxSxS),1,2, 1(211njdMMMjjjjjj(7.10)其中 , 1 ,0,111njhhhhhh
26、jjjjjjjj52.,6,611111jjjjjjjjjjxxxfhhxxfxxfd(7.11)对第一种边界条件(7.3),可导出两个方程 ).,(62),(62111010010nnnnnnxxffhMMfxxfhMM(7.12)53如果令 ),(6, 1010000fxxfhd),(6, 111nnnnnnxxffhd.222211011011110nnnnnnnddddMMMM(7.13)那么(7.10)及(7.12)可写成矩阵形式 54 对第二种边界条件(7.4),直接得端点方程 ,00fM .nnfM (7.14)如果令 , nnnfdfd 2,2,0000也可以写成(7.13)的
27、矩阵形式.则(7.10)和(7.14)55 对于第三种边界条件(7.5),可得 其中 ,010hhhnn,601110hhxxfxxfdnnnn,0nMM,211nnnnndMMM(7.15)(7.10)和(7.15)可以写成矩阵形式,1011hhhnnnn56.2222121121112211nnnnnnnnddddMMMM(7.16) (7.13)和(7.16)是关于 的三对角方程组, 在力学上解释为细梁在 截面处的弯矩,称为 的矩,), 1 ,0(njMjjMjx)(xS (7.13)和(7.16)的系数矩阵为严格对角占优阵,有唯一解, 求解方法可见第5章第4节追赶法,将解得结果代入(7
28、.8)的表达式即可.方程组(7.13)和(7.16)称为三弯矩方程三弯矩方程.57 设 为定义在 上的函数,在节点 )(xf30,7.27)3,2, 1 ,0(ixi.0.3)30()(, 1.4)29()(,3.4)28()(, 1.4)7.27()(3210fxffxffxffxf试求三次样条函数 ,使它满足边界条件 )(xS,0.3)7.27(S.0.4)30(S例例5 5上的值如下:58,666.46),(6,21,1310, 101000210fxxfhd.4.17),(632323xxffhd,70000.2,6,00002.4,632122101xxxfdxxxfd由此得矩阵形式
29、的方程组(7.13)为 , 1,21,133, 1,30.0321210hhh 解解由(7.11)及(7.12)59.4000.177000.200002.4666.464000.17212122113102133123210MMMM求解得 .115.9,830.0,395.0,531.233210MMMM60代入(7.8)得)29(51917.4)29(51917.1)30(96167.3)30(13833.0),28(96167.3)28(13833.0)29(23417.4)29(06583.0),7 .27(31358.14)7 .27(21944.0)28(84322.14)28(0
30、7278.13)(333333xxxxxxxxxxxxxS,28,7.27x,29,28x,30,29x(曲线见图2-6)61图2-662 给定函数 节点 ,55,11)(2xxxf),10, 1 ,0(5kkxk用三次样条插值求).(10 xS 取)()(10kkxfxS),5()5(),10, 1 ,0(10fSk).5()5(10fS直接上机计算可求出 在表2-6所列各点的值. )(10 xS例例6 66325376. 013971. 013793. 05 . 20000. 10000. 10000. 1019837. 011366. 011312. 08 . 294090. 09275
31、4. 091743. 03 . 010000. 010000. 010000. 00 . 384340. 082051. 080000. 05 . 010832. 008426. 008410. 03 . 364316. 062420. 060976. 08 . 022620. 007606. 007547. 05 . 350000. 050000. 050000. 00 . 120130. 006556. 006477. 08 . 331650. 036133. 037175. 03 . 105882. 005882. 005882. 00 . 423535. 029744. 030769.
32、 05 . 188808. 004842. 005131. 03 . 418878. 023154. 023585. 08 . 157872. 104248. 004706. 05 . 420000. 020000. 020000. 00 . 280438. 103758. 004160. 08 . 424145. 016115. 015898. 03 . 203846. 003846. 003846. 00 . 5)()(11)()(111010210102xLxSxxxLxSxx6表264下图是用Matlab完成的样条插值(附程序):65附:样条插值程序n=11; m=61;x=-5:10
33、/(m-1):5;y=1./(1+x.2);z=0*x;x0=-5:10/(n-1):5;y0=1./(1+x0.2);y1=interp1(x0, y0, x, spline);plot(x, z, r, x, y, k:, x, y1, r)gtext(Spline), gtext(y=1/(1+x2)title(Spline)注:interp1(x0, y0, x, spline)为Matlab中现成的样条插值程序.66也可以将三种插值结果画在一起:67 2.7.3 2.7.3 误差界与收敛性误差界与收敛性 定理定理4 4设 为满足第一种或第二,)(4baCxf)( xS则有估计式),1,1 ,0(,max110nixxhhhiiiini,)(max)()(max4)4()()(kbxakkkbxahxfCxSxf,2,1 ,0k(7.17)其中.83,241,3845210CCC种边界条件(7.3)或(7.4)的三次样条函数,令68 这个定理不但给出了三次样条插值函数 的误差估计.)( xS还说明了当 时, 及其一阶导数 和0h)( xS)( xS 二阶导数 均分别一致收敛于 , 及 )( xS )( xf)(xf ).( xf