數值分析第二章 插值

1、第二章 插值方法/* Interpolation */Interpolation_introduction1 問題提出問題提出函數逼近函數逼近 / /* *problem formulation-function approximationproblem formulation-function approximation* */ /用用Interpolation_introduction函數逼近的方法有很多，例如函數逼近的方法有很多，例如Taylor級數，級數，Fourier級數，有限元方法、邊界元方法，小級數，有限元方法、邊界元方法，小波分析等，大學科叫波分析等，大學科叫逼近論逼近論。本書

2、討論連續(xù)函數的逼近，主要介紹本書討論連續(xù)函數的逼近，主要介紹插值法插值法(chapter 2)和和最佳一直逼近、最小平方逼近最佳一直逼近、最小平方逼近離散數據擬合離散數據擬合(chapter 3)Interpolation_introduction插值節(jié)點插值節(jié)點插值條件插值條件-插值問題插值問題多項式插值是數值分析的基本工具，常用來計算被插函數多項式插值是數值分析的基本工具，常用來計算被插函數的近似的近似函數值函數值，零、極點零、極點，導數、積分導數、積分（第四章（第四章 數值積數值積分和數值微分），分和數值微分），解微分方程解微分方程（第五章）、（第五章）、積分方程積分方程插值插值多項式插

3、值多項式插值-polynomial interpolationProblem I. 給定給定y=f(x)的函數表的函數表, xi a,bniyxPiin,., 0,)(= = =求求 次數不超過次數不超過 n 的多項式的多項式 使得使得nnnxaxaaxP = =10)(條件：條件：無重合節(jié)點無重合節(jié)點，即，即jixx ji Interpolation intervalInterpolation conditionInterpolation polynomialInterpolation pointsInterpolation polynomial (2.1)(2.2)x0 x1x2x3x4x

4、Pn(x) f(x)多項式插值的幾何意義多項式插值的幾何意義Interpolation polynomial 求求插值多項式的唯一性插值多項式的唯一性 提問：提問：Problem I 中的中的Pn(x)是否存在？是否存在？ 若存在，是否唯一？如何求？若存在，是否唯一？如何求？Interpolation polynomial Interpolation polynomial 如何求？解線性方程如何求？解線性方程組（組（2.3）-待定系待定系數法數法Interpolation polynomial 2 拉格朗日多項式拉格朗日多項式 /* Lagrange Polynomial */niyxPiin

5、,., 0,)(= = =求求 n 次多項式次多項式 使得使得nnnxaxaaxP = =10)(條件：條件：無重合節(jié)點，即無重合節(jié)點，即jixx ji n = 1已知已知 x0 , x1 ; y0 , y1 ，求，求xaaxP101)( = =使得使得111001)(,)(yxPyxP= = =可見可見 P1(x) 是過是過 ( x0 , y0 ) 和和 ( x1, y1 ) 兩點的直線。兩點的直線。)()(0010101xxxxyyyxP- - - - = =101xxxx- - -010 xxxx- - -= y0 + y1l0(x)l1(x) = = =10)(iiiyxl稱為稱為拉氏

6、基函數拉氏基函數 /* Lagrange Basis */，滿足條件滿足條件 li(xj)= ij /* Kronecker Delta */2 Lagrange Polynomialn 1希望找到希望找到li(x)，i = 0, , n 使得使得 li(xj)= ij ；然后令；然后令，則顯然有，則顯然有Pn(xi) = yi 。li(x)每個每個 li 有有 n 個根個根 x0 xi xn= = =jiC0 = =- -nj i jxx)(- - - -inxxixxxxC0).().(ixl)( - -= =j i jxixiC)(1= =iixl1)(= = - - -= =njijj

7、ijixxxxxl0)()()( = = =niiinyxlxL0)()(Lagrange Polynomial與節(jié)點有關，而與與節(jié)點有關，而與 f無關無關 = = =niinxlxP0)()(yi 基函數法（基函數法（n=1情形的推廣情形的推廣）2 Lagrange Polynomial定理定理 (唯一性唯一性) 滿足滿足 的的 n 階插值多階插值多項式是唯一存在的。項式是唯一存在的。niyxPii,., 0,)(= = =證明證明： ( 前面已利用前面已利用Vandermonde 行列式行列式論證論證)反證：若不唯一，則除了反證：若不唯一，則除了Ln(x) 外還有另一外還有另一 n 階多項

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

9、在內存在, 考察截斷誤差考察截斷誤差)()()(xLxfxRnn- -= =, baCfn bxxxan 10，且，且 f 滿足條件滿足條件 ,Rolles Theorem: 若若 充分光滑，充分光滑， ，則，則存在存在 使得使得 。)(x 0)()(10= = =xx ),(10 xx 0)(= = 推廣：推廣：若若0)()()(210= = = =xxx ),(),(211100 xxxx 使得使得0)()(10= = = = ),(10 使得使得0)(= = 0)()(0= = = =nxx 存在存在),(ba 使得使得0)()(= = nRn(x) 至少有至少有 個根個根n+1 = =

10、- -= =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 Lagrange Polynomial1 La

11、grange Polynomial注：注： 通常不能確定通常不能確定 x , 而是估計而是估計 , x (a,b) 將將 作為誤差估計上限。作為誤差估計上限。1)1()( nnMxf= = - - niinxxnM01|)!1(當當 f(x) 為任一個次數為任一個次數 n 的的多項式多項式時，時， , 可知可知 ，即插值多項式對于次數，即插值多項式對于次數 n 的的多項多項式是式是精確精確的。的。0)()1( xfn0)( xRnQuiz: 給定給定 xi = i +1, i = 0, 1, 2, 3, 4, 5. 下面哪個是下面哪個是 l2(x)的圖像？的圖像？ y 0 - - - 1 0.

12、5 -0.5 1 2 3 4 5 6 x y 0 - - - 1 0.5 -0.5 1 2 3 4 5 6 x y 0 - - - 1 0.5 -0.5 1 2 3 4 5 6 x ABC 1 Lagrange Polynomial例：例：已知已知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/

13、)(1 - - - - - -= = xxxL這里這里)3,6(,sin)(,sin)()2( - -= = =xxxfxxf而而)4)(6(!2)()(,23sin21)2(1 - - -= = xxfxRxx00762. 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內插內插 /* interpol

14、ation */ 的實際誤差的實際誤差 0.005960.00596內插通常優(yōu)于外推。選擇內插通常優(yōu)于外推。選擇要計算的要計算的 x 所在的區(qū)間的所在的區(qū)間的端點，插值效果較好。端點，插值效果較好。1 Lagrange Polynomialn = 223)()(21)()(21)()()(4363463464363646342 - - - - - - - - - - - - - - -= = xxxxxxxL)185(50sin20 L0.7654323cos21;)3)(4)(6(!3cos)(2 - - - - -= =xxxxxxR 00077. 018500044. 02 Rsin 5

15、0 = 0.76604442次插值的實際誤差次插值的實際誤差 0.000610.00061高次插值通常優(yōu)于高次插值通常優(yōu)于低次插值低次插值但絕對不是次數越但絕對不是次數越高就越好，嘿高就越好，嘿嘿嘿 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

16、 more interpolating points into account.Right. Then all the Lagrange basis, li(x), will have to be re-calculated. Excellent point !We will come to discuss this problemnext time.3 逐次線性插值逐次線性插值 /* Lagrange Polynomial */ 0,1010011100,100100,202200,200200,20,10,1 20,1121:1Lagrange,;:1LagrangeIxx xxf xx

17、f xf xf xIxf xxxxxIxx xf xf xIxf xxxxxIxIxIxIxxxxx-=-=-=-，以 ，為節(jié)點的 次插值公式，實際上是過，點的直線，采用點斜式以 ，為節(jié)點的 次插值公式,同理有仿上，令易證 0,200,100,1 200,10010,100210,210,110,1 210,11110,111210,220,120,1 220,12210,222210,1 20122IxIxIxIxxxIxf xxxIxIxIxIxxxIxf xxxIxIxIxIxxxIxf xxxIxx x x-=-=-=-=-=-=-，由插值公式的唯一性知，是以 ，, 為節(jié)點的 次Lag

18、range插值多項式。以此 0,12,0,110,10,111101kkkkkkkkkIxIxIxIxxxxxx xxk-=-，類推，是以 ，, , 為節(jié)點的 次Lagrange插值多項式。 10,10,110,12,11kkkkkkkkkkx xx xIxIxIxxxxx-=-，實際上實際上，是對兩個低次插值的線性插值，這種通過低次插值再作是對兩個低次插值的線性插值，這種通過低次插值再作線性插值生成高次插值的方法稱為線性插值生成高次插值的方法稱為逐次線性插值逐次線性插值。 Aitken法法(按下表計算按下表計算) 0,12,0,110,10,1111kkkkkkkkIxIxIxIxxxxx-

19、=-，線性插值基函數線性插值基函數 0,0,1,0,1,2,0,1,2,3,00110,10220,20,1,2 kkkkkkxf xIxIxIxIxxf xxf xIxx xxf xIxIx- 1330,30,1,30,1,2,32440,40,1,40,1,2,40,1,2,3,43 x xxf xIxIxIxx xxf xIxIxIxIxx x-0123 x xx xx xx x-增加增加如果精度不構，增加節(jié)點如果精度不構，增加節(jié)點x4, ,同時表中增加一行，三同時表中增加一行，三角形斜邊上即為所要求的各次插值多項式。角形斜邊上即為所要求的各次插值多項式。k1k0k2k3k4 Nevil

20、le法法(按下表計算按下表計算) 1,1, ,11, ,11, ,1,200110,10221,20,1,2 kkkkkk kkk kkk kkxf xIxIxIxIxxf xxf xIxx xxf xIxIx- 0332,31,2,30,1,2,30443,42,3,40,1,2,40,1 x xxf xIxIxIxx xxf xIxIxIxI- ,2,3,401111 kkkkxx xx xx xx xx x-增加增加如果精度不構，增加節(jié)點如果精度不構，增加節(jié)點x4, ,同時表中增加一行，三角形斜邊上同時表中增加一行，三角形斜邊上即為所要求的各次插值多項式。即為所要求的各次插值多項式。k1

21、k0k1k1k1 11,0,110,10,1100kkkkkkIxIxIxIxxxxx-=-，HW: 用類似于前面的方法用類似于前面的方法構造構造Neville計算公式計算公式注：注：Atkin方法和方法和Neville方法與方法與Lagrange公式相比，當公式相比，當需要增加節(jié)點時，很容易由低次插值構造高次插值，而需要增加節(jié)點時，很容易由低次插值構造高次插值，而Lagrange插值公式中，每個基函數都需要作適當變化。插值公式中，每個基函數都需要作適當變化。誤差估計：誤差估計：由插值多項式的存在唯一性知，仍有由插值多項式的存在唯一性知，仍有= = - - = =niixnnxxnfxR0)1

22、()(! ) 1()()( 但這里可采用一種更簡便的方法。當但這里可采用一種更簡便的方法。當f (n+1)(x)在插在插值區(qū)間變化不大時，設值區(qū)間變化不大時，設f (n+1)(x)L，則有則有 0,110-10,12,0-2-!-!kkkkkkLf xIxx xx xkLf xIxx xx xx xk-， 10,110,12,0,1110,12,0,11 if kkkkkkkkkkxxf xIxIxIxxxIxIx-，可認為可認為 滿足精度要求。滿足精度要求。 0,11kf xIx-， 0,11-10,12,-kkkkkf xIxx xf xIxx x-，根據前面的計算結果估計當根據前面的計算

23、結果估計當前的誤差：事后誤差估計前的誤差：事后誤差估計（實用），前面給出的誤差（實用），前面給出的誤差估計（事先誤差估計）不實估計（事先誤差估計）不實用用HW: p.58-59#1-84 牛頓插值牛頓插值 /* Newtons Interpolation */Lagrange 插值雖然易算，但若要增加一個節(jié)點時，插值雖然易算，但若要增加一個節(jié)點時，全部基函數全部基函數 li(x) 都需重新算過。公式不具有繼承性，都需重新算過。公式不具有繼承性，不利于編程。不利于編程。將將 Ln(x) 改寫成改寫成.)()(102010 - - - - - xxxxaxxaa).(10- - - - nnxxx

24、xa的形式，希望每加一個節(jié)點時，的形式，希望每加一個節(jié)點時，只附加一項只附加一項上去即可。上去即可。? 差商差商( (亦稱均差亦稱均差) ) /* divided difference */),()()(,jijijijixxjixxxfxfxxf - - -= =1階差商階差商 /* the 1st divided difference of f w.r.t. xi and xj */)(,kixxxxfxxfxxxfkikjjikji - - -= =2階差商階差商f(x0)1010()()f xf xxx-2010201021()()()()f xf xf xf xxxxxxx-1階差商

25、的幾何階差商的幾何意義：弦截線意義：弦截線的斜率的斜率4 Newtons Interpolation4 Newtons Interpolation11101010111010,.,.,.,.,., - - - - - -= =- - -= =kkkkkkkkkkkxxxxxfxxxfxxxxxfxxxfxxf(k+1)階差商階差商： = = = =kiikikxxfxxf010)()(,., 事實上事實上其中其中,)()(01= = - -= =kiikxxx = = - -= = kijjjiikxxx01)()( Warning: my head is explodingWhat is t

26、he point of this formula?差商的值與差商的值與 xi 的順序無關！的順序無關！1.1.線性線性：1201020( )( )( ) , , ,kkkf xaf xbf xf xxaf xxbf xx=2.2.差商差商可以表示為可以表示為函數值的線性組合函數值的線性組合：0000-111( )( ),.,( )kkiikiiiiiiiikkif xf xf xxxxxxxxxxx=-,)()(01= = - -= =kiikxxx = = - -= = kijjjiikxxx01)()( 3.3. 對稱性：對稱性：由由2 2知，知，差商的值與節(jié)點的順序無關！差商的值與節(jié)點的

27、順序無關！4.4. 差商的另一種定義：由差商的另一種定義：由2，3及均差定義可得及均差定義可得10100 ,.,.,.,kkkkf xxf xxf xxxx-=-4 Newtons Interpolation6.5.5. 差商與導數的關系：差商與導數的關系：f(x)C ka,b,則則 a,b, s.t. 0,.,!kkff xxk=HW: 證明差商的性質證明差商的性質2，44 Newtons Interpolation4 Newtons Interpolation 牛頓插值牛頓插值 /* Newtons Interpolation */,)()()(000 xxfxxxfxf- - = =,)

28、(,101100 xxxfxxxxfxxf- - = =,.,)(,.,.,0010nnnnxxxfxxxxfxxxf- - = =- -).(.)()()(10102010- - - - - - - - - = =nnnxxxxaxxxxaxxaaxN .)(,)(,)()(102100100 - - - - - = =xxxxxxxfxxxxfxfxf).(,.,100- - - - nnxxxxxxf)().(,.,100nnnxxxxxxxxxf- - - - - -Nn(x)n次多次多項式，滿足：項式，滿足： Nn(xi)= f(xi)Rn(x)插值余項，插值余項，滿足滿足Rn(xi

29、)=0，i=0,n ai = f x0, , xi (1)(2)(n)(n)(n-1) (2) (1)4 Newtons Interpolation注：注： 由由唯一性可知唯一性可知 Nn(x) Ln(x)， 只是算法不同，表達只是算法不同，表達形式不同，故其余項也相同，即形式不同，故其余項也相同，即(1)011() ,.,( )( )(1)!nxnnnff x xxxxn=),(,!)(,.,maxmin)(0 xxkfxxfkk = = 實際計算過程為實際計算過程為f (x0)f (x1)f (x2)f (xn- -1)f (xn)f x0, x1f x1, x2 f xn- -1, xn

30、f x0, x1 , x2 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+1增加增加如果精度不構，增加節(jié)點如果精度不構，增加節(jié)點xn+1, ,同時表中增加一行，三角形斜邊同時表中增加一行，三角形斜邊上即為所要求的各次項系數。上即為所要求的各次項系數。4 Newtons Interpolation 等距節(jié)點公式等距節(jié)點公式 /* Formulae with Equal Spacing */向前差分向前差分 /* forward difference */i

31、iifff- -= = 1ikikikikffff1111)(- - - - - - - = = = = 向后差分向后差分 /* backward difference */111- - - - - - = = ikikikfffi- -1iifff- -= = 中心差分中心差分 /* centered difference */212111- - - - - -= =ikikikfff 其中其中)(221hiixff = = 當節(jié)點當節(jié)點等距等距分布時分布時:),.,0(0nihixxi= = = =fi= f(xi) 差分計算可通過構造差分表得到差分計算可通過構造差分表得到 2342345

32、0000000234111111232222233 = kkkkkkkx f xfffffxffffffxfffffxffffxf23344455 ffxffxf增加增加 23450011122222333 = kkkkkkkkxf xffffffxfxffxfffxff233323444444423455555555 ffxfffffxffffff1112211112211122identity calculus translation calculus: -, kkkkkkkkkkIffEffE ffE ffEffE II EEE-=-：， 差分的重要性質：差分的重要性質： 線性：例如線性

33、：例如gbfaxgbxfa = = )()( 各階差分可用函數值表示：各階差分可用函數值表示：!)1).(1(jjnnnjn - - -= = 其中其中/* binomial coefficients */4 Newtons Interpolation0011nnnjjnnjkkkn kjjjnnfEIfEffjj- -=-=-=-1()0011nnnnjnjnnjkkkkj njjnnfIEfEffjj- -=-=-=- 函數值可用各階差分表示：函數值可用各階差分表示：0e.g. Ennnjn kkkkjnffIffj= = 差商與差分的關系：差商與差分的關系：111121211222211

34、-1,-,-,1,-22, 1,!11,!kkkkkkkkkkkkkkkkkkkkkmkkkmkmmmkkkmkmkmmffff xxfxxxhf xxf xxfff xxxfxxhhgenerally we havef xxxfm hf xxxffm hm h-=- = 若若 f (x)是是 n 次多項式，則次多項式，則 是是 次多次多項式，而項式，而 ( ) (0)mf xmnnm-( )0 ()mf xmn= 差分與導數的關系差分與導數的關系( (由差分與差商、差商與導數的關系得）：由差分與差商、差商與導數的關系得）： ()1()1,!mmkkkmkmmmkmff xxxfmm hffh

35、=4 Newtons Interpolation等距節(jié)點牛頓公式等距節(jié)點牛頓公式 牛頓前差公式牛頓前差公式 /* Newtons forward-difference formula */001100010010-12000001( )( )( ),01, , ,( )( -)1( )(),( -),( -)( -)1111( -1)2!( ),nnjjnnnjjnnnnnf xNxRxxxthtxxjhxxtj hxx xht ttnNxf xf xxx xf xxxx xx xfftft tft ttnnRxf x xxx=-=-=-= -=01(1)( -)( -)1( )(1)!nnn

36、nx xx xt ttnhfn-=注：注：一般當一般當 x 靠近靠近 x0 時用前插，靠近時用前插，靠近 xn 時用后插，故兩時用后插，故兩種公式亦稱為種公式亦稱為表初公式表初公式和和表末公式表末公式。 牛頓后差公式牛頓后差公式 /* Newtons forward-difference formula */nn110nnn 1nnn 10n12nnnnnnn 10, 10, , , ( )( -)1( )(),( -),( -)( -)1111( -1)2!( ),( -jjnnnjjnnxxthtxxjhxxtj hxx xht ttnNxf xf x xx xf x xxx xx xff

37、tft tft ttnnR xf x x xxx x=-=- =-= =節(jié)點倒排，n01(1)( -)1( )(1)!nnx xt ttnhfn=HW: p. 59#9-16 5 厄米插值厄米插值 /* Hermite Interpolation */不僅要求函數值重合，而且要求若干階不僅要求函數值重合，而且要求若干階導數導數也重合。也重合。即：要求插值函數即：要求插值函數 (x) 滿足滿足 (xi) = f (xi), (xi) = f (xi), (mi) (xi) = f (mi) (xi).注：注： N 個條件可以確定個條件可以確定 階多項式。階多項式。N - - 1要求在要求在1個節(jié)

38、點個節(jié)點 x0 處直到處直到m0 階導數都重合的插階導數都重合的插值多項式即為值多項式即為Taylor多項式多項式00)(!)(.)()()(000)(000mmxxmxfxxxfxfx- - - - = = 其余項為其余項為)1(00)1(00)()!1()()()()(-=-=mmxxmfxxfxR 一般只考慮一般只考慮 f 與與f 的值。的值。厄米插值厄米插值3 Hermite Interpolation例：例：設設 x0 x1 x2, 已知已知 f(x0)、 f(x1)、 f(x2) 和和 f (x1), 求多項式求多項式 P(x) 滿足滿足 P(xi) = f (xi)，i = 0,

39、 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- - - - -= =h2(x)h1(x)有根有根 x0, x2 )()()(201xxxxBAxxh- - - = =由余

40、下條件由余下條件 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)()()4(xfxK = =與與 Lagrange 分析完全類分析完全類似似3 Hermite Interpola

41、tion一般地，已知一般地，已知 x0 , , xn 處有處有 y0 , , yn 和和 y0 , , yn ，求，求 H2n+1(x) 滿足滿足 H2n+1(xi) = yi ， H2n+1(xi) = yi。解：解：設設=ni)()()(=0iix x yixH2n+1n=0iyi其中其中 i(xj) = ij , i(xj) = 0, (xj) = 0, (xj) = ij i i(x)有根有根 x0 , , xi , , xn且都是且都是2重根重根 - - -= =ijjijixxxxxl)()()(由余下條件由余下條件 i (xi) = 1 和和 i(xi) = 0 可解可解Ai 和

42、和 Bi 2( )1 2 ( )()( )iiiiixl xxxlx=- (x) i有根有根 x0 , , xn, 除了除了xi 外都是外都是2重根重根 i)()(iili2(x)xxCx- -= =又又 i (xi) = 1 Ci = 1 i)(x)(ili2(x)xx- -= =設設,.210baCfbxxxann = = = =則則20)22()()!22()()( - - = = = niixnnxxnfxR 這樣的這樣的Hermite 插值唯插值唯一一 i)()()(2xlBxAxiii = = i類似的，類似的，Th. Th. 設設f(x)C 2n+2a,b,則則 a,b, s.t

43、. 滿足下面插值條件滿足下面插值條件3 Hermite InterpolationQuiz: 給定給定 xi = i +1, i = 0, 1, 2, 3, 4, 5. 下面哪個是下面哪個是 i(x)的圖像？的圖像？x0-10.5123456yxy0-10.5123456斜率斜率=1 求求Hermite多項式的基本步驟：多項式的基本步驟： 寫出相應于條件的寫出相應于條件的 i(x)、 i(x) 的組合形式；的組合形式； 對每一個對每一個 i(x)、 i(x)找出盡可能多的條件給出的根；找出盡可能多的條件給出的根； 根據多項式的總階數和根的個數寫出表達式；根據多項式的總階數和根的個數寫出表達式；

44、 根據尚未利用的條件解出表達式中的待定系數；根據尚未利用的條件解出表達式中的待定系數； 最后完整寫出最后完整寫出H2n+1(x)。HW: p.59#17-19注：待定系數法仍適用，但插值節(jié)點多注：待定系數法仍適用，但插值節(jié)點多時比較麻煩。時比較麻煩。4 分段低次插值分段低次插值 /* piecewise polynomial approximation */Remember what I have said? Increasing the degree of interpolating polynomial will NOT guarantee a good result, since hig

45、h-degree polynomials 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) 分段分段低次低次插值插值4 Piecewise Polynomial Approximation 分段線性插值分段線性插值 /* piecewise linear interpolatio

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

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

48、= =.10,)(2baCxS ,1 iixx),., 0(),()(nixfxSii= = =注：注：三次樣條與分段三次樣條與分段 Hermite 插值的根本區(qū)別在于插值的根本區(qū)別在于S(x)自自身光滑身光滑，不需要知道，不需要知道 f 的導數值（除了在的導數值（除了在2個端點可能需個端點可能需要）；而要）；而Hermite插值依賴于插值依賴于f 在所有插值點的導數值。在所有插值點的導數值。f(x)H(x)S(x)5 Cubic Spline 構造三次樣條插值函數的構造三次樣條插值函數的三彎矩法三彎矩法 /* method of bending moment */在在 上，記上，記,1jjx

49、x- -,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)

50、(21121Sj(x) =jjjjjjjjBxAhxxMhxxM - - - - - -6)(6)(3131Sj(x) =利用已知利用已知Sj(xj- -1) = yj- -1 Sj(xj) = yj 可解可解5 Cubic SplinejjjjjjjhMMhyyA611- - - - - -= =jjjjjjjjjjjjhxxhMyhxxhMyBxA12211)6()6(- - - - - - - - -= = 下面解決下面解決 Mj ： 利用利用S 在在 xj 的的連續(xù)性連續(xù)性xj- -1, xj : Sj(x) =jjjjjjjjjjjhMMxxfhxxMhxxM6,2)(2)(1121

51、21- - - - - - - - - - - -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 = =

52、 - - - -110111122nnnnggMMl lm ml lm m還需還需2個個邊界條件邊界條件 /* boundary conditions */5 Cubic Spline 第第1類邊條件類邊條件 /* clamped boundary */： S(a) = y0， S(b) = yna , x1 : S1(x) =1011012112106,2)(2)(hMMxxfhaxMhxxM- - - - - - - -010110),(62gy0 xxfhMM= =- -= = nnnnnngxxfynhMM= =- -= = - - -),(6211類似地利用類似地利用 xn- -1,

53、 b 上的上的 Sn(x) 第第2類邊條件：類邊條件： S”(a) = y0” = M0， S”(b) = yn” = Mn這時：這時：nnnygyg = = = = = =2,0;2,0000m ml l特別地，特別地，M0 = Mn = 0 稱為稱為自由邊界自由邊界 /* free boundary */，對應的對應的樣條函數稱為樣條函數稱為自然樣條自然樣條 /* Natural Spline */。 第第3類邊條件類邊條件 /* periodic boundary */ ： 當當 f 為為周期函數周期函數時，時， yn = y0 ， S(a+) = S(b- -) M0 = Mn = =

54、 - - -nnnnnnggMM111122112222m ml ll lm ml lm mm ml l5 Cubic Spline注：注：另有另有三轉角法三轉角法(p.49-53)得到樣條函數，即設得到樣條函數，即設 Sj(xj) = mj，則易知，則易知xj- -1, xj 上的上的Sj(x) 就是就是Hermite函數函數。再。再利用利用S”的連續(xù)性，可導出關于的連續(xù)性，可導出關于mj 的方程組，加上邊的方程組，加上邊界條件即可解。界條件即可解。Cubic Spline 由由boundary conditions 唯一唯一確定。確定。收斂性：收斂性：若若 ，且，且 ，則，則,baCf Chhiiminmax一致一致S(x) f(x)0maxihas即即:提高精度只須提高精度只須增加節(jié)點增加節(jié)點, 而無須提高樣條階數。而無須提高樣條階數。穩(wěn)定性：穩(wěn)定性：只要邊條件保證只要邊條件保證 | m m0 |, | l l0 |, | m mn |, | l l n | 2，則方程組系數陣為則方程組系數陣為SDD陣陣，保證數值穩(wěn)定。保證數值穩(wěn)定。HW: p.59-60,#19-255 Cubic Spline Sketch of the Algorithm: Cubic Spline 計算計算 m mj , l l j , gj ; 計算計算 Mj (