第二章插值法

1、第二章第二章 插值法插值法 /* Chapter 2 Interpolation */当精确函数当精确函数 y = f(x) 非常复杂或未知时，在一非常复杂或未知时，在一系列节点系列节点 x0 xn 处测得函数值处测得函数值 y0 = f(x0), yn = f(xn)，由此构造一个简单易算的近似函，由此构造一个简单易算的近似函数数 g(x) f(x)，满足条件，满足条件g(xi) = f(xi) (i = 0, n)。这里的。这里的 g(x) 称为称为f(x) 的的插值函数插值函数。最常。最常用的插值函数是用的插值函数是 ?多项式多项式x0 x1x2x3x4xg(x) f(x)2.1 多项式

2、插值多项式插值 /* Polynomial Interpolation */2.2 拉格朗日多项式拉格朗日多项式 /* Lagrange Polynomial */niyxPiin,., 0,)(= = =求求 n 次多项式次多项式 使得使得nnnxaxaaxP = =10)(条件：条件：无重合节点，即无重合节点，即jixx ji n = 1已知已知 x0 , x1 ; y0 , y1 ，求，求xaaxP101)( = =使得使得111001)(,)(yxPyxP= = =可见可见 P1(x) 是过是过 ( x0 , y0 ) 和和 ( x1, y1 ) 两点的直线。两点的直线。)()(001

3、0101xxxxyyyxP- - - - = =101xxxx- - -010 xxxx- - -= y0 + y1l0(x)l1(x) = = =10)(iiiyxl称为称为拉氏基函数拉氏基函数 /* Lagrange Basis */，满足条件满足条件 li(xj)= ij /* Kronecker Delta */线性插值线性插值n = 2已知已知 x0 , x1 , x2; y0 , y1 , y2 ，满足，满足()y00,xf= =抛物线插值抛物线插值，写出二次拉格朗日插值，写出二次拉格朗日插值11)(yxf= =,22)(yxf= =()()( )()()xxxxlxxxxx-=-

4、02110120122021()()()()()xxxxlxxxxx- - -= =- - -多项式：多项式： = = =20)(iiiyxl)(xL22.2 Lagrange Polynomial()()( )()()xxxxlxxxxx-=-12001022.2 Lagrange Polynomialn 1希望找到希望找到li(x)，i = 0, , n 使得使得 li(xj)= ij ；然后令；然后令 = = =niiinyxlxP0)()(，则显然有，则显然有Pn(xi) = yi 。li(x)每个每个 li 有有 n 个根个根 x0 xi xn = =- -= =- - - -= =

5、njj i jiniiixxCxxxxxxCxl00)().().()( - -= = =j i jiiiixxCxl)(11)(= = - - -= =njijjijixxxxxl0)()()( = = =niiinyxlxL0)()(Lagrange Polynomial与与 有关，而与有关，而与 无关无关节点节点f定理定理 (唯一性唯一性) 满足满足 的的 n 阶插值多阶插值多项式是唯一存在的。项式是唯一存在的。niyxPii,., 0,)(= = =证明：证明：反证：若不唯一，则除了反证：若不唯一，则除了Ln(x) 外还有另一外还有另一 n 阶多项阶多项式式 Pn(x) 满足满足 Pn

6、(xi) = yi 。考察考察 则则 Qn 的阶数的阶数, )()()(xLxPxQnnn- -= = n而而 Qn 有有 个不同的根个不同的根n + 1x0 xn注：注：若不将多项式次数限制为若不将多项式次数限制为 n ，则插值多项式，则插值多项式不唯一不唯一。例如例如 也是一个插值也是一个插值多项式，其中多项式，其中 可以是任意多项式。可以是任意多项式。= =- - = =niinxxxpxLxP0)()()()()(xp2.2 Lagrange Polynomial 插值余项插值余项 /* Remainder */设节点设节点)1( nf在在a , b内存在内存在, 考察截断误差考察截断

7、误差)()()(xLxfxRnn- -= =, baCfn bxxxan 10，且，且 f 满足条件满足条件 ,Rolles Theorem: 若若 充分光滑，充分光滑， ，则，则存在存在 使得使得 。)(x 0)()(10= = =xx ),(10 xx 0)(= = 推广：推广：若若0)()()(210= = = =xxx ),(),(211100 xxxx 使得使得0)()(10= = = = ),(10 使得使得0)(= = ),(, 0)()1(baxxn = = 0)()(0= = = =nxx 2.2 Lagrange PolynomialRn(x) 至少有至少有 个根个根n+1

8、 = =- -= =niinxxxKxR0)()()(任意固定任意固定 x xi (i = 0, , n), 考察考察 = =- - -= =niixtxKtRnt0)()()()( (x)有有 n+2 个不同的根个不同的根 x0 xn x),(, 0)()1(baxxn = = !)1()()()1(-nxKRxnn 注意这里是对注意这里是对 t 求导求导= = - - - !)1)()()()1()1(nxKLfxnnxn !)1()()()1( = = nfxKxn = = - - = =niixnnxxnfxR0)1()(! ) 1()()( 2.2 Lagrange Polynomi

9、al注：注： 通常不能确定通常不能确定 x , 而是估计而是估计 , x (a,b) 将将 作为误差估计上限。作为误差估计上限。1)1()( nnMxf= = - - niinxxnM01|)!1(当当 f(x) 为任一个次数为任一个次数 n 的的多项式多项式时，时， , 可知可知 ，即插值多项式对于次数，即插值多项式对于次数 n 的的多项多项式是式是精确精确的。的。0)()1( xfn0)( xRn2.2 Lagrange Polynomial注：注： 小的区间上插值有利于减少误差；小的区间上插值有利于减少误差； 依靠增多插值节点不一定能减少误差；依靠增多插值节点不一定能减少误差； 多项式插

10、值，外推误差可能要比内插误差大。多项式插值，外推误差可能要比内插误差大。例：例：已知已知233sin,214sin,216sin= = = = 分别利用分别利用 sin x 的的1次、次、2次次 Lagrange 插值计算插值计算 sin 50 并估计误差。并估计误差。 解：解：0 x1x2x185500 =n = 1分别利用分别利用x0, x1 以及以及 x1, x2 计算计算4,610 =xx利用利用216/4/6/214/6/4/)(1 - - - - - -= = xxxL这里这里(2)( )sin ,( )sin,( , )6 4xxxf xx f =-而而(2)1()12sin,(

11、 )()()222!64xxfR xxx=-00762. 0)185(01319. 01- - - - Rsin 50 = 0.7660444)185(50sin10 L0.77614外推外推 /* extrapolation */ 的实际误差的实际误差 - -0.010010.010013,421 = = =xx利用利用sin 50 0.76008, 00660. 018500538. 01 R内插内插 /* interpolation */ 的实际误差的实际误差 0.00596 0.00596内插通常优于外推。选择内插通常优于外推。选择要计算的要计算的 x 所在的区间的所在的区间的端点，插

12、值效果较好。端点，插值效果较好。2.2 Lagrange Polynomialn = 223)()(21)()(21)()()(4363463464363646342 - - - - - - - - - - - - - - -= = xxxxxxxL)185(50sin20 L0.7654323cos21;)3)(4)(6(!3cos)(2 - - - - -= =xxxxxxR 00077. 018500044. 02 Rsin 50 = 0.76604442次插值的实际误差次插值的实际误差 0.00061 0.00061高次插值通常优于高次插值通常优于低次插值低次插值但绝对不是次数越但绝对

13、不是次数越高就越好高就越好2.2 Lagrange Polynomial When you start writing the program, you will find how easy it is to calculate the Lagrange polynomial.Oh yeah? What if I find the current interpolation not accurate enough? Then you might want to take more interpolating points into account.Right. Then all the Lag

14、range basis, li(x), will have to be re-calculated. Excellent point !We will come to discuss this problemnext time.2.3 牛顿插值牛顿插值 /* Newtons Interpolation */Lagrange 插值虽然易算，但若要增加一个节点时，插值虽然易算，但若要增加一个节点时，全部基函数全部基函数 li(x) 都需重新算过。都需重新算过。的形式，希望每加一个节点时，的形式，希望每加一个节点时，将将 Ln(x) 改写成改写成.)()(102010 - - - - - xxxxa

15、xxaa).(10- - - - nnxxxxa只附加一项只附加一项上去即可。上去即可。? 差商差商( (亦称均差亦称均差) ) /* divided difference */),()()(,jijijijixxjixxxfxfxxf - - -= =1阶差商阶差商 /* the 1st divided difference of f w.r.t. xi and xj */)(,kixxxxfxxfxxxfkikjjikji - - -= =2阶差商阶差商2.3 Newtons Interpolation11101010111010,.,.,.,.,., - - - - - -= =- -

16、-= =kkkkkkkkkkkxxxxxfxxxfxxxxxfxxxfxxf(k+1)阶差商：阶差商： = = = =kiikikxxfxxf010)()(,., 事实上事实上其中其中,)()(01= = - -= =kiikxxx = = - -= = kijjjiikxxx01)()( Warning: my head is explodingWhat is the point of this formula?差商的值与差商的值与 xi 的顺序无关！的顺序无关！P30均差的性质均差的性质 牛顿插值牛顿插值 /* Newtons Interpolation */,)()()(000 xxfx

17、xxfxf- - = =,)(,101100 xxxfxxxxfxxf- - = =,.,)(,.,.,0010nnnnxxxfxxxxfxxxf- - = =- -).(.)()()(10102010- - - - - - - - - = =nnnxxxxaxxxxaxxaaxN12 n- -11+ (x - - x0) 2+ + (x - - x0)(x - - xn- -1) n- -1.)(,)(,)()(102100100 - - - - - = =xxxxxxxfxxxxfxfxf).(,.,100- - - - nnxxxxxxf)().(,.,100nnnxxxxxxxxxf-

18、 - - - - -Nn(x)-牛顿牛顿插值多项式插值多项式Rn(x)ai = f x0, , xi 2.3 Newtons Interpolation注：注： 由由唯一性可知唯一性可知 Nn(x) Ln(x)， 只是算法不同，故其只是算法不同，故其余项也相同，即余项也相同，即)(!)1()()(,.,1)1(10 xnfxxxxfkxnkn = = ),(,!)(,.,maxmin)(0 xxkfxxfkk = = 实际计算过程为实际计算过程为f (x0)f (x1)f (x2)f (xn- -1)f (xn)f x0, x1f x1, x2 f xn- -1, xnf x0, x1 , x

19、2 f xn- -2, xn- -1, xnf x0, , xn f (xn+1) f xn, xn+1 f xn- -1, xn, xn+1 f x1, , xn+1 f x0, , xn+12.3 Newtons Interpolation2.3 Newtons Interpolation例例4 (P32): 给出给出f (x)的函数表，求的函数表，求4次牛顿插值多项式，并由此次牛顿插值多项式，并由此计算计算f (0.596)的近似值。的近似值。N4(x)=0.41075+1.116(x-0.4)+0.28(x-0.4)(x-0.55)+0.19733(x-0.4)(x-0.55)(x-0

20、.65) +0.03134 (x-0.4)(x-0.55)(x-0.65)(x-0.8)f(0.596)N4(0.596)=0.63192R4(x)=fx,x0, xn (x-0.4)(x-0.55)(x-0.65)(x-0.8)(x-0.9) R4(0.596)=? 差分形式差分形式等距节点公式等距节点公式 /* Formulae with Equal Spacing */向前差分向前差分 iiifff- -= = 1ikikikikffff1111)(- - - - - - - = = = = 向后差分向后差分 111- - - - - - = = ikikikfffi- -1iifff-

21、 -= = 中心差分中心差分 212111- - - - - -= =ikikikfff 其中其中)(221hiixff = = 当节点当节点等距等距分布时分布时:),.,0(0nihixxi= = = =2.3 Newtons Interpolation 差分的重要性质：差分的重要性质： 线性：例如线性：例如gbfaxgbxfa = = )()( 若若 f (x)是是 m 次多项式，则次多项式，则 是是 次多项次多项式，而式，而 )0()(mkxfk km -)(0)(mkxfk = = 差分值可由函数值算出：差分值可由函数值算出： = =- - - -= = njjknjknfjnf0)1

22、( = =- - - - -= = njnjkjnknfjnf0) 1(!)1).(1(jjnnnjn - - -= = 其中其中/* binomial coefficients */ 函数值可由差分值算出：函数值可由差分值算出：kjnjknfjnf =0kkkhkfxxf!,.,00 = =knkknnnhkfxxxf!,.,1 = =- - -kkkhff0)()( = = 由由 Rn 表达式表达式2.3 Newtons Interpolation 牛顿公式牛顿公式 ).(,.,.)(,)()(1000100- - - - - - = =nnnxxxxxxfxxxxfxfxN 牛顿前插公式

23、牛顿前插公式 /* Newtons forward-difference formula */ 牛顿后插公式牛顿后插公式 /* Newtons backward-difference formula */将节点顺序倒置：将节点顺序倒置：).(,.,.)(,)()(101xxxxxxfxxxxfxfxNnnnnnnn- - - - - = =- -设设htxx = =0，则，则)()()(000 xfkthtxNxNknknn = = = = = =),(,).(1()!1()()(01)1(nnnnxxhntttnfxR - - - = = 设设htxxn = =，则，则)() 1()()(0

24、nknkknnnxfkthtxNxN - - -= = = = = =注：注：一般当一般当 x 靠近靠近 x0 时用前插，靠近时用前插，靠近 xn 时用后插，故两时用后插，故两种公式亦称为种公式亦称为表初公式表初公式和和表末公式表末公式。2.3 Newtons Interpolation2.4 埃尔米特插值埃尔米特插值 /* Hermite Interpolation */不仅要求函数值重合，而且要求若干阶不仅要求函数值重合，而且要求若干阶导数导数也重合。也重合。即：要求插值函数即：要求插值函数 (x) 满足满足 (xi) = f (xi), (xi) = f (xi), (mi) (xi)

25、= f (mi) (xi).注：注： N 个条件可以确定个条件可以确定 阶多项式。阶多项式。N - - 1要求在要求在1个节点个节点 x0 处直到处直到m0 阶导数都重合的插阶导数都重合的插值多项式即为值多项式即为Taylor多项式多项式00)(!)(.)()()(000)(000mmxxmxfxxxfxfx- - - - = = 其余项为其余项为)1(00)1(00)()!1()()()()(-=-=mmxxmfxxfxR 一般只考虑一般只考虑 f 与与f 的值。的值。2.4 Hermite Interpolation例：例：设设 x0 x1 x2, 已知已知 f(x0)、 f(x1)、 f

26、(x2) 和和 f (x1), 求多项式求多项式 P(x) 满足满足 P(xi) = f (xi)，i = 0, 1, 2，且且 P(x1) = f (x1), 并估计误差。并估计误差。模仿模仿 Lagrange 多项式的思想，设多项式的思想，设解：解：首先，首先，P 的阶数的阶数 = 3=213)()()()()(=0iiixhx1f xhxfxP h0(x)有根有根x1, x2，且且 h0(x1) = 0 x1 是重根。是重根。)()()(22100 xxxxCxh- - -= =又又: h0(x0) = 1 C0 )()()()()(202102210 xxxxxxxxxh- - - -

27、 -= =h2(x)h1(x)有根有根 x0, x2 )()()(201xxxxBAxxh- - - = =由余下条件由余下条件 h1(x1) = 1 和和 h1(x1) = 0 可解。可解。与与h0(x) 完全类似。完全类似。 (x) h1有根有根 x0, x1, x2 h1)()()(2101xxxxxxCx- - - -= = h1又又: (x1) = 1 C1 可解。可解。其中其中 hi(xj) = ij , hi(x1) = 0, (xi) = 0, (x1) = 1 h1 h1),()()()()()(221033xxxxxxxKxPxfxR- - - -= =- -= =!4)(

28、)()4(xfxK = =与与 Lagrange 分析分析完全类似完全类似一般地，已知一般地，已知 x0 , , xn 处有处有 y0 , , yn 和和 y0 , , yn ，求，求 H2n+1(x) 满足满足 H2n+1(xi) = yi ， H2n+1(xi) = yi。解：解：设设=ni)()()(=0iixhxhyixH2n+1 n=0iyi其中其中 hi(xj) = ij , hi(xj) = 0, (xj) = 0, (xj) = ij hi hihi(x)有根有根 x0 , , xi , , xn且都是且都是2重根重根 )()()(2xlBxAxhiiii = = - - -=

29、 =ijjijixxxxxl)()()(由余下条件由余下条件 hi(xi) = 1 和和 hi(xi) = 0 可解可解Ai 和和 Bi )()(21 )(2xlxxxlxhiiiii- - - -= = (x) hi有根有根 x0 , , xn, 除了除了xi 外都是外都是2重根重根 hi)()(iili2(x)xxCx- -= = hi又又: (xi) = 1 Ci = 1 hi)(x)(ili2(x)xx- -= =设设,.210baCfbxxxann = = = =则则20)22()()!22()()( - - = = = niixnnxxnfxR 这样的这样的Hermite 插值唯插

30、值唯一一 埃尔米特插值埃尔米特插值构造基函数的方法构造基函数的方法2.4 Hermite InterpolationQuiz: 给定给定 xi = i +1, i = 0, 1, 2, 3, 4, 5. 下面哪个是下面哪个是 h2(x)的图像？的图像？ x0-10.5123456yxy0-10.5123456斜率斜率=1 求求Hermite多项式的基本步骤：多项式的基本步骤： 写出相应于条件的写出相应于条件的hi(x)、 hi(x) 的组合形式；的组合形式； 对每一个对每一个hi(x)、 hi(x) 找出尽可能多的条件给出的根；找出尽可能多的条件给出的根； 根据多项式的总阶数和根的个数写出表达

31、式；根据多项式的总阶数和根的个数写出表达式； 根据尚未利用的条件解出表达式中的待定系数；根据尚未利用的条件解出表达式中的待定系数； 最后完整写出最后完整写出H(x)。2.4 Hermite Interpolation2.5 分段低次插值分段低次插值 /* piecewise polynomial approximation */Remember what I have said? Increasing the degree of interpolating polynomial will NOT guarantee a good result, since high-degree polyno

32、mials are oscillating.例：例：在在- -5, 5上考察上考察 的的Ln(x)。取。取211)(xxf=),., 0(105niinxi= = - -= = -5 -4 -3 -2 -1 0 1 2 3 4 5 -0.5 0 0.5 1 1.5 2 2.5 n 越大，越大，端点附近抖动端点附近抖动越大，称为越大，称为Runge 现象现象Ln(x) f (x) 分段分段低次低次插值插值2.5 Piecewise Polynomial Approximation 分段线性插值分段线性插值 /* piecewise linear interpolation */在每个区间在每个区

33、间 上，用上，用1阶多项式阶多项式 (直线直线) 逼近逼近 f (x):,1 iixx11111)()( - - - - - -= = iiiiiiiiyxxxxyxxxxxPxf,for 1 iixxx记记 ，易证：当，易证：当 时，时，|max1iixxh- -= = 0h)()(1xfxPh一致一致失去了原函数的光滑性。失去了原函数的光滑性。 分段三次分段三次Hermite插值插值 /* Hermite piecewise polynomials */给定给定nnnyyyyxx ,.,;,.,;,.,000在在 上利用两点的上利用两点的 y 及及 y 构造构造3次次Hermite函数函数

34、,1 iixx导数一般不易得到。导数一般不易得到。How can we make a smooth interpolation without asking too much from f ?Headache 2.6 三次样条插值三次样条插值 /* Cubic Spline */定义定义设设 。三次样条函数三次样条函数 , 且在每个且在每个 上为上为三次多项式三次多项式 /* cubic polynomial */。若它同。若它同时还满足时还满足 ，则称为，则称为 f 的的三次样条插值函三次样条插值函数数 /* cubic spline interpolant */.bxxxan= = = =

35、.10,)(2baCxS ,1 iixx),., 0(),()(nixfxSii= = =注：注：三次样条与分段三次样条与分段 Hermite 插值的根本区别在于插值的根本区别在于S(x)自自身光滑身光滑，不需要知道，不需要知道 f 的导数值（除了在的导数值（除了在2个端点可能需个端点可能需要）；而要）；而Hermite插值依赖于插值依赖于f 在所有插值点的导数值。在所有插值点的导数值。f(x)H(x)S(x)2.6 Cubic Spline 构造三次样条插值函数的构造三次样条插值函数的三弯矩法三弯矩法 /* method of bending moment */在在 上，记上，记,1jjxx

36、- -,1- - -= =jjjxxh,for )()(1jjjxxxxSxS- - = =对每个对每个j, 此为此为3次多项式次多项式则则 Sj”(x) 为为 次多项式，需次多项式，需 个点的值确定之。个点的值确定之。12设设 Sj”(xj- -1) = Mj- -1， Sj”(xj) = Mj 对应力学中的对应力学中的梁弯矩梁弯矩，故名，故名对于对于x xj- -1, xj 可可得到得到Sj”(x) =jjjjjjhxxMhxxM11- - - - - -积分积分2次，可得次，可得 Sj(x) 和和 Sj(x) ：jjjjjjjAhxxMhxxM - - - - - - - -2)(2)(

37、21121Sj(x) =jjjjjjjjBxAhxxMhxxM - - - - - -6)(6)(3131Sj(x) =利用已知利用已知Sj(xj- -1) = yj- -1 Sj(xj) = yj 可解可解jjjjjjjhMMhyyA611- - - - - -= =jjjjjjjjjjjjhxxhMyhxxhMyBxA12211)6()6(- - - - - - - - -= = 下面解决下面解决 Mj ： 利用利用S 在在 xj 的的连续性连续性xj- -1, xj : Sj(x) =jjjjjjjjjjjhMMxxfhxxMhxxM6,2)(2)(112121- - - - - - -

38、 - - - - -1111211216,2)(2)( - - - - - - - -jjjjjjjjjjjhMMxxfhxxMhxxMxj , xj+1: Sj+1(x) =利用利用Sj(xj) = Sj+1(xj)，合并关于，合并关于Mj- -1、 Mj、 Mj+1的同类项，并的同类项，并记记 , , , 整理整理后得到：后得到：11jjjjhhh=l l1jj-=l lm m),(6111jjjjjjjxxfxxfhhg-=211gMMMjjjjjj=-l lm m j 1n- -1即：有即：有 个未知数，个未知数， 个方程。个方程。n- -1n+1 = = - - - -110111122nnnnggMMl lm ml lm m还需还需2个个边界条件边界条件 /* boundary conditions */2.6 Cubic Spline 第第1类边条件类边条件 /* clamped boundary */： S(a) = y0， S(b) = yna , x1 : S1(x) =1011012112106,2)(2)(hMMxxfhaxMhxxM- - - -