1、一、數值求積的基本思想一、數值求積的基本思想)()()(aFbFdxxfba 積分積分 只要找到被積函數只要找到被積函數 f (x)原函數原函數F(x)，便有，便有牛頓牛頓萊布尼茲萊布尼茲(NewtonLeibniz)公式公式 baxxfId)(實際困難實際困難：大量的被積函數（：大量的被積函數（ , sin x2 等）等）, 找不到用初等函找不到用初等函數表示的原函數數表示的原函數；另外；另外, f (x)是（測量或數值計算出的）一張數是（測量或數值計算出的）一張數據表時，據表時，牛頓牛頓萊布尼茲公式萊布尼茲公式也也不能直接運用不能直接運用。xxsin 積分中值定理：在積分中值定理：在a,

2、b內存在一點內存在一點 ，有，有 f( )成立。成立。 )(d)(abxxfba 1 引言引言第四章第四章 數值積分與數值微分數值積分與數值微分 就是說就是說, 底為底為b- -a 而高為而高為f( )的的矩形面積矩形面積恰恰等于所求等于所求曲邊梯形的面積曲邊梯形的面積 .問題問題 在于點在于點的具體位置一般是不知道的，因而的具體位置一般是不知道的，因而 難以準確算出難以準確算出 f( )的值的值我們將我們將f ( )稱為區間稱為區間a, b上的平均高度這樣上的平均高度這樣,只要對只要對平均高度平均高度f( )提供一種算法提供一種算法,相應地便獲得一種數值求積方法相應地便獲得一種數值求積方法

3、如果用兩端點的如果用兩端點的“高度高度”f(a)與與f(b)的算術平均作為平均高度的算術平均作為平均高度f ( ) 的近似值，這樣導出的求積公式的近似值，這樣導出的求積公式 : 便是我們所熟悉的便是我們所熟悉的梯形公式梯形公式 . )()(2bfafabT 2)(bafabR2bac 而如果改用區間中點而如果改用區間中點 的的“高度高度”f (c)近似地取代平近似地取代平均高度均高度f ( )，則又可導出所謂，則又可導出所謂中矩形公式中矩形公式(今后簡稱矩形公式今后簡稱矩形公式)：(1.1)(1.2) 更一般地，我們可以在區間更一般地，我們可以在區間a，b上適當選取某些節點上適當選取某些節點

4、xk ，然后然后用用 f (xk )加權平均得到平均高度加權平均得到平均高度 f ()的近似值的近似值，這樣構造出的，這樣構造出的求積公式具有下列形式求積公式具有下列形式式中式中 xk 稱為稱為求積節點求積節點；Ak 稱為稱為求積系數求積系數，亦稱為伴隨節點，亦稱為伴隨節點 xk 的的權權權權Ak 僅僅與節點僅僅與節點xk 的選取有關，而不依賴于被積函數的選取有關，而不依賴于被積函數 f(x)的的具體形式具體形式 ban0kkkxfAdxxf)()(使積分公式具有通用性使積分公式具有通用性 這類數值積分方法通常稱作能這類數值積分方法通常稱作能機械求積機械求積， 其特點是將積分其特點是將積分求值

5、問題歸結為函數值的計算，這就避開了牛頓求值問題歸結為函數值的計算，這就避開了牛頓萊布尼茲公萊布尼茲公式需要尋求原函數的困難式需要尋求原函數的困難(1.3)二、二、代數精度的概念代數精度的概念 數值求積方法是近似方法，為要保證精度，我們自然希望求積數值求積方法是近似方法，為要保證精度，我們自然希望求積公式能對公式能對“盡可能多盡可能多”的函數準確地成立，這就提出了所謂代數精的函數準確地成立，這就提出了所謂代數精度的概念度的概念 定義定義 1 如果某個求積公式對于次數如果某個求積公式對于次數m的多項式均能準確地成的多項式均能準確地成立，但對于立，但對于m+1次多項式就不一定準確，則稱該求積公式具有

6、次多項式就不一定準確，則稱該求積公式具有m次代次代數精度數精度 一般地，欲使求積公式一般地，欲使求積公式 具有具有m次代數次代數精度，只要令它對于精度，只要令它對于f (x) = 1，x，xm 都能準確成立，這就要求都能準確成立，這就要求 bankkkxfAxxf0)(d)( . )(11;)(21;1122mmmkkkkkabmxAabxAabA例例1： 考察其代數精度。考察其代數精度。 f(x)abf(a)f(b)梯形公式梯形公式解：解：逐次檢查公式是否精確成立逐次檢查公式是否精確成立代入代入 P0 = 1： baabdx111 2 ab=代入代入 P1 = x ：=代入代入 P2 = x

7、2 ： 222abbadxx 2baab 3233abbadxx 222baab 代數精度代數精度 = 1)()(2)(bfafabdxxfba 例例2 試構造形如試構造形如 f(x)dx A0f(0)+ A1f(h)+ A2f(2h) 的數值的數值求積公式求積公式,使其代數精度盡可能高使其代數精度盡可能高,并指出其代數精度的階數并指出其代數精度的階數.3h0解解: 令公式對令公式對 f(x)=1,x, x2 均準確成立均準確成立,則有則有3h=A0+ A1+ A2h2=0 + A1h+ A22h9h3=0 + A1h2+ A24h229故求積公式的形式為故求積公式的形式為解之得解之得 A0=

8、 h, A1=0, A2= h. 94 34 f(x)dx f(0) + f(2h)3h43h43h0由公式的構造知由公式的構造知,公式公式至少至少具有具有2次代數精度次代數精度; 而當而當f(x)=x3時時,公式的左邊公式的左邊= h4, 右邊右邊=18h4, 公式的左邊公式的左邊 右邊右邊,說明說明此公式對此公式對 f(x)=x3不能準確成立不能準確成立.因此因此,公式只具有公式只具有2次代數次代數精度精度.814三、三、求積公式的收斂性與穩定性求積公式的收斂性與穩定性 定理定理3表明，只要求積系數表明，只要求積系數Ak0 (k0,1,n)，就能保證，就能保證計算的穩定性計算的穩定性 定義

9、定義2 在求積公式在求積公式 中，若中，若 其中其中 ，則稱求積公式是收斂的，則稱求積公式是收斂的 由于計算由于計算 f (xk)可能有誤差可能有誤差,實際得到實際得到 定義定義3 對任給對任給 e e 0，若，若 (k=0,1, ,n), 就有就有 , 則稱求積公式是穩定的則稱求積公式是穩定的. bankkkxfAxxf0)(d)(e e |)(|00nkkknkkkfAxfA)(max11 iinixxhkkfxf)(0 ，只只要要 .)(,kkkkfxff 即即 bankkkhnxxfxfAd)()(lim00 定理定理3 若求積公式若求積公式(13)中系數中系數Ak0 (k0，1，n)

10、，則此求積公式是穩定的則此求積公式是穩定的近似近似計算計算 badxxfI)(思思路路利用利用插值多項式插值多項式 則積分易算。則積分易算。)()(xfxPn 在在a, b上取上取 a x0 x1 xn b，做，做 f 的的 n 次插值次插值多項式多項式 ，即得到，即得到 nkkknxlxfxL0)()()( babaknkkdxxlxfdxxf)()()(0Ak bakjxxxxkdxAjkj)()(由由 決定，決定，與與 無關。無關。節點節點 f (x)插值型積分公式插值型積分公式bankkxnbanbanbankkkdxxxnfdxxRdxxLxfxfAdxxffR0)1(0)()!1(

11、)()()()()()(誤差誤差bandxxP)(2 插值型的求積公式插值型的求積公式與與Newton-Cotes 公式公式關鍵是關鍵是f(x)如果求積公式是插值型的如果求積公式是插值型的, 按余項式按余項式, 對于次數對于次數 n的多項式的多項式 f (x)，其余項其余項R f 等于等于0，因而這時求積公式至少具有，因而這時求積公式至少具有n次代數精度次代數精度定理定理1：形如形如 的求積公式至少有的求積公式至少有 n 次代數精度次代數精度 該該公式為公式為插值型插值型（即：（即： ） nkkkxfA0)( bakkdxxlA)(為便于計算，一般取為便于計算，一般取等距離節點等距離節點得到近

12、似公式得到近似公式:一、一、Newton-Cotes 公式公式2、 把把a, b二等分，作二等分，作2次插值，有次插值，有)()(4)( )(66bffafdxxfbaabba此公式稱為此公式稱為辛普森（辛普森（Simpson）公式）公式。badxxL)(21、 對于對于a, b上上1次插值，有次插值，有)()()(1bfafxLabaxbabx )()()(2221bfafdxxfAAabbaab 此即此即梯形公式梯形公式。 節點節點等距分布等距分布：ninabhhiaxi,., 1, 0, dxxxxxAnxxijjiji 0)()( njiinnjidtjtininabdthhjihjt

13、00)()!( !)1)()()(令令htax Cotes系數系數)(niC注：注：Cotes 系數僅取決于系數僅取決于 n 和和 i，可查表得到。與可查表得到。與 f (x) 及區及區間間a, b均無關。均無關。 3、 把把a, b n 等分，用插值等分，用插值Ln(x)近似近似 f(x)積分，有積分，有當當n=4時時, 牛頓牛頓-柯特斯公式特別稱作柯特斯公式特別稱作柯特斯公式柯特斯公式,其形式為其形式為 )(7)(32)(12)(32)(79043210 xfxfxfxfxfabC 21,21)1(1)1(0 CCn = 1:)()(2)(bfafabdxxfba Trapezoidal

14、RuledxbxaxffRbax)(!2)( /* 令令 x = a+th, h = b a, 用中用中值定理值定理 */1, , )(1213abhbafh 代數精度代數精度 = 1n = 2:61,32,61)2(2)2(1)2(0 CCC)()(4)(6)(2bffafabdxxfbaba Simpsons Rule代數精度代數精度 = 32,),(,)(901)4(5abhbafhfR n = 3: Simpsons 3/8-Rule, 代數精度代數精度 = 3,)(803)5(5 fhfR n = 4: Cotes Rule, 代數精度代數精度 = 5,)(9458)6(7 fhfR

15、 n 為為偶數階偶數階的的Newton-Cotes 公式至少有公式至少有 n+1 次代數精度。次代數精度。二、幾種低階求積公式的余項二、幾種低階求積公式的余項三、三、偶階求積公式的代數精度偶階求積公式的代數精度 作為插值型的求積公式，作為插值型的求積公式，n 階的牛頓階的牛頓-柯特斯公式至柯特斯公式至少具有少具有n 次的插值精度（定理次的插值精度（定理1）。實際的代數精度還可）。實際的代數精度還可進一步提高，一般地，可以證明下述定理進一步提高，一般地，可以證明下述定理: 定理定理 2 當階當階 n 為偶數時，牛頓為偶數時，牛頓-柯特斯公式柯特斯公式至少有至少有 n+1 次代數精度次代數精度 .

16、 nkknknxfCabI0)()()(注：注：由公式知，當由公式知，當n8時，柯特斯系數出現負值，這時時，柯特斯系數出現負值，這時，初始數據誤差將會引起計算結果誤差增大，即計算不，初始數據誤差將會引起計算結果誤差增大，即計算不穩定。因此，實際計算不用穩定。因此，實際計算不用n8的牛頓的牛頓-柯特斯公式柯特斯公式 .估計估計截斷誤差截斷誤差為為解解 用用梯形公式梯形公式計算計算:=2.1835估計估計截斷誤差截斷誤差為為=0.6796用用Simpson公式公式計算：計算：=2. 0263例例3 試分別使用梯形公式和試分別使用梯形公式和Simpson公式計算積分公式計算積分 的近的近似值，并估計

17、截斷誤差似值，并估計截斷誤差.=198.4306890. 0)(max2880)12()4(2152 xfRx=0.068903 3 復化求積公式復化求積公式高次插值有高次插值有Runge 現象現象，故采用分段低次插值，故采用分段低次插值 分段低次合成的分段低次合成的 Newton-Cotes 復合復合求積公式。求積公式。一、復化梯形公式一、復化梯形公式:),., 0(,nkhkaxnabhk 在每個在每個 上用梯形公式：上用梯形公式：,1kkxx 11)()(2)(2nkkbfxfafh bankkkxfxfhdxxf11)()(2)(=Tn),(),()(12)()(12)(1221213

18、bafabhnfabhfhfRnkknkk /*中值定理中值定理*/nkxfxfxxdxxfkkxxkkkk,., 1,)()(2)(111 二、復化辛普森公式二、復化辛普森公式:),., 0(,nkhkaxnabhk )()(4)(6)(1211 kkkxxxfxfxfhdxxfkkkx21 kx1 kx44444 )()(2)(4)(6)(1010121 nknkkkbabfxfxfafhdxxf= Sn)(2180)4(4 fhabfR 注：注：為方便編程，可采用另一記法：令為方便編程，可采用另一記法：令 n = 2n 為偶數，為偶數， 這時這時 ，有，有hkaxhnabhk ,2 )(

19、)(2)(4)(3 koddkevenkknbfxfxfafhS三、收斂速度與誤差估計：三、收斂速度與誤差估計：定義定義 若一個積分公式的誤差滿足若一個積分公式的誤差滿足 且且C 0，則則稱該公式是稱該公式是 p 階收斂階收斂的。的。 ChfRphlim0)(,)(,)(642hOChOShOTnnn例例4：計算計算dxx 10142 解：解： )1()(2)0(161718fxffTkk8kxk 其中其中= 3.138988494 )1()(2)(4)0(241oddeven4fxfxffSkk8kxk 其中其中= 3.141592502運算量基運算量基本相同本相同問題問題: 給定精度給定精

20、度 e e，如何取，如何取 n ?例如：要求例如：要求 ，如何判斷，如何判斷 n = ？e e |nTI)()(122 fabhfR ? nkkhfh12)(12 )()(12)(1222afbfhdxxfhba 上述上述例例4中若要求中若要求 , 則則610| nTI622106| )0() 1 (|12| | hffhfRn00244949.0 h即：取即：取 n = 409通常采取將區間通常采取將區間不斷對分不斷對分的方法，即取的方法，即取 n = 2k上述上述例例4中中2k 409 k = 9 時，時，T512 = 3.14159202例例4中：中：S4 = 3.141592502注意

21、到區間再次對分時注意到區間再次對分時412)()(12122fRhafbffRnn 412 nnTITI)(3122nnnTTTI 可用來判斷迭代可用來判斷迭代是否停止。是否停止。(1)(2)(3)事后誤差估計事后誤差估計一、梯形法的遞推化一、梯形法的遞推化逐次分半法逐次分半法 上一節介紹的復化求積方法對提高精度是行之有效的，但上一節介紹的復化求積方法對提高精度是行之有效的，但在使用求積公式之前必須給出合適的步長，在使用求積公式之前必須給出合適的步長，步長步長取得取得太大精度太大精度難以保證難以保證，步長太小步長太小則會導致則會導致計算量計算量的的增加增加，而事先給出一個，而事先給出一個恰當的

22、步長又往往是困難的恰當的步長又往往是困難的 實際計算中常常實際計算中常常采用變步長的計算方案采用變步長的計算方案，即在步長，即在步長逐次分逐次分半半(即步長二分即步長二分)的過程中，反復利用復化求積公式進行計算，的過程中，反復利用復化求積公式進行計算，直至所求得的積分值滿足精度要求為止直至所求得的積分值滿足精度要求為止 設將求積區間設將求積區間a，b分成分成n等分，則一共有等分，則一共有n+1個分點，按個分點，按梯形公式計算積分值梯形公式計算積分值Tn，需要提供，需要提供n+1個函數值如果將求積個函數值如果將求積區間再二分一次，則分點增至區間再二分一次，則分點增至2n+1個，我們來個，我們來考

23、察考察二分二分前后兩前后兩個積分值個積分值之間的之間的聯系聯系4 4 龍貝格求積公式龍貝格求積公式逐次分半逐次分半計算計算方案方案的實現的實現: 注意到每個子區間注意到每個子區間xk，xk+1經過二分只增加了一個分經過二分只增加了一個分點點 xk+1/2( xk+xk+1)/2，用復化梯形公式求得該子區間上的，用復化梯形公式求得該子區間上的積分值為積分值為 101021102110122)12(221)(221)(2)()(4nknnkknnkknkkknhkafhTxfhTxfhxfxfhT)()(2)(4121 kkkxfxfxfh這里這里 代表二分前后的步長代表二分前后的步長. .將每個

24、子區間上的積分將每個子區間上的積分值相加得值相加得nabh 二、龍貝格算法二、龍貝格算法).,()(212);,()(12)(222bafhabTIbafhabTIfRnnn 有有：根據復化梯形公式的余項表達式根據復化梯形公式的余項表達式. )(31.41)()(222nnnnnTTTITITIff 整理后可得：整理后可得：，則有，則有假定假定 可見，可見，利用利用兩種步長兩種步長計算的結果能估計截斷誤差計算的結果能估計截斷誤差.若將該截斷若將該截斷誤差加到計算結果中誤差加到計算結果中，nnnnnTTTTTT3134)(31222 就得出就得出“改進的梯形求積公式改進的梯形求積公式”：事后誤差

25、事后誤差估計估計例：例：計算計算dxx 10142 已知對于已知對于e e = 10 6 須將區間對分須將區間對分 9 次，得到次，得到 T512 = 3.14159202由由 來計算來計算 I 效果是否好些？效果是否好些？nnnnTTTTI313414422 483134TT = 3.141592502 = S4改進梯形求積公式改進梯形求積公式的右邊實際是的右邊實際是nnknkkknkknkknkknnnkknnnSbfxfxfafhxfhbfxfafhxfhTTxfhTTT 101121102111102110212)()(2)(4)(6)(2)()(2)(231)(231)(221431

26、)4(31這就是說用這就是說用梯形法二分前后的兩個積分值梯形法二分前后的兩個積分值Tn與與T2n的的線性組合線性組合的結果的結果得到得到復化辛普森法求積公式復化辛普森法求積公式nnnnnTTTTS141144313422 類似的情況，用辛普森法二分前后的兩個積分值類似的情況，用辛普森法二分前后的兩個積分值Sn與與S2n的線性組合的結果可得到的線性組合的結果可得到復化柯特斯求積公式復化柯特斯求積公式nnnnnSSSSC151151614114422222 重復同樣的手續，用柯特斯法二分前后的兩個積分值重復同樣的手續，用柯特斯法二分前后的兩個積分值Cn與與C2n的線性組合的結果可得到的線性組合的結

27、果可得到龍貝格龍貝格(Romberg)求積公式求積公式nnnnnCCCCR631636414114423233 我們在變步長的過程中運用加速公式，就能將粗糙的梯我們在變步長的過程中運用加速公式，就能將粗糙的梯形值形值Tn逐步加工成精度較高的辛普森值逐步加工成精度較高的辛普森值Sn 、柯特斯值、柯特斯值Cn和龍和龍貝格值貝格值Rn .一般有：一般有：nnnSTT 1442nnnCSS 144222nnnRCC 144323 Romberg 算法：算法： e e ? e e ? e e ? T1 =)0(0T T8 =)3(0T T4 =)2(0T T2 =)1(0T S1 =)0(1T R1 =

28、)0(3T S2 =)1(1T C1 =)0(2T C2 =)1(2T S4 =)2(1TRomberg 序列序列kk2kT212 kS22 kC32 kR0 20=1 T11 21=2 T2 S12 22=4 T4 S2 C13 23=8 T8 S4 C2 R14 24=16 T16 S8 C4 R25 25=32 T32 S16 C8 R4 區間等分數區間等分數 梯形序列梯形序列 辛普森序列辛普森序列 柯特斯序列柯特斯序列 龍貝格序列龍貝格序列 龍貝格求積算法可用下表來表示：龍貝格求積算法可用下表來表示： 例例5 用龍貝格方法計算橢圓用龍貝格方法計算橢圓 x2/4 + y2 l 的周長，使

29、結果的周長，使結果具有五位有效數字具有五位有效數字 分析分析 為便于計算，先將橢圓方程采用參數形式表示為便于計算，先將橢圓方程采用參數形式表示, ,再根再根據弧長公式將橢圓周長用積分形式表示由于計算結果要求具據弧長公式將橢圓周長用積分形式表示由于計算結果要求具有五位有效數字，因此需要估計所求積分值有幾位整數，從而有五位有效數字，因此需要估計所求積分值有幾位整數，從而確定所求積分值的絕對誤差限最后再應用龍貝格方法計算積確定所求積分值的絕對誤差限最后再應用龍貝格方法計算積分分 解解 令令 x 2cosq q，y sinq q ， 則橢圓的周長為則橢圓的周長為Iyxl4d sin314d420220

30、22 q qq qq qq qq q.10125. 01081)(1021)(4422d sin3124451202 fRIfRIlI的的截截斷斷誤誤差差為為故故計計算算，的的截截斷斷誤誤差差為為則則需需結結果果有有五五位位有有效效數數字字，有有一一位位整整數數，要要求求，因因此此由由于于 q qq q 下表給出了用龍貝格方法計算積分下表給出了用龍貝格方法計算積分I= 1+1+3sin2q q dx 的過程的過程. /20kk2kT212 kS22 kC32 kR4322 kkRR0 1 2.356 1941 2 2.419 921 2.441 1632 4 2.422 103 2.422 8

31、30 2.421 608 3 8 2.422 112 2.422 115 2.422 067 2.422 074 4 16 2.422 112 2.422 112 2.422 112 2.422 113 0.000 0395 32 2.422 112 2.422 112 2.422 112 2.422 112 0.000 001 0.125 10- -4 故積分故積分I 2.422112, 橢圓周長的近似值為橢圓周長的近似值為l = 4I 9.6884。三、理查森三、理查森(Richardson)外推加速法外推加速法 上面討論說明由梯形公式出發上面討論說明由梯形公式出發, 將區間將區間a, b

32、逐次二分逐次二分可提高求積公式的精度可提高求積公式的精度, 上述加速過程還可繼續下去上述加速過程還可繼續下去. 下面我們討論其下面我們討論其理論依據理論依據. ,)(24221 llhhhIhT .)(122nabhbafhabTIn ， 22hTTn若記若記Tn = T(h), 當區間當區間a, b分為分為2n等分時等分時, 有有 , 則則可見可見I = T(h)的誤差為的誤差為O(h2). llhhhIhT2422121642 3)(24)(1hThThT 若記若記 ，則，則 將梯形公式按余項展開將梯形公式按余項展開. 由誤差公式有由誤差公式有 62411)(hhIhT 641626241

33、1hhIhT 顯然顯然T1(h)與與 I 近似的階為近似的階為O(h4) . 就是就是辛普森公式辛普森公式序列序列Sn, S2n, . .， 2),(11hThT這樣構造的這樣構造的 )(1412144)(11hThThTmmmmmm 則又可進一步從余項中則又可進一步從余項中消去消去 h4 項，這樣構造出的項，這樣構造出的 ，其實就是，其實就是柯特斯公式柯特斯公式序序列，它與列，它與 I 的逼近階為的逼近階為O(h6) . )(2hT)(151 21516)(112hThThT 若令若令 , 一般地，若記一般地，若記T0(h) = T(h)，經過，經過m (m =1,2,)次加速次加速后，則有

34、后，則有如此繼續下去，每加速一次，誤差的量級便提高如此繼續下去，每加速一次，誤差的量級便提高2階階. )21(141144)(1)1(1)()(0)()(0，次次加加速速值值，可可得得的的序序列列表表示示以以次次后后求求得得的的梯梯形形值值，且且表表示示二二分分設設以以 kTTTmTTkTkmmkmmmkmkkmk. ., 321.數數表表來來計計算算構構造造一一個個三三角角形形數數表表根根據據公公式式可可以以逐逐行行龍龍貝貝格格序序列列公公式式辛辛普普森森、柯柯特特斯斯、即即可可得得到到加加速速、若若取取算算法法上上式式也也稱稱為為龍龍貝貝格格求求積積Tm Romberg 算法算法 可以證明

35、，如果可以證明，如果 f (x) 充分光滑，那么充分光滑，那么T 數表每一列的數表每一列的元素及對角線元素均收斂到所求的積分值元素及對角線元素均收斂到所求的積分值 I ，即，即ITmITkmmkmk )()(lim)(lim，固固定定 理查德森理查德森外推法外推法利用利用低低階公式產生階公式產生高高精度的結果。精度的結果。設對于某一設對于某一 h 0，有公式，有公式 T0(h) 近似計算某一未知值近似計算某一未知值 I。由。由Taylor展開得到：展開得到： T0(h) I = 1 h + 2 h2 + 3 h3 + i 與與 h 無關無關現將現將 h 對分，得：對分，得：( () )( ()

36、 )( () ).)(3232222120 hhhhIT Q：如何將公式精度由如何將公式精度由 O(h) 提高到提高到 O(h2) ？.432112)()(23322020 hhIhTTh 即：即：.12)()(2)(32210201 hhIhTThTh .)(42312 hhIhT 12)()(221212 hTTh.)(2211 mmmhhIhT 12)()(2121 mmhmmhTTHW: p.159 6,7,8.3) 在構造在構造Newton-Cotes公式公式時，限定用積分區間的時，限定用積分區間的等分點等分點作為求積節點作為求積節點，這樣做雖然使問題的處理過程得以簡化，但，這樣做雖

37、然使問題的處理過程得以簡化，但同時也同時也限制了精度限制了精度。 求積公式含有求積公式含有2n+2個待定參數個待定參數xk、Ak(k0,1,n)若用若用待定系數法確定它們待定系數法確定它們, 則最好需要則最好需要2n+2個獨立的條件聯立方個獨立的條件聯立方程組求解程組求解, 從而易知求積公式的從而易知求積公式的最大代數精度最大代數精度可達到可達到2n+1次次. 在節點數目固定為在節點數目固定為n 的條件下，能否通過的條件下，能否通過適當選取求積適當選取求積節點節點xk的位置以及相應的求積系數的位置以及相應的求積系數Ak，使求積公式，使求積公式具有盡可能高具有盡可能高(最高最高)的代數精度？的代

38、數精度？ bankkkxfAxxf0)(d)(這類求積公式稱為這類求積公式稱為高斯高斯(Gauss)求積公式求積公式。4 高斯求積公式高斯求積公式 將節點將節點 x0 xn 以及系數以及系數 A0 An 都作為待定系數。都作為待定系數。令令 f (x) = 1, x, x2, , x2n+1 代入可求解，得到的公式代入可求解，得到的公式具有具有2n+1 次代數精度。這樣的節點稱為次代數精度。這樣的節點稱為Gauss 點點，公式稱為公式稱為Gauss 型求積公式型求積公式。 baxxxfId)()( 為使問題更具一般性為使問題更具一般性,我們研究帶權積分我們研究帶權積分 (x)為權函數為權函數,

39、 Ak(k0,1,n)為不依賴于為不依賴于f (x)的求積系數的求積系數, xk (k0,1,n)為求積節點為求積節點. bamnkmkknmxxxxA)2.5(.12,1,0d)(0 要使要使(5.1)具有具有2n+1次代數精度，則需要滿足次代數精度，則需要滿足 bankkkxfAdxxfx0)()()( 構造具有構造具有2n+1次代數精度的求積公式次代數精度的求積公式(5.1) 從例中可看到求解非線性方程組從例中可看到求解非線性方程組(5.2)較復雜，通常較復雜，通常n2就很難求解故一般不通過解方程就很難求解故一般不通過解方程(5.2)求求 xk 及及 Ak (k0,1, , n)例：例：

40、求求 的的 2 點點 Gauss 公式。公式。dxxfx)(10 解：解：設設 ，應有，應有 3 次代數精度。次代數精度。 101100)()()(xfAxfAdxxfx代入代入 f (x) = 1, x, x2, x3 31130092211200721100521032xAxAxAxAxAxAAA2776. 03891. 02899. 08212. 01010 AAxx不是線性方程組，不是線性方程組，不易求解。不易求解。 而從研究而從研究高斯點的基本特性高斯點的基本特性來著手解決來著手解決Gauss 求積公式求積公式的構造問題的構造問題由插值型公式構由插值型公式構造知造知,關鍵求關鍵求xk

41、，一、高斯點的基本特性一、高斯點的基本特性0)(d)()()(1 banxxxxxP 證明：證明： “” x0 xn 為為 Gauss 點點, 則公式則公式 至少有至少有 2n+1 次代數精度。次代數精度。 bankkkxfAdxxfx0)()()( 對任意次數對任意次數不大于不大于n 的多項式的多項式 Pm(x)， Pm(x) w(x)的次數的次數不大于不大于2n+1，則代入公式應則代入公式應精確成立精確成立： nkkkmkbamxwxPAdxxwxPx0)()()()()( 0= 0 “” 要證明要證明 x0 xn 為為 Gauss 點，即要證公式對任意次點，即要證公式對任意次數數不大于不

42、大于2n+1 的多項式的多項式 Pm(x) 精確成立，即證明：精確成立，即證明： nkkmkbamxPAdxxPx0)()()( 設設)()()()(xrxqxwxPm bababamdxxrxdxxqxwxdxxPx)()()()()()()( 0 nkkkxrA0)( nkkmkxPA0)( x0 xn 為為 Gauss 點點 與任意次數與任意次數不大于不大于n 的多項式的多項式 P(x) （帶權）正交（帶權）正交。 nkkxxxw0)()(定理定理求求 Gauss 點點 求求w(x)的的零點零點 Gauss 公式的余項：公式的余項： bankkkxfAdxxffR0)()(/* 設設P為

43、為f 的過的過x0 xn的插值多項式的插值多項式 */ bankkkxPAdxxf0)()(/*只要只要P 的階數不大于的階數不大于2n+1，則下一步，則下一步等式成立等式成立*/dxxPxfdxxPdxxfbababa)()()()( 插值多項式插值多項式的余項的余項Q：什么樣的什么樣的插值多項式插值多項式在在 x0 xn 上有上有 2n+1 階？階？A：Hermite 多項式！多項式！ 滿足滿足)()(),()(kkkkxfxHxfxH badxxHxffR)()(),(,)()!22()()()!22()(2)12(2)12(badxxwnfdxxwnfbanbaxn 二、高斯求積公式的

44、余項二、高斯求積公式的余項三、高斯求積公式的穩定性與收斂性三、高斯求積公式的穩定性與收斂性 定理定理6 高斯求積公式高斯求積公式(5.1)的求積系數的求積系數 Ak (k0,1,n)全是正的全是正的 由本定理及定理由本定理及定理2，則得，則得 推論推論 高斯求積公式高斯求積公式(5.1)是穩定的是穩定的. 定理定理7 設設 f (x)C a，b，則高斯求積公式，則高斯求積公式(5.1)是收斂是收斂 的，即的，即 nkbakknxxxfxfA0.d)()()(lim 正交多項式族正交多項式族 0, 1, , n, 有性質：任意次數不大有性質：任意次數不大于于n 的多項式的多項式 P(x) 必與必

45、與 n+1 正交。正交。若取若取 w(x) 為其中的為其中的 n+1，則，則 n+1的根的根就是就是 Gauss 點。點。再解上例：再解上例： 101100)()()(xfAxfAdxxfxStep 1：構造正交多項式構造正交多項式 2設設cbxxxaxxx 2210)(,)(, 1)( 53 a0)(10 dxaxx0),(10 1021102100)(53(0),(0)(0),(dxcbxxxxdxcbxxx 215910 cb即：即：215910)(22 xxx 四、常用的高斯型求積公式四、常用的高斯型求積公式Step 2：求求 2 = 0 的的 2 個根，即為個根，即為 Gauss 點

46、點 x0 ，x1221/20)9/10(9/1021;0 xStep 3：代入代入 f (x) = 1, x 以求解以求解 A0 ，A1解解線性線性方程組，方程組，簡單。簡單。結果與前一方法相同：結果與前一方法相同：2776. 0,3891. 0,2899. 0,8212. 01010 AAxx 利用此公式計算利用此公式計算 的值的值 10dxexx2555. 1 10dxexx2899. 08212. 0102776. 03891. 010eeeAeAxx 注：注：構造正交多項式也可以利用最小二乘數據擬合中介紹構造正交多項式也可以利用最小二乘數據擬合中介紹過的遞推式進行。過的遞推式進行。 特

47、殊正交多項式族：特殊正交多項式族： Legendre 多項式族：多項式族：1)( x 定義在定義在 1, 1上，上，kkkkkxdxdkxP)1(!21)(2 滿足：滿足： lklkPPklk1220),(xPP 10, 1由由 有遞推有遞推11)12()1( kkkkPxPkPk以以 Pn+1 的根為節點的求積公式稱為的根為節點的求積公式稱為Gauss-Legendre 公式公式。 Chebyshev 多項式族：多項式族：211)(xx 定義在定義在 1, 1上，上，) arccos( cos)(xkxTk Tn+1 的根為的根為 2212cosnkxkk = 0, , n以此為節點構造公式

48、以此為節點構造公式 1102)()(11nkkkxfAdxxfx稱為稱為 Gauss-Chebyshev 公式公式。注意到積分端點注意到積分端點 1 可能是積分可能是積分的的奇點奇點，用普通，用普通Newton-Cotes公公式在端點會出問題。而式在端點會出問題。而Gauss公公式可能避免此問題的發生。式可能避免此問題的發生。其它公式見教材其它公式見教材p.144-148注：注：一般一般a,b上的積分可化為上的積分可化為-1,1上特殊高斯公式進行計算。上特殊高斯公式進行計算。5 5 數值微分數值微分 數值微分的數值微分的概念概念 數值微分的數值微分的計算方法計算方法 原始概念近似原始概念近似:

49、 :中點法及外推法中點法及外推法 函數近似函數近似: :插值型的求導公式插值型的求導公式 函數相互關系轉化函數相互關系轉化: :利用數值積分求導利用數值積分求導 數值微分的數值微分的誤差分析誤差分析 泰勒展式估計泰勒展式估計 事后誤差估計事后誤差估計 基本關系轉化基本關系轉化 數值微分數值微分就是就是用函數值的線性組合近似函數在某點用函數值的線性組合近似函數在某點的導數值的導數值一、中點法和外推法一、中點法和外推法 按導數定義按導數定義 , 是差商是差商 當當 時的極限時的極限取取差商差商作為作為導數導數的近似值的近似值，建立簡單的數值微分方法，建立簡單的數值微分方法:)(0 xf hxfhx

50、f)()(00 0h( () )( () )( () )hxfhxfxf000 (6.1)向后差商近似導數向后差商近似導數(6.2)(6.3)中心差商近似導數中心差商近似導數( () )( () )( () )hhxfxfxf 000( () )( () )( () )hhxfhxfxf2000 ( ( ) )( () )( () )hhxfhxfhD200 容易看出，就精度而言，以（容易看出，就精度而言，以（6.3）式更為可取，稱）式更為可取，稱(6.4)為為 的的中點公式中點公式, 其中其中h為一增量，稱為為一增量，稱為步長步長 這種數值這種數值微分方法稱為微分方法稱為中點方法中點方法,

51、它是前兩種方法的算術平均它是前兩種方法的算術平均)(0 xf 分別將分別將在在 x=a 處做處做Taylor展開有展開有)(haf )(! 5)(! 4)(! 3)(! 2)()()() 5(5) 4(432afhafhafhafhafhafhaf代入代入D(h)得得 )(! 5)(! 3)()()5(42afhafhafhD,6)()(2MhhDaf 其中其中)(maxxfMhax 現在來考慮中點公式現在來考慮中點公式 的截斷誤差，的截斷誤差，hhafhafhD2)()()( (6.5)所以所以截斷誤差截斷誤差 )(! 5)(! 3)()() 5(42afhafhafhD(6.6)從截斷誤差

52、的角度來看，步長從截斷誤差的角度來看，步長h越小，計算結果越準確。且越小，計算結果越準確。且 但從計算角度看，但從計算角度看，h 越小，越小， f (a+h)與與 f (a- -h) 越接近，直越接近，直接相減會造成有效數字的嚴重損失接相減會造成有效數字的嚴重損失(參看第參看第1章第章第4節節)。因此。因此, 從舍入誤差的角度來看，步長從舍入誤差的角度來看，步長 h 不宜太小。不宜太小。 所以所以, 在在實際計算時實際計算時，通常，通常采用采用二分步長二分步長及及誤差事后估誤差事后估計法計法, 在變步長的過程中實現步長的自動選擇，在保證截斷在變步長的過程中實現步長的自動選擇，在保證截斷誤差滿足

53、的精度要求的前提下選取取盡可能大的步長。誤差滿足的精度要求的前提下選取取盡可能大的步長。kD(h)01239103.017652.791352.736442.722812.718282.71828 解解 這里采用的計算公式是這里采用的計算公式是 計算結果見表計算結果見表6.1，表中，表中 k 代表二分的次代表二分的次數，步長數，步長 。二分。二分 9 次得結果次得結果 D= 2.71828，它的每一數字都是有效數字，它的每一數字都是有效數字(所所求導數的準確值為求導數的準確值為e=2.7182818)。( ( ) )heehDhh211 kh28 . 0M 例例6.1 用用變步長變步長的中點方

54、法求的中點方法求 在在x=1處的導數值處的導數值,設取設取h= 0.8起算。起算。xe表表6.1 計算結果計算結果 210)()(hxfhD 4)()2()()(00 xfhDxfhD21041)()2(hxfhD )()2(34)()(0hDhDhDxf 事后誤差估計事后誤差估計 我們看到，中點公式具有如下形式我們看到，中點公式具有如下形式 (6.7)式中的系數均與步長無關。式中的系數均與步長無關。 6342210)()(hhhxfhD 若將步長二分，則有若將步長二分，則有 (6.8) 634221064116141)()2(hhhxfhD 41h 若令若令 (6.10)則進一步消去則進一步

55、消去D1(h)誤差主項誤差主項 , 有有 )(151)2(1516)(112hDhDhD 6102)()(hxfhD2141h則可消去誤差主項則可消去誤差主項 ，得，得 )(31)2(34)(1hDhDhD 624101)()(hhxfhD取取(6.7)與與(6.8)加加權平均權平均(6.9)重復同樣的手續重復同樣的手續，再導出下列加速公式，再導出下列加速公式 這種加速過程還可繼續下去。這種加速過程還可繼續下去。這種加速方法通常稱作這種加速方法通常稱作Richardson外推加速法外推加速法。 )(631)2(6364)(223hDhDhD 例例6.2 運用加速公式和加工例運用加速公式和加工例

56、6.1的結果。的結果。表表6.2 Richardson外推加速法計算結果外推加速法計算結果)(1hD)(2hD)(3hDh D(h)0.80.40.20.13.017652.791352.736442.722812.7159172.7181372.7182672.7182852.7182762.71828 解解 計算結果計算結果見表見表6.2。這里，。這里，加速的效果同樣加速的效果同樣是相當顯著的。是相當顯著的。hhaGafafe ee ee e 2)()()(21最小最小,步長步長h不宜太大不宜太大,也不宜太小也不宜太小. 其其最優步長最優步長應為應為 它表明它表明h 越小越小, 舍入誤差舍

57、入誤差)(af 越大越大, 故它是病態的故它是病態的. 用中點用中點公式計算公式計算)(af 的的誤差上界誤差上界為為,6)(2hMhhEe e 要使誤差要使誤差E(h)3opt/3Mhe e 則計算則計算 當當 f (a+h) 及及 f (a-h) 分別有舍入誤差分別有舍入誤差 e e1 及及 e e2 時，若令時，若令 e e)(af 的舍入誤差上界為的舍入誤差上界為 21,maxe ee e 注注: 中點心公式及其加速方法適合用表達式表示的函數。對中點心公式及其加速方法適合用表達式表示的函數。對于列表函數，則宜使用插值方法等導出數值求導公式。于列表函數，則宜使用插值方法等導出數值求導公式

58、。 二、插值型的求導公式二、插值型的求導公式 x x0 x1 x2 xn y y0 y1 y2 yn 對于列表函數對于列表函數 y = f (x)： 插值多項式插值多項式 y = Pn(x)作為它的近似作為它的近似, 我們取我們取統稱統稱插值型的求導公式插值型的求導公式)()(xPxfn )(xPn )(xf 的近似值，的近似值，作為作為建立的數值公式建立的數值公式 依據插值余項定理，依據插值余項定理，求導公式求導公式(6.11)的的余項余項為為)(dd)!1()()()!1()()()()1(11)1( nnnnnfxnxxnfxPxf式中式中. )()(01 niinxxx (6.11)

59、我們限定我們限定:求某個節點求某個節點 xk 上的導數值，上的導數值，上面的第二項變為零上面的第二項變為零，這時有余項公式，這時有余項公式)()!1()()()(1)1(knnknkxnfxPxf (6.12) 1 1兩點公式兩點公式 已給兩節點已給兩節點 x0, x1 上的函數值上的函數值 f (x0), f (x1)，做線性插值，做線性插值)()()(101001011xfxxxxxfxxxxxP 記記 x1 x0 = h，對上式兩端求導，有，對上式兩端求導，有)()(1)(101xfxfhxP )()(1)( )()(1)(01110101xfxfhxPxfxfhxP 于是有下列求導公式

60、：于是有下列求導公式：)(2)()(1)()(2)()(1)(011010 fhxfxfhxffhxfxfhxf 而利用余項公式而利用余項公式(6. 4)知，帶知，帶余項的兩點公式余項的兩點公式是：是： 下面我們僅僅考察下面我們僅僅考察節點處的導數值節點處的導數值為簡化討論為簡化討論, 假定假定所給的節點是等距的所給的節點是等距的 設已給出三節點設已給出三節點x0, xl=x0+h, x2=x0+2h上的函數值上的函數值,做二次插值做二次插值)()()()()()()()()()(2120210121012002010212xfxxxxxxxxxfxxxxxxxxxfxxxxxxxxxP 令令