第二章插值法

1、第二章第二章 插值法插值法 /* Chapter 2 Interpolation */當精確函數(shù)當精確函數(shù) y = f(x) 非常復雜或未知時，在一非常復雜或未知時，在一系列節(jié)點系列節(jié)點 x0 xn 處測得函數(shù)值處測得函數(shù)值 y0 = f(x0), yn = f(xn)，由此構造一個簡單易算的近似函，由此構造一個簡單易算的近似函數(shù)數(shù) g(x) f(x)，滿足條件，滿足條件g(xi) = f(xi) (i = 0, n)。這里的。這里的 g(x) 稱為稱為f(x) 的的插值函數(shù)插值函數(shù)。最常。最常用的插值函數(shù)是用的插值函數(shù)是 ?多項式多項式x0 x1x2x3x4xg(x) f(x)2.1 多項式

2、插值多項式插值 /* Polynomial Interpolation */2.2 拉格朗日多項式拉格朗日多項式 /* Lagrange Polynomial */niyxPiin,., 0,)(= = =求求 n 次多項式次多項式 使得使得nnnxaxaaxP = =10)(條件：條件：無重合節(jié)點，即無重合節(jié)點，即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稱為稱為拉氏基函數(shù)拉氏基函數(shù) /* 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與與 有關，而與有關，而與 無關無關節(jié)點節(jié)點f定理定理 (唯一性唯一性) 滿足滿足 的的 n 階插值多階插值多項式是唯一存在的。項式是唯一存在的。niyxPii,., 0,)(= = =證明：證明：反證：若不唯一，則除了反證：若不唯一，則除了Ln(x) 外還有另一外還有另一 n 階多項階多項式式 Pn(x) 滿足滿足 Pn

6、(xi) = yi ?？疾炜疾?則則 Qn 的階數(shù)的階數(shù), )()()(xLxPxQnnn- -= = n而而 Qn 有有 個不同的根個不同的根n + 1x0 xn注：注：若不將多項式次數(shù)限制為若不將多項式次數(shù)限制為 n ，則插值多項式，則插值多項式不唯一不唯一。例如例如 也是一個插值也是一個插值多項式，其中多項式，其中 可以是任意多項式?？梢允侨我舛囗検?。= =- - = =niinxxxpxLxP0)()()()()(xp2.2 Lagrange Polynomial 插值余項插值余項 /* Remainder */設節(jié)點設節(jié)點)1( nf在在a , b內(nèi)存在內(nèi)存在, 考察截斷誤差考察截斷

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) 為任一個次數(shù)為任一個次數(shù) n 的的多項式多項式時，時， , 可知可知 ，即插值多項式對于次數(shù)，即插值多項式對于次數(shù) n 的的多項多項式是式是精確精確的。的。0)()1( xfn0)( xRn2.2 Lagrange Polynomial注：注： 小的區(qū)間上插值有利于減少誤差；小的區(qū)間上插值有利于減少誤差； 依靠增多插值節(jié)點不一定能減少誤差；依靠增多插值節(jié)點不一定能減少誤差； 多項式插

10、值，外推誤差可能要比內(nèi)插誤差大。多項式插值，外推誤差可能要比內(nèi)插誤差大。例：例：已知已知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內(nèi)插內(nèi)插 /* interpolation */ 的實際誤差的實際誤差 0.00596 0.00596內(nèi)插通常優(yōu)于外推。選擇內(nèi)插通常優(yōu)于外推。選擇要計算的要計算的 x 所在的區(qū)間的所在的區(qū)間的端點，插

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高次插值通常優(yōu)于高次插值通常優(yōu)于低次插值低次插值但絕對不是次數(shù)越但絕對

13、不是次數(shù)越高就越好高就越好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 插值雖然易算，但若要增加一個節(jié)點時，插值雖然易算，但若要增加一個節(jié)點時，全部基函數(shù)全部基函數(shù) li(x) 都需重新算過。都需重新算過。的形式，希望每加一個節(jié)點時，的形式，希望每加一個節(jié)點時，將將 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均差的性質(zhì)均差的性質(zhì) 牛頓插值牛頓插值 /* 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)的函數(shù)表，求的函數(shù)表，求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)=? 差分形式差分形式等距節(jié)點公式等距節(jié)點公式 /* Formulae with Equal Spacing */向前差分向前差分 iiifff- -= = 1ikikikikffff1111)(- - - - - - - = = = = 向后差分向后差分 111- - - - - - = = ikikikfffi- -1iifff-

21、 -= = 中心差分中心差分 212111- - - - - -= =ikikikfff 其中其中)(221hiixff = = 當節(jié)點當節(jié)點等距等距分布時分布時:),.,0(0nihixxi= = = =2.3 Newtons Interpolation 差分的重要性質(zhì)：差分的重要性質(zhì)： 線性：例如線性：例如gbfaxgbxfa = = )()( 若若 f (x)是是 m 次多項式，則次多項式，則 是是 次多項次多項式，而式，而 )0()(mkxfk km -)(0)(mkxfk = = 差分值可由函數(shù)值算出：差分值可由函數(shù)值算出： = =- - - -= = njjknjknfjnf0)1

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

23、牛頓前插公式 /* Newtons forward-difference formula */ 牛頓后插公式牛頓后插公式 /* Newtons backward-difference formula */將節(jié)點順序倒置：將節(jié)點順序倒置：).(,.,.)(,)()(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 */不僅要求函數(shù)值重合，而且要求若干階不僅要求函數(shù)值重合，而且要求若干階導數(shù)導數(shù)也重合。也重合。即：要求插值函數(shù)即：要求插值函數(shù) (x) 滿足滿足 (xi) = f (xi), (xi) = f (xi), (mi) (xi)

25、= f (mi) (xi).注：注： N 個條件可以確定個條件可以確定 階多項式。階多項式。N - - 1要求在要求在1個節(jié)點個節(jié)點 x0 處直到處直到m0 階導數(shù)都重合的插階導數(shù)都重合的插值多項式即為值多項式即為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 的階數(shù)的階數(shù) = 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、值唯一一 埃爾米特插值埃爾米特插值構造基函數(shù)的方法構造基函數(shù)的方法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) 找出盡可能多的條件給出的根；找出盡可能多的條件給出的根； 根據(jù)多項式的總階數(shù)和根的個數(shù)寫出表達

31、式；根據(jù)多項式的總階數(shù)和根的個數(shù)寫出表達式； 根據(jù)尚未利用的條件解出表達式中的待定系數(shù)；根據(jù)尚未利用的條件解出表達式中的待定系數(shù)； 最后完整寫出最后完整寫出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 現(xiàn)象現(xiàn)象Ln(x) f (x) 分段分段低次低次插值插值2.5 Piecewise Polynomial Approximation 分段線性插值分段線性插值 /* piecewise linear interpolation */在每個區(qū)間在每個區(qū)

33、間 上，用上，用1階多項式階多項式 (直線直線) 逼近逼近 f (x):,1 iixx11111)()( - - - - - -= = iiiiiiiiyxxxxyxxxxxPxf,for 1 iixxx記記 ，易證：當，易證：當 時，時，|max1iixxh- -= = 0h)()(1xfxPh一致一致失去了原函數(shù)的光滑性。失去了原函數(shù)的光滑性。 分段三次分段三次Hermite插值插值 /* Hermite piecewise polynomials */給定給定nnnyyyyxx ,.,;,.,;,.,000在在 上利用兩點的上利用兩點的 y 及及 y 構造構造3次次Hermite函數(shù)函數(shù)

34、,1 iixx導數(shù)一般不易得到。導數(shù)一般不易得到。How can we make a smooth interpolation without asking too much from f ?Headache 2.6 三次樣條插值三次樣條插值 /* Cubic Spline */定義定義設設 。三次樣條函數(shù)三次樣條函數(shù) , 且在每個且在每個 上為上為三次多項式三次多項式 /* cubic polynomial */。若它同。若它同時還滿足時還滿足 ，則稱為，則稱為 f 的的三次樣條插值函三次樣條插值函數(shù)數(shù) /* cubic spline interpolant */.bxxxan= = = =

35、.10,)(2baCxS ,1 iixx),., 0(),()(nixfxSii= = =注：注：三次樣條與分段三次樣條與分段 Hermite 插值的根本區(qū)別在于插值的根本區(qū)別在于S(x)自自身光滑身光滑，不需要知道，不需要知道 f 的導數(shù)值（除了在的導數(shù)值（除了在2個端點可能需個端點可能需要）；而要）；而Hermite插值依賴于插值依賴于f 在所有插值點的導數(shù)值。在所有插值點的導數(shù)值。f(x)H(x)S(x)2.6 Cubic Spline 構造三次樣條插值函數(shù)的構造三次樣條插值函數(shù)的三彎矩法三彎矩法 /* 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 的的連續(xù)性連續(xù)性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即：有即：有 個未知數(shù)，個未知數(shù)， 個方程。個方程。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- - - -