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

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

3、地便獲得一種數(shù)值求積方法 如果用兩端點的如果用兩端點的“高度高度”f(a)與與f(b)的算術(shù)平均作為平均高度的算術(shù)平均作為平均高度f ( ) 的近似值，這樣導(dǎo)出的求積公式的近似值，這樣導(dǎo)出的求積公式 : 便是我們所熟悉的便是我們所熟悉的梯形公式梯形公式(trapezoidal rule). )()(2bfafabT 2)(bafabR2bac 而如果改用區(qū)間中點而如果改用區(qū)間中點 的的“高度高度”f (c)近似地取代平近似地取代平均高度均高度f ( )，則又可導(dǎo)出所謂，則又可導(dǎo)出所謂中矩形公式中矩形公式(今后簡稱矩形公式今后簡稱矩形公式)：(1.1)(1.2) 更一般地，我們可以在區(qū)間更一般地

4、，我們可以在區(qū)間a，b上適當(dāng)選取某些節(jié)點上適當(dāng)選取某些節(jié)點 xk ，然后然后用用 f (xk )加權(quán)平均得到平均高度加權(quán)平均得到平均高度 f ()的近似值的近似值，這樣構(gòu)造出的，這樣構(gòu)造出的求積公式具有下列形式求積公式具有下列形式式中式中 xk 稱為稱為求積節(jié)點求積節(jié)點；Ak 稱為稱為求積系數(shù)求積系數(shù)，亦稱為伴隨節(jié)點，亦稱為伴隨節(jié)點 xk 的的權(quán)權(quán)權(quán)權(quán)Ak 僅僅與節(jié)點僅僅與節(jié)點xk 的選取有關(guān)，而不依賴于被積函數(shù)的選取有關(guān)，而不依賴于被積函數(shù) f(x)的的具體形式具體形式 ban0kkkxfAdxxf)()(使積分公式具有通用性使積分公式具有通用性 這類數(shù)值積分方法通常稱作這類數(shù)值積分方法通常

5、稱作機械求積機械求積， 其特點是將積分求其特點是將積分求值問題歸結(jié)為函數(shù)值的計算，這就避開了牛頓值問題歸結(jié)為函數(shù)值的計算，這就避開了牛頓萊布尼茲公式萊布尼茲公式需要尋求原函數(shù)的困難需要尋求原函數(shù)的困難(1.3)二、二、代數(shù)精度的概念代數(shù)精度的概念 數(shù)值求積方法是近似方法，為要保證精度，我們自然希望求積數(shù)值求積方法是近似方法，為要保證精度，我們自然希望求積公式能對公式能對“盡可能多盡可能多”的函數(shù)準(zhǔn)確地成立，這就提出了所謂代數(shù)精的函數(shù)準(zhǔn)確地成立，這就提出了所謂代數(shù)精度的概念度的概念 定義定義 1 如果某個求積公式對于次數(shù)如果某個求積公式對于次數(shù)m的多項式均能準(zhǔn)確地成的多項式均能準(zhǔn)確地成立，但對于

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

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

8、故求積公式的形式為故求積公式的形式為解之得解之得 A0= h, A1=0, A2= h. 94 34 f(x)dx f(0) + f(2h)3h49h43h0由公式的構(gòu)造知由公式的構(gòu)造知,公式公式至少至少具有具有2次代數(shù)精度次代數(shù)精度; 而當(dāng)而當(dāng)f(x)=x3時時,公式的左邊公式的左邊= h4, 右邊右邊=18h4, 公式的左邊公式的左邊 右邊右邊,說明說明此公式對此公式對 f(x)=x3不能準(zhǔn)確成立不能準(zhǔn)確成立.因此因此,公式只具有公式只具有2次代數(shù)次代數(shù)精度精度.814三、三、求積公式的收斂性與穩(wěn)定性求積公式的收斂性與穩(wěn)定性 定理定理3表明，只要求積系數(shù)表明，只要求積系數(shù)Ak0 (k0,1

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

10、式若求積公式(13)中系數(shù)中系數(shù)Ak0 (k0，1，n)，則此求積公式是穩(wěn)定的則此求積公式是穩(wěn)定的近似近似計算計算 badxxfI)(思思路路利用利用插值多項式插值多項式 則積分易算。則積分易算。)()(xfxPn 在在a, b上取上取 a x0 x1 xn b，做，做 f 的的 n 次插值次插值多項式多項式 ，即得到，即得到 nkkknxlxfxL0)()()( babaknkkdxxlxfdxxf)()()(0Ak bakjxxxxkdxAjkj)()(由由 決定，決定，與與 無關(guān)。無關(guān)。節(jié)點節(jié)點 f (x)插值型積分公式插值型積分公式Interpolatory quadratureban

11、kkxnbanbanbankkkdxxxnfdxxRdxxLxfxfAdxxffR0)1(0)()!1()()()()()()(誤差誤差bandxxP)(2 插值型的求積公式插值型的求積公式與與Newton-Cotes 公式公式關(guān)鍵是關(guān)鍵是f(x)如果求積公式是插值型的如果求積公式是插值型的, 按余項式按余項式, 對于次數(shù)對于次數(shù) n的多項式的多項式 f (x)，其余項其余項R f 等于等于0，因而這時求積公式至少具有，因而這時求積公式至少具有n次代數(shù)精度次代數(shù)精度定理定理1：形如形如 的求積公式至少有的求積公式至少有 n 次代數(shù)精度次代數(shù)精度 該該公式為公式為插值型插值型（即：（即： ） n

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

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

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

15、 fhfR n = 4: Cotes Rule, 代數(shù)精度代數(shù)精度 = 5,)(9458)6(7 fhfR n 為為偶數(shù)階偶數(shù)階的的Newton-Cotes 公式至少有公式至少有 n+1 次代數(shù)精度。次代數(shù)精度。二、幾種低階求積公式的余項二、幾種低階求積公式的余項三、三、偶階求積公式的代數(shù)精度偶階求積公式的代數(shù)精度 作為插值型的求積公式，作為插值型的求積公式，n 階的牛頓階的牛頓-柯特斯公式至柯特斯公式至少具有少具有n 次的插值精度（定理次的插值精度（定理1）。實際的代數(shù)精度還可）。實際的代數(shù)精度還可進一步提高，一般地，可以證明下述定理進一步提高，一般地，可以證明下述定理: 定理定理 2 當(dāng)階

16、當(dāng)階 n 為偶數(shù)時，牛頓為偶數(shù)時，牛頓-柯特斯公式柯特斯公式至少有至少有 n+1 次代數(shù)精度次代數(shù)精度 . nkknknxfCabI0)()()(注：注：由公式知，當(dāng)由公式知，當(dāng)n8時，柯特斯系數(shù)出現(xiàn)負值，這時時，柯特斯系數(shù)出現(xiàn)負值，這時，初始數(shù)據(jù)誤差將會引起計算結(jié)果誤差增大，即計算不，初始數(shù)據(jù)誤差將會引起計算結(jié)果誤差增大，即計算不穩(wěn)定。因此，實際計算不用穩(wěn)定。因此，實際計算不用n8的牛頓的牛頓-柯特斯公式柯特斯公式 .估計估計截斷誤差截斷誤差為為解解 用用梯形公式梯形公式計算計算:=2.1835估計估計截斷誤差截斷誤差為為=0.6796用用Simpson公式公式計算：計算：=2. 0263例

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

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

19、程，可采用另一記法：令 n = 2n 為偶數(shù)，為偶數(shù)， 這時這時 ，有，有hkaxhnabhk ,2 )()(2)(4)(3 koddkevenkknbfxfxfafhS三、收斂速度與誤差估計：三、收斂速度與誤差估計：定義定義 若一個積分公式的誤差滿足若一個積分公式的誤差滿足 且且C 0，則則稱該公式是稱該公式是 p 階收斂階收斂的。的。 ChfRphlim0)(,)(,)(642hOChOShOTnnn例例4：計算計算dxx 10142 解：解： )1()(2)0(161718fxffTkk8kxk 其中其中= 3.138988494 )1()(2)(4)0(241oddeven4fxfxf

20、fSkk8kxk 其中其中= 3.141592502運算量基運算量基本相同本相同問題問題: 給定精度給定精度 e e，如何取，如何取 n ?例如：要求例如：要求 ，如何判斷，如何判斷 n = ？e e |nTI)()(122 fabhfR ? nkkhfh12)(12 )()(12)(1222afbfhdxxfhba 上述上述例例4中若要求中若要求 , 則則610| nTI622106| )0() 1 (|12| | hffhfRn00244949.0 h即：取即：取 n = 409通常采取將區(qū)間通常采取將區(qū)間不斷對分不斷對分的方法，即取的方法，即取 n = 2k上述上述例例4中中2k 409

21、 k = 9 時，時，T512 = 3.14159202例例4中：中：S4 = 3.141592502注意到區(qū)間再次對分時注意到區(qū)間再次對分時412)()(12122fRhafbffRnn 412 nnTITI)(3122nnnTTTI 可用來判斷迭代可用來判斷迭代是否停止。是否停止。(1)(2)(3)事后誤差估計事后誤差估計一、梯形法的遞推化一、梯形法的遞推化逐次分半法逐次分半法 上一節(jié)介紹的復(fù)化求積方法對提高精度是行之有效的，但上一節(jié)介紹的復(fù)化求積方法對提高精度是行之有效的，但在使用求積公式之前必須給出合適的步長，在使用求積公式之前必須給出合適的步長，步長步長取得取得太大精度太大精度難以保

22、證難以保證，步長太小步長太小則會導(dǎo)致則會導(dǎo)致計算量計算量的的增加增加，而事先給出一個，而事先給出一個恰當(dāng)?shù)牟介L又往往是困難的恰當(dāng)?shù)牟介L又往往是困難的 實際計算中常常實際計算中常常采用變步長的計算方案采用變步長的計算方案，即在步長，即在步長逐次分逐次分半半(即步長二分即步長二分)的過程中，反復(fù)利用復(fù)化求積公式進行計算，的過程中，反復(fù)利用復(fù)化求積公式進行計算，直至所求得的積分值滿足精度要求為止直至所求得的積分值滿足精度要求為止 設(shè)將求積區(qū)間設(shè)將求積區(qū)間a，b分成分成n等分，則一共有等分，則一共有n+1個分點，按個分點，按梯形公式計算積分值梯形公式計算積分值Tn，需要提供，需要提供n+1個函數(shù)值如果

23、將求積個函數(shù)值如果將求積區(qū)間再二分一次，則分點增至區(qū)間再二分一次，則分點增至2n+1個，我們來個，我們來考察考察二分二分前后兩前后兩個積分值個積分值之間的之間的聯(lián)系聯(lián)系4 4 龍貝格求積公式龍貝格求積公式逐次分半逐次分半計算計算方案方案的實現(xiàn)的實現(xiàn): 注意到每個子區(qū)間注意到每個子區(qū)間xk，xk+1經(jīng)過二分只增加了一個分經(jīng)過二分只增加了一個分點點 xk+1/2( xk+xk+1)/2，用復(fù)化梯形公式求得該子區(qū)間上的，用復(fù)化梯形公式求得該子區(qū)間上的積分值為積分值為 101021102110122)12(221)(221)(2)()(4nknnkknnkknkkknhkafhTxfhTxfhxfxf

24、hT)()(2)(4121 kkkxfxfxfh這里這里 代表二分前的步長代表二分前的步長. .將每個子區(qū)間上的積分值將每個子區(qū)間上的積分值相加得相加得nabh 二、龍貝格算法二、龍貝格算法).,()(212);,()(12)(222bafhabTIbafhabTIfRnnn 有有：根據(jù)復(fù)化梯形公式的余項表達式根據(jù)復(fù)化梯形公式的余項表達式. )(31.41)()(222nnnnnTTTITITIff 整理后可得：整理后可得：，則有，則有假定假定 可見，可見，利用利用兩種步長兩種步長計算的結(jié)果能估計截斷誤差計算的結(jié)果能估計截斷誤差.若將該截斷若將該截斷誤差加到計算結(jié)果中誤差加到計算結(jié)果中，nnn

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

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

27、積分值重復(fù)同樣的手續(xù)，用柯特斯法二分前后的兩個積分值Cn與與C2n的線性組合的結(jié)果可得到的線性組合的結(jié)果可得到龍貝格龍貝格(Romberg)求積公式求積公式nnnnnCCCCR631636414114423233 我們在變步長的過程中運用加速公式，就能將粗糙的梯我們在變步長的過程中運用加速公式，就能將粗糙的梯形值形值Tn逐步加工成精度較高的辛普森值逐步加工成精度較高的辛普森值Sn 、柯特斯值、柯特斯值Cn和龍和龍貝格值貝格值Rn .一般有：一般有：nnnSTT 1442nnnCSS 144222nnnRCC 144323 Romberg 算法：算法： e e ? e e ? e e ? T1

28、=)0(0T T8 =)3(0T T4 =)2(0T T2 =)1(0T S1 =)0(1T R1 =)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 區(qū)間等分?jǐn)?shù)區(qū)間等分?jǐn)?shù) 梯形序列梯形序列 辛普森序列辛普森序列 柯特斯序列柯特斯序列 龍貝格序列龍貝格序列 龍貝格求積算法可用下表來表示：龍貝格求積算法

29、可用下表來表示： 例例5 用龍貝格方法計算橢圓用龍貝格方法計算橢圓 x2/4 + y2 l 的周長，使結(jié)果的周長，使結(jié)果具有五位有效數(shù)字具有五位有效數(shù)字 分析分析 為便于計算，先將橢圓方程采用參數(shù)形式表示為便于計算，先將橢圓方程采用參數(shù)形式表示, ,再根再根據(jù)弧長公式將橢圓周長用積分形式表示由于計算結(jié)果要求具據(jù)弧長公式將橢圓周長用積分形式表示由于計算結(jié)果要求具有五位有效數(shù)字，因此需要估計所求積分值有幾位整數(shù)，從而有五位有效數(shù)字，因此需要估計所求積分值有幾位整數(shù)，從而確定所求積分值的絕對誤差限最后再應(yīng)用龍貝格方法計算積確定所求積分值的絕對誤差限最后再應(yīng)用龍貝格方法計算積分分 解解 令令 x 2c

30、osq q，y sinq q ， 則橢圓的周長為則橢圓的周長為Iyxl4d sin314d42022022 q qq qq qq qq q.10125. 01081)(1021)(4422d sin3124451202 fRIfRIlI的的截截斷斷誤誤差差為為故故計計算算，的的截截斷斷誤誤差差為為則則需需結(jié)結(jié)果果有有五五位位有有效效數(shù)數(shù)字字，有有一一位位整整數(shù)數(shù)，要要求求，因因此此由由于于 q qq q 下表給出了用龍貝格方法計算積分下表給出了用龍貝格方法計算積分I= 1+1+3sin2q q dx 的過程的過程. /20kk2kT212 kS22 kC32 kR4322 kkRR0 1 2.

31、356 1941 2 2.419 921 2.441 1632 4 2.422 103 2.422 830 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。三、理查森三、理查森(Richardso

32、n)外推加速法外推加速法 上面討論說明由梯形公式出發(fā)上面討論說明由梯形公式出發(fā), 將區(qū)間將區(qū)間a, b逐次二分逐次二分可提高求積公式的精度可提高求積公式的精度, 上述加速過程還可繼續(xù)下去上述加速過程還可繼續(xù)下去. 下面我們討論其下面我們討論其理論依據(jù)理論依據(jù). ,)(24221 llhhhIhT .)(122nabhbafhabTIn ， 22hTTn若記若記Tn = T(h), 當(dāng)區(qū)間當(dāng)區(qū)間a, b分為分為2n等分時等分時, 有有 , 則則可見可見I = T(h)的誤差為的誤差為O(h2). llhhhIhT2422121642 3)(24)(1hThThT 若記若記 ，則，則 將梯形公式按

33、余項展開將梯形公式按余項展開. 由誤差公式有由誤差公式有 62411)(hhIhT 6416262411hhIhT 顯然顯然T1(h)與與 I 近似的階為近似的階為O(h4) . 就是就是辛普森公式辛普森公式序列序列Sn, S2n, . .， 2),(11hThT這樣構(gòu)造的這樣構(gòu)造的 )(1412144)(11hThThTmmmmmm 則又可進一步從余項中則又可進一步從余項中消去消去 h4 項，這樣構(gòu)造出的項，這樣構(gòu)造出的 ，其實就是，其實就是柯特斯公式柯特斯公式序序列，它與列，它與 I 的逼近階為的逼近階為O(h6) . )(2hT)(151 21516)(112hThThT 若令若令 ,

34、一般地，若記一般地，若記T0(h) = T(h)，經(jīng)過，經(jīng)過m (m =1,2,)次加速次加速后，則有后，則有如此繼續(xù)下去，每加速一次，誤差的量級便提高如此繼續(xù)下去，每加速一次，誤差的量級便提高2階階. )21(141144)(1)1(1)()(0)()(0，次次加加速速值值，可可得得的的序序列列表表示示以以次次后后求求得得的的梯梯形形值值，且且表表示示二二分分設(shè)設(shè)以以 kTTTmTTkTkmmkmmmkmkkmk. ., 321.數(shù)數(shù)表表來來計計算算構(gòu)構(gòu)造造一一個個三三角角形形數(shù)數(shù)表表根根據(jù)據(jù)公公式式可可以以逐逐行行龍龍貝貝格格序序列列公公式式辛辛普普森森、柯柯特特斯斯、即即可可得得到到加加

35、速速、若若取取算算法法上上式式也也稱稱為為龍龍貝貝格格求求積積Tm Romberg 算法算法 可以證明，如果可以證明，如果 f (x) 充分光滑，那么充分光滑，那么T 數(shù)表每一列的數(shù)表每一列的元素及對角線元素均收斂到所求的積分值元素及對角線元素均收斂到所求的積分值 I ，即，即ITmITkmmkmk )()(lim)(lim，固固定定 理查德森理查德森外推法外推法利用利用低低階公式產(chǎn)生階公式產(chǎn)生高高精度的結(jié)果。精度的結(jié)果。設(shè)對于某一設(shè)對于某一 h 0，有公式，有公式 T0(h) 近似計算某一未知值近似計算某一未知值 I。由。由Taylor展開得到：展開得到： T0(h) I = 1 h + 2

36、 h2 + 3 h3 + i 與與 h 無關(guān)無關(guān)現(xiàn)將現(xiàn)將 h 對分，得：對分，得：( () )( () )( () ).)(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) 在構(gòu)造在構(gòu)造Newton-Cotes公式公式時，限定用積分區(qū)間的時，

37、限定用積分區(qū)間的等分點等分點作為求積節(jié)點作為求積節(jié)點，這樣做雖然使問題的處理過程得以簡化，但，這樣做雖然使問題的處理過程得以簡化，但同時也同時也限制了精度限制了精度。 求積公式含有求積公式含有2n+2個待定參數(shù)個待定參數(shù)xk、Ak(k0,1,n)若用若用待定系數(shù)法確定它們待定系數(shù)法確定它們, 則最好需要則最好需要2n+2個獨立的條件聯(lián)立方個獨立的條件聯(lián)立方程組求解程組求解, 從而易知求積公式的從而易知求積公式的最大代數(shù)精度最大代數(shù)精度可達到可達到2n+1次次. 在節(jié)點數(shù)目固定為在節(jié)點數(shù)目固定為n 的條件下，能否通過的條件下，能否通過適當(dāng)選取求積適當(dāng)選取求積節(jié)點節(jié)點xk的位置以及相應(yīng)的求積系數(shù)的

38、位置以及相應(yīng)的求積系數(shù)Ak，使求積公式，使求積公式具有盡可能高具有盡可能高(最高最高)的代數(shù)精度？的代數(shù)精度？ bankkkxfAxxf0)(d)(這類求積公式稱為這類求積公式稱為高斯高斯(Gauss)求積公式求積公式。4 高斯求積公式高斯求積公式 將節(jié)點將節(jié)點 x0 xn 以及系數(shù)以及系數(shù) A0 An 都作為待定系數(shù)。都作為待定系數(shù)。令令 f (x) = 1, x, x2, , x2n+1 代入可求解，得到的公式代入可求解，得到的公式具有具有2n+1 次代數(shù)精度。這樣的節(jié)點稱為次代數(shù)精度。這樣的節(jié)點稱為Gauss 點點，公式稱為公式稱為Gauss 型求積公式型求積公式。 baxxxfId)(

39、)( 為使問題更具一般性為使問題更具一般性,我們研究帶權(quán)積分我們研究帶權(quán)積分 (x)為權(quán)函數(shù)為權(quán)函數(shù), Ak(k0,1,n)為不依賴于為不依賴于f (x)的求積系數(shù)的求積系數(shù), xk (k0,1,n)為求積節(jié)點為求積節(jié)點. bamnkmkknmxxxxA)2.5(.12,1,0d)(0 要使要使(5.1)具有具有2n+1次代數(shù)精度，則需要滿足次代數(shù)精度，則需要滿足 bankkkxfAdxxfx0)()()( 構(gòu)造具有構(gòu)造具有2n+1次代數(shù)精度的求積公式次代數(shù)精度的求積公式(5.1) 從例中可看到求解非線性方程組從例中可看到求解非線性方程組(5.2)較復(fù)雜，通常較復(fù)雜，通常n2就很難求解故一般不

40、通過解方程就很難求解故一般不通過解方程(5.2)求求 xk 及及 Ak (k0,1, , n)例：例：求求 的的 2 點點 Gauss 公式。公式。dxxfx)(10 解：解：設(shè)設(shè) ，應(yīng)有，應(yīng)有 3 次代數(shù)精度。次代數(shù)精度。 101100)()()(xfAxfAdxxfx代入代入 f (x) = 1, x, x2, x3 31130092211200721100521032xAxAxAxAxAxAAA2776. 03891. 02899. 08212. 01010 AAxx不是線性方程組，不是線性方程組，不易求解。不易求解。 而從研究而從研究高斯點的基本特性高斯點的基本特性來著手解決來著手解決

41、Gauss 求積公式求積公式的構(gòu)造問題的構(gòu)造問題由插值型公式構(gòu)由插值型公式構(gòu)造知造知,關(guān)鍵求關(guān)鍵求xk，一、高斯點的基本特性一、高斯點的基本特性0)(d)()()(1 banxxxxxP 證明：證明： “” x0 xn 為為 Gauss 點點, 則公式則公式 至少有至少有 2n+1 次代數(shù)精度。次代數(shù)精度。 bankkkxfAdxxfx0)()()( 對任意次數(shù)對任意次數(shù)不大于不大于n 的多項式的多項式 Pm(x)， Pm(x) w(x)的次數(shù)的次數(shù)不大于不大于2n+1，則代入公式應(yīng)則代入公式應(yīng)精確成立精確成立： nkkkmkbamxwxPAdxxwxPx0)()()()()( 0= 0 “”

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

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

44、()()!22()(2)12(2)12(badxxwnfdxxwnfbanbaxn 二、高斯求積公式的余項二、高斯求積公式的余項三、高斯求積公式的穩(wěn)定性與收斂性三、高斯求積公式的穩(wěn)定性與收斂性 定理定理6 高斯求積公式高斯求積公式(5.1)的求積系數(shù)的求積系數(shù) Ak (k0,1,n)全是正的全是正的 由本定理及定理由本定理及定理2，則得，則得 推論推論 高斯求積公式高斯求積公式(5.1)是穩(wěn)定的是穩(wěn)定的. 定理定理7 設(shè)設(shè) f (x)C a，b，則高斯求積公式，則高斯求積公式(5.1)是收斂是收斂 的，即的，即 nkbakknxxxfxfA0.d)()()(lim 正交多項式族正交多項式族 0

45、, 1, , n, 有性質(zhì)：任意次數(shù)不大有性質(zhì)：任意次數(shù)不大于于n 的多項式的多項式 P(x) 必與必與 n+1 正交。正交。若取若取 w(x) 為其中的為其中的 n+1，則，則 n+1的根的根就是就是 Gauss 點。點。再解上例：再解上例： 101100)()()(xfAxfAdxxfxStep 1：構(gòu)造正交多項式構(gòu)造正交多項式 2設(shè)設(shè)cbxxxaxxx 2210)(,)(, 1)( 53 a0)(10 dxaxx0),(10 1021102100)(53(0),(0)(0),(dxcbxxxxdxcbxxx 215910 cb即：即：215910)(22 xxx 四、常用的高斯型求積公式

46、四、常用的高斯型求積公式Step 2：求求 2 = 0 的的 2 個根，即為個根，即為 Gauss 點點 x0 ，x1221/20)9/10(9/1021;0 xStep 3：代入代入 f (x) = 1, x 以求解以求解 A0 ，A1解解線性線性方程組，方程組，簡單。簡單。結(jié)果與前一方法相同：結(jié)果與前一方法相同：2776. 0,3891. 0,2899. 0,8212. 01010 AAxx 利用此公式計算利用此公式計算 的值的值 10dxexx2555. 1 10dxexx2899. 08212. 0102776. 03891. 010eeeAeAxx 注：注：構(gòu)造正交多項式也可以利用最

47、小二乘數(shù)據(jù)擬合中介紹構(gòu)造正交多項式也可以利用最小二乘數(shù)據(jù)擬合中介紹過的遞推式進行。過的遞推式進行。 特殊正交多項式族：特殊正交多項式族： Legendre 多項式族：多項式族：1)( x 定義在定義在 1, 1上，上，kkkkkxdxdkxP)1(!21)(2 滿足：滿足： lklkPPklk1220),(xPP 10, 1由由 有遞推有遞推11)12()1( kkkkPxPkPk以以 Pn+1 的根為節(jié)點的求積公式稱為的根為節(jié)點的求積公式稱為Gauss-Legendre 公式公式。 Chebyshev 多項式族：多項式族：211)(xx 定義在定義在 1, 1上，上，) arccos( co

48、s)(xkxTk Tn+1 的根為的根為 2212cosnkxkk = 0, , n以此為節(jié)點構(gòu)造公式以此為節(jié)點構(gòu)造公式 1102)()(11nkkkxfAdxxfx稱為稱為 Gauss-Chebyshev 公式公式。注意到積分端點注意到積分端點 1 可能是積分可能是積分的的奇點奇點，用普通，用普通Newton-Cotes公公式在端點會出問題。而式在端點會出問題。而Gauss公公式可能避免此問題的發(fā)生。式可能避免此問題的發(fā)生。其它公式見教材其它公式見教材p.144-148注：注：一般一般a,b上的積分可化為上的積分可化為-1,1上特殊高斯公式進行計算。上特殊高斯公式進行計算。5 5 數(shù)值微分?jǐn)?shù)

49、值微分 數(shù)值微分的數(shù)值微分的概念概念 數(shù)值微分的數(shù)值微分的計算方法計算方法 原始概念近似原始概念近似: :中點法及外推法中點法及外推法 函數(shù)近似函數(shù)近似: :插值型的求導(dǎo)公式插值型的求導(dǎo)公式 函數(shù)相互關(guān)系轉(zhuǎn)化函數(shù)相互關(guān)系轉(zhuǎn)化: :利用數(shù)值積分求導(dǎo)利用數(shù)值積分求導(dǎo) 數(shù)值微分的數(shù)值微分的誤差分析誤差分析 泰勒展式估計泰勒展式估計 事后誤差估計事后誤差估計 基本關(guān)系轉(zhuǎn)化基本關(guān)系轉(zhuǎn)化 數(shù)值微分?jǐn)?shù)值微分就是就是用函數(shù)值的線性組合近似函數(shù)在某點用函數(shù)值的線性組合近似函數(shù)在某點的導(dǎo)數(shù)值的導(dǎo)數(shù)值一、中點法和外推法一、中點法和外推法 按導(dǎo)數(shù)定義按導(dǎo)數(shù)定義 , 是差商是差商 當(dāng)當(dāng) 時的極限時的極限取取差商差商作為

50、作為導(dǎo)數(shù)導(dǎo)數(shù)的近似值的近似值，建立簡單的數(shù)值微分方法，建立簡單的數(shù)值微分方法:)(0 xf hxfhxf)()(00 0h( () )( () )( () )hxfhxfxf000 (6.1)向后差商近似導(dǎo)數(shù)向后差商近似導(dǎo)數(shù)(6.2)(6.3)中心差商近似導(dǎo)數(shù)中心差商近似導(dǎo)數(shù)( () )( () )( () )hhxfxfxf 000( () )( () )( () )hhxfhxfxf2000 ( ( ) )( () )( () )hhxfhxfhD200 容易看出，就精度而言，以（容易看出，就精度而言，以（6.3）式更為可取，稱）式更為可取，稱(6.4)為為 的的中點公式中點公式, 其中其

51、中h為一增量，稱為為一增量，稱為步長步長 這種數(shù)值這種數(shù)值微分方法稱為微分方法稱為中點方法中點方法, 它是前兩種方法的算術(shù)平均它是前兩種方法的算術(shù)平均)(0 xf 分別將分別將在在 x=a 處做處做Taylor展開有展開有)(haf )(! 5)(! 4)(! 3)(! 2)()()() 5(5) 4(432afhafhafhafhafhafhaf代入代入D(h)得得 )(! 5)(! 3)()()5(42afhafhafhD,6)()(2MhhDaf 其中其中)(maxxfMhax 現(xiàn)在來考慮中點公式現(xiàn)在來考慮中點公式 的截斷誤差，的截斷誤差，hhafhafhD2)()()( (6.5)所以

52、所以截斷誤差截斷誤差 )(! 5)(! 3)()() 5(42afhafhafhD(6.6)從截斷誤差的角度來看，步長從截斷誤差的角度來看，步長h越小，計算結(jié)果越準(zhǔn)確。且越小，計算結(jié)果越準(zhǔn)確。且 但從計算角度看，但從計算角度看，h 越小，越小， f (a+h)與與 f (a- -h) 越接近，直越接近，直接相減會造成有效數(shù)字的嚴(yán)重損失接相減會造成有效數(shù)字的嚴(yán)重損失(參看第參看第1章第章第4節(jié)節(jié))。因此。因此, 從舍入誤差的角度來看，步長從舍入誤差的角度來看，步長 h 不宜太小。不宜太小。 所以所以, 在在實際計算時實際計算時，通常，通常采用采用二分步長二分步長及及誤差事后估誤差事后估計法計法,

53、 在變步長的過程中實現(xiàn)步長的自動選擇，在保證截斷在變步長的過程中實現(xiàn)步長的自動選擇，在保證截斷誤差滿足的精度要求的前提下選取取盡可能大的步長。誤差滿足的精度要求的前提下選取取盡可能大的步長。kD(h)01239103.017652.791352.736442.722812.718282.71828 解解 這里采用的計算公式是這里采用的計算公式是 計算結(jié)果見表計算結(jié)果見表6.1，表中，表中 k 代表二分的次代表二分的次數(shù)，步長數(shù)，步長 。二分。二分 9 次得結(jié)果次得結(jié)果 D= 2.71828，它的每一數(shù)字都是有效數(shù)字，它的每一數(shù)字都是有效數(shù)字(所所求導(dǎo)數(shù)的準(zhǔn)確值為求導(dǎo)數(shù)的準(zhǔn)確值為e=2.7182

54、818)。( ( ) )heehDhh211 kh28 . 0M 例例6.1 用用變步長變步長的中點方法求的中點方法求 在在x=1處的導(dǎo)數(shù)值處的導(dǎo)數(shù)值,設(shè)取設(shè)取h= 0.8起算。起算。xe表表6.1 計算結(jié)果計算結(jié)果 210)()(hxfhD 4)()2()()(00 xfhDxfhD21041)()2(hxfhD )()2(34)()(0hDhDhDxf 事后誤差估計事后誤差估計 我們看到，中點公式具有如下形式我們看到，中點公式具有如下形式 (6.7)式中的系數(shù)均與步長無關(guān)。式中的系數(shù)均與步長無關(guān)。 6342210)()(hhhxfhD 若將步長二分，則有若將步長二分，則有 (6.8) 63

55、4221064116141)()2(hhhxfhD 41h 若令若令 (6.10)則進一步消去則進一步消去D1(h)誤差主項誤差主項 , 有有 )(151)2(1516)(112hDhDhD 6102)()(hxfhD2141h則可消去誤差主項則可消去誤差主項 ，得，得 )(31)2(34)(1hDhDhD 624101)()(hhxfhD取取(6.7)與與(6.8)加加權(quán)平均權(quán)平均(6.9)重復(fù)同樣的手續(xù)重復(fù)同樣的手續(xù)，再導(dǎo)出下列加速公式，再導(dǎo)出下列加速公式 這種加速過程還可繼續(xù)下去。這種加速過程還可繼續(xù)下去。這種加速方法通常稱作這種加速方法通常稱作Richardson外推加速法外推加速法。

56、 )(631)2(6364)(223hDhDhD 例例6.2 運用加速公式和加工例運用加速公式和加工例6.1的結(jié)果。的結(jié)果。表表6.2 Richardson外推加速法計算結(jié)果外推加速法計算結(jié)果)(1hD)(2hD)(3hDh D(h)0.80.40.20.13.017652.791352.736442.722812.7159172.7181372.7182672.7182852.7182762.71828 解解 計算結(jié)果計算結(jié)果見表見表6.2。這里，。這里，加速的效果同樣加速的效果同樣是相當(dāng)顯著的。是相當(dāng)顯著的。hhaGafafe ee ee e 2)()()(21最小最小,步長步長h不宜太大

57、不宜太大,也不宜太小也不宜太小. 其其最優(yōu)步長最優(yōu)步長應(yīng)為應(yīng)為 它表明它表明h 越小越小, 舍入誤差舍入誤差)(af 越大越大, 故它是病態(tài)的故它是病態(tài)的. 用中點用中點公式計算公式計算)(af 的的誤差上界誤差上界為為,6)(2hMhhEe e 要使誤差要使誤差E(h)3opt/3Mhe e 則計算則計算 當(dāng)當(dāng) f (a+h) 及及 f (a-h) 分別有舍入誤差分別有舍入誤差 e e1 及及 e e2 時，若令時，若令 e e)(af 的舍入誤差上界為的舍入誤差上界為 21,maxe ee e 注注: 中點心公式及其加速方法適合用表達式表示的函數(shù)。對中點心公式及其加速方法適合用表達式表示的

58、函數(shù)。對于列表函數(shù)，則宜使用插值方法等導(dǎo)出數(shù)值求導(dǎo)公式。于列表函數(shù)，則宜使用插值方法等導(dǎo)出數(shù)值求導(dǎo)公式。 二、插值型的求導(dǎo)公式二、插值型的求導(dǎo)公式 x x0 x1 x2 xn y y0 y1 y2 yn 對于列表函數(shù)對于列表函數(shù) y = f (x)： 插值多項式插值多項式 y = Pn(x)作為它的近似作為它的近似, 我們?nèi)∥覀內(nèi)〗y(tǒng)稱統(tǒng)稱插值型的求導(dǎo)公式插值型的求導(dǎo)公式)()(xPxfn )(xPn )(xf 的近似值，的近似值，作為作為建立的數(shù)值公式建立的數(shù)值公式 依據(jù)插值余項定理，依據(jù)插值余項定理，求導(dǎo)公式求導(dǎo)公式(6.11)的的余項余項為為)(dd)!1()()()!1()()()()1

59、(11)1( nnnnnfxnxxnfxPxf式中式中. )()(01 niinxxx (6.11) 我們限定我們限定:求某個節(jié)點求某個節(jié)點 xk 上的導(dǎo)數(shù)值，上的導(dǎo)數(shù)值，上面的第二項變?yōu)榱闵厦娴牡诙椬優(yōu)榱?，這時有余項公式，這時有余項公式)()!1()()()(1)1(knnknkxnfxPxf (6.12) 1 1兩點公式兩點公式 已給兩節(jié)點已給兩節(jié)點 x0, x1 上的函數(shù)值上的函數(shù)值 f (x0), f (x1)，做線性插值，做線性插值)()()(101001011xfxxxxxfxxxxxP 記記 x1 x0 = h，對上式兩端求導(dǎo)，有，對上式兩端求導(dǎo)，有)()(1)(101xfxf

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