數(shù)值第4章數(shù)值積分與數(shù)值微分

1、數(shù)理學院數(shù)理學院SCHOOL OF MATHEMATICS AND PHYSICS4 .1 引言引言4 .2 牛頓牛頓-柯特斯公式柯特斯公式4 .3 復化求積公式復化求積公式4 .4 龍貝格求積公式龍貝格求積公式4 .5 高斯求積公式高斯求積公式4 .6 數(shù)值微分數(shù)值微分第四章第四章 數(shù)值積分與數(shù)值微分數(shù)值積分與數(shù)值微分一、數(shù)值求積的基本思想一、數(shù)值求積的基本思想)()()(aFbFdxxfba 積分積分 只要找到被積函數(shù)只要找到被積函數(shù) f (x)原函數(shù)原函數(shù)F(x)，便有，便有牛頓牛頓萊布尼茲萊布尼茲(NewtonLeibniz)公式公式 baxxfId)(實際困難實際困難：大量的被積函數(shù)

2、（：大量的被積函數(shù)（ , sin x2 等）等）, 找不到用初等函找不到用初等函數(shù)表示的原函數(shù)數(shù)表示的原函數(shù)；另外；另外, 當當f (x)是（測量或數(shù)值計算出的）一張是（測量或數(shù)值計算出的）一張數(shù)據(jù)表時，數(shù)據(jù)表時，牛頓牛頓萊布尼茲公式萊布尼茲公式也也不能直接運用不能直接運用。xxsin 積分中值定理：在積分中值定理：在a, b內(nèi)存在一點內(nèi)存在一點 ，有，有 f( )成立。成立。 )(d)(abxxfba 4.1 引言引言 就是說就是說, 底為底為b- -a 而高為而高為f( )的的矩形面積矩形面積恰恰等于所求等于所求曲邊梯形的面積曲邊梯形的面積 .問題問題 在于點在于點的具體位置一般是不知道的

3、，因而的具體位置一般是不知道的，因而 難以準確算出難以準確算出 f( )的值的值我們將我們將f ( )稱為區(qū)間稱為區(qū)間a, b上的平均高度這樣上的平均高度這樣,只要對只要對平均高度平均高度f( )提供一種算法提供一種算法,相應地便獲得一種數(shù)值求積方法相應地便獲得一種數(shù)值求積方法 如果用兩端點的如果用兩端點的“高度高度”f(a)與與f(b)的算術平均作為平均高度的算術平均作為平均高度f ( ) 的近似值，這樣導出的求積公式的近似值，這樣導出的求積公式 : 便是我們所熟悉的便是我們所熟悉的梯形公式梯形公式 . )()(2bfafabT 2)(bafabR2bac 而如果改用區(qū)間中點而如果改用區(qū)間中

4、點 的的“高度高度”f (c)近似地取代平近似地取代平均高度均高度f ( )，則又可導出所謂，則又可導出所謂中矩形公式中矩形公式(今后簡稱矩形公式今后簡稱矩形公式)：(1.1)(1.2) 更一般地，我們可以在區(qū)間更一般地，我們可以在區(qū)間a，b上適當選取某些節(jié)點上適當選取某些節(jié)點 xk ，然后然后用用 f (xk )加權平均得到平均高度加權平均得到平均高度 f ()的近似值的近似值，這樣構造出的，這樣構造出的求積公式具有下列形式求積公式具有下列形式式中式中 xk 稱為稱為求積節(jié)點求積節(jié)點；Ak 稱為稱為求積系數(shù)求積系數(shù)，亦稱為伴隨節(jié)點，亦稱為伴隨節(jié)點 xk 的的權權權權Ak 僅僅與節(jié)點僅僅與節(jié)點

5、xk 的選取有關，而不依賴于被積函數(shù)的選取有關，而不依賴于被積函數(shù) f(x)的的具體形式具體形式 ban0kkkxfAdxxf)()(使積分公式具有通用性使積分公式具有通用性 這類數(shù)值積分方法通常稱作能這類數(shù)值積分方法通常稱作能機械求積機械求積， 其特點是將積分其特點是將積分求值問題歸結為函數(shù)值的計算，這就避開了牛頓求值問題歸結為函數(shù)值的計算，這就避開了牛頓萊布尼茲公萊布尼茲公式需要尋求原函數(shù)的困難式需要尋求原函數(shù)的困難(1.3)二、二、代數(shù)精度的概念代數(shù)精度的概念 數(shù)值求積方法是近似方法，為要保證精度，我們自然希望求積數(shù)值求積方法是近似方法，為要保證精度，我們自然希望求積公式能對公式能對“盡

6、可能多盡可能多”的函數(shù)準確地成立，這就提出了所謂代數(shù)精的函數(shù)準確地成立，這就提出了所謂代數(shù)精度的概念度的概念 定義定義 1 如果某個求積公式對于次數(shù)如果某個求積公式對于次數(shù)m的多項式均能準確地成的多項式均能準確地成立，但對于立，但對于m+1次多項式就不準確成立，則稱該求積公式具有次多項式就不準確成立，則稱該求積公式具有m次代次代數(shù)精度數(shù)精度(或代數(shù)精確度或代數(shù)精確度) 一般地，欲使求積公式一般地，欲使求積公式 具有具有m次代數(shù)次代數(shù)精度，只要令它對于精度，只要令它對于f (x) = 1，x，xm 都能準確成立，這就要求都能準確成立，這就要求 bankkkxfAxxf0)(d)( . )(11;

7、)(21;1122mmmkkkkkabmxAabxAabA例例1： 考察其代數(shù)精度。考察其代數(shù)精度。 f(x)abf(a)f(b)梯形公式梯形公式解：解：逐次檢查公式是否精確成立逐次檢查公式是否精確成立代入代入 f(x) = 1： baabdx111 2 ab=代入代入 f(x) = x ：=代入代入 f(x) = x2 ： 222abbadxx 2baab 3233abbadxx 222baab 代數(shù)精度代數(shù)精度 = 1)()(2)(bfafabdxxfba 例例2 試構造形如試構造形如 f(x)dx A0f(0)+ A1f(h)+ A2f(2h) 的數(shù)值的數(shù)值求積公式求積公式,使其代數(shù)精度

8、盡可能高使其代數(shù)精度盡可能高,并指出其代數(shù)精度的階數(shù)并指出其代數(shù)精度的階數(shù).3h0解解: 令公式對令公式對 f(x)=1,x, x2 均準確成立均準確成立,則有則有3h=A0+ A1+ A2h2=0 + A1h+ A22h9h3=0 + A1h2+ A24h229故求積公式的形式為故求積公式的形式為解之得解之得 A0= h, A1=0, A2= h. 34 34 f(x)dx f(0) + f(2h)3h43h43h0由公式的構造知由公式的構造知,公式公式至少至少具有具有2次代數(shù)精度次代數(shù)精度; 而當而當f(x)=x3時時,公式的左邊公式的左邊= h4, 右邊右邊=18h4, 公式的左邊公式的

9、左邊 右邊右邊,說明說明此公式對此公式對 f(x)=x3不能準確成立不能準確成立.因此因此,公式只具有公式只具有2次代數(shù)次代數(shù)精度精度.814三、三、求積公式的收斂性與穩(wěn)定性求積公式的收斂性與穩(wěn)定性 定理定理3表明，只要求積系數(shù)表明，只要求積系數(shù)Ak0 (k0,1,n)，就能保證，就能保證計算的穩(wěn)定性計算的穩(wěn)定性 定義定義2 在求積公式在求積公式 中，若中，若 其中其中 ，則稱求積公式是收斂的，則稱求積公式是收斂的 由于計算由于計算 f (xk)可能有誤差可能有誤差,實際得到實際得到 定義定義3 對任給對任給 e e 0，若，若 (k=0,1, ,n), 就有就有 , 則稱求積公式是穩(wěn)定的則稱

10、求積公式是穩(wěn)定的. bankkkxfAxxf0)(d)(e e |)(|00nkkknkkkfAxfA)(max11 iinixxhkkfxf)(0 ，只只要要 .)(,kkkkfxff 即即 bankkkhnxxfxfAd)()(lim00 定理定理3 若求積公式若求積公式(13)中系數(shù)中系數(shù)Ak0 (k0，1，n)，則此求積公式是穩(wěn)定的則此求積公式是穩(wěn)定的近似近似計算計算 badxxfI)(思思路路利用利用插值多項式插值多項式 則積分易算。則積分易算。)()(xfxPn 在在a, b上取上取 a x0 x1 x=1,2； y=exp(1./x) I=1/2*(y(1)+y(2)*1I =

11、2.1835 R1=(2-1)3/12*8.1548R1 = 0.6796估計估計截斷誤差截斷誤差為為用用Simpson公式公式計算：計算：=2. 0263=198.4306890. 0)(max2880) 12()4(2152xfRx=0.06890 x=1:0.5:2; y=exp(1./x)y = 2.7183 1.9477 1.6487 I2=(2-1)/6*(y(1)+4*y(2)+y(3)I2 = 2.0263 f41=exp(1/x1)*(1./x18+12./x17+36./x16+24./x15)f41 = 198.4346 R2=(2-1)5/2880*f41R2 = 0.

12、06894.3 4.3 復合求積公式復合求積公式高次插值有高次插值有Runge 現(xiàn)象現(xiàn)象，故采用分段低次插值，故采用分段低次插值 分段低次合成的分段低次合成的 Newton-Cotes 復合復合求積公式。求積公式。一、復合梯形公式一、復合梯形公式:),., 0(,nkhkaxnabhk 在每個在每個 上用梯形公式：上用梯形公式：,1 kkxx 11)()(2)(2nkkbfxfafh bankkkxfxfhdxxf101)()(2)(=Tn),(),()(12)()(12)(122102103bafabhnfabhfhfRnkknkk /*中值定理中值定理*/1,., 0,)()(2)(111

13、 nkxfxfxxdxxfkkxxkkkk二、復合辛普森公式二、復合辛普森公式:),., 0(,nkhkaxnabhk )()(4)(6)(1211 kkkxxxfxfxfhdxxfkkkx21 kx1 kx44444 )()(2)(4)(6)(1010121 nknkkkbabfxfxfafhdxxf= Sn)(2180)4(4 fhabfR 注：注：為方便編程，可采用另一記法：令為方便編程，可采用另一記法：令 n = 2n 為偶數(shù)，為偶數(shù)， 這時這時 ，有，有hkaxhnabhk ,2 )()(2)(4)(3 koddkevenkknbfxfxfafhS三、收斂速度與誤差估計：三、收斂速度

14、與誤差估計：定義定義 若一個積分公式的誤差滿足若一個積分公式的誤差滿足 且且C 0，則則稱該公式是稱該公式是 p 階收斂階收斂的。的。 ChfRphlim0)(,)(,)(642hOChOShOTnnn例例4：計算計算dxx 10142 解：解： )1()(2)0(161718fxffTkk8kxk 其中其中= 3.138988494 )1()(2)(4)0(241oddeven4fxfxffSkk8kxk 其中其中= 3.141592502運算量基運算量基本相同本相同function t=rctrap(fun,a,b,n)%復化梯形公式%n等分%a,b區(qū)間的左右端點h=(b-a)/n;t=0

15、;for i=1:n-1 t=t+h*feval(fun,i*h+a);endt=t+0.5*h*(feval(fun,a+eps)+feval(fun,b)function f=fun4(x)f=4/(1+x2) t8=rctrap(fun4,0,1,8)function s=sptrap(fun,a,b,n)% n，對應的等分點為2nh=(b-a)/(2*n);s=0;for i=1:ns=s+feval(fun,(2*i-1)*h+a)*4;endfor i=1:n-1s=s+feval(fun,(2*i)*h+a)*2;ends=s+feval(fun,a+eps)+feval(fun

16、,b);s=s*h/3 s4=sptrap(fun4,0,1,2)問題問題: 給定精度給定精度 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 k = 9

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

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

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

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

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

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

23、重復同樣的手續(xù)，用柯特斯法二分前后的兩個積分值Cn與與C2n的線性組合的結果可得到的線性組合的結果可得到龍貝格龍貝格(Romberg)求積公式求積公式nnnnnCCCCR631636414114423233 我們在變步長的過程中運用加速公式，就能將粗糙的梯我們在變步長的過程中運用加速公式，就能將粗糙的梯形值形值Tn逐步加工成精度較高的辛普森值逐步加工成精度較高的辛普森值Sn 、柯特斯值、柯特斯值Cn和龍和龍貝格值貝格值Rn . Romberg 算法：算法： e e ? e e ? e e ? T1 =)0(0T T8 =)3(0T T4 =)2(0T T2 =)1(0T S1 =)0(1T R

24、1 =)0(3T S2 =)1(1T C1 =)0(2T C2 =)1(2T S4 =)2(1T一般有：一般有：nnnSTT 1442nnnCSS 144222nnnRCC 144323Romberg 序列序列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ū)間等分數(shù)區(qū)間等分數(shù) 梯形序列梯形序列 辛普森序列辛普森序列 柯特斯序列柯特斯序列 龍貝格序列龍貝格序列 龍貝格求積算法可用下表來表示：龍貝格求積算法可用下

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

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

27、 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。三、理查森三、理查森(Richardson)外

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

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

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

31、若若取取算算法法上上式式也也稱稱為為龍龍貝貝格格求求積積Tm Romberg 算法算法 可以證明，如果可以證明，如果 f (x) 充分光滑，那么充分光滑，那么T 數(shù)表每一列的數(shù)表每一列的元素及對角線元素均收斂到所求的積分值元素及對角線元素均收斂到所求的積分值 I ，即，即ITmITkmmkmk )()(lim)(lim，固固定定例用龍貝格算法計算積分例用龍貝格算法計算積分.102/3dxxI function s,n=rbg1(fun,a,b,eps) if nargineps) h=(b-a)/2(k-1); w=0; if(h=0) for i=1:(2(k-1)-1) w=w+f (a+

32、i*h); end t(k,1)=h*(f(a)/2+w+f(b)/2) for l=2:kfor i=1:(k-l+1) t(i,l)=(4(l-1)*t(i+1,l-1)-t(i,l-1)/(4(l-1)-1)endends=t(1,k); s0=(t(1,k-1); k=k+1; n=k; else s=s0; n=-k; end end function f=fun(x)f=x(3/2) R=rbg1(fun,0,1,1e-6)自適應積分方法自適應積分方法(補充內(nèi)容補充內(nèi)容) 復合求積方法通常適用于被積函數(shù)變化不太大的積分復合求積方法通常適用于被積函數(shù)變化不太大的積分,如果在求積區(qū)間中

33、被積函數(shù)變化很大如果在求積區(qū)間中被積函數(shù)變化很大,有的部分函數(shù)值變化有的部分函數(shù)值變化劇烈劇烈,另一部分變化平緩另一部分變化平緩.這時統(tǒng)一將區(qū)間等分用復合求積公這時統(tǒng)一將區(qū)間等分用復合求積公式計算積分工作量大式計算積分工作量大. 要達到誤差要求對變化劇烈部分必須將區(qū)間細分要達到誤差要求對變化劇烈部分必須將區(qū)間細分,而平而平緩部分則可用大步長緩部分則可用大步長.針對被積函數(shù)在區(qū)間不同情形采用不針對被積函數(shù)在區(qū)間不同情形采用不用的步長用的步長,使得在滿足精度前提下積分計算工作量盡可能小使得在滿足精度前提下積分計算工作量盡可能小. 在不同區(qū)間上預測被積函數(shù)變化的劇烈程度確定相應在不同區(qū)間上預測被積函

34、數(shù)變化的劇烈程度確定相應步長步長,這種方法稱為自適應積分方法這種方法稱為自適應積分方法設給定精度要求設給定精度要求, 0 e e計算積分計算積分 badxxffI)()(的近似值的近似值.先取步長先取步長h=b-a,應用辛普森公式有應用辛普森公式有),(),()2(180),()()()4(4bafhabbaSdxxffIba (5.1)其中其中).()2(4)(6),(bfbafafhbaS 若把區(qū)間若把區(qū)間a,b對分對分,步長步長,222abhh 在每個小區(qū)間上用辛普森在每個小區(qū)間上用辛普森公式公式,則得則得),(),()2(180),()()4(422bafhabbaSfI (5.2)其

35、中其中).()43(4)2(6),2(),2()4(4)(6)2,(),2()2,(),(222bfhafhafhbbaShafhafafhbaaSbbaSbaaSbaS 實際上實際上(5.2)式即為式即為),(),()4(180),()()4(42bafhabbaSfI 與與(5.1)式比較，若式比較，若)()4(xf在在(a,b)上變化不大上變化不大,可假定可假定)()()4()4( ff 從而可得從而可得).()2(180),(),(1516)4(42 fhabbaSbaS 與與(5.2)比較比較,則得則得|,|151| ),(),(|151| ),()(|2122SSbaSbaSbaS

36、fI 如如果果有有).,(),(221baSSbaSS 這里這里,15|21e e SS則可得到則可得到,| ),()(|2e e baSfI此時此時,可取可取S2(a,b)作為作為I(f)的近似的近似,則可達到給定的誤差精度則可達到給定的誤差精度.(5.3)若不等式若不等式(5.3)不成立不成立,則應分別對子區(qū)間則應分別對子區(qū)間,22,bbabaa 及及再用辛普森公式再用辛普森公式,此時步長此時步長h3=0.5h2,得到得到)2()2,(33bbaSbaaS，及及 只要分別考察只要分別考察2| )2,()(|3e e baaSfI及及2| ),2()(|3e e bbaSfI是否成立是否成立

37、.對滿足要求的區(qū)間不再細分對滿足要求的區(qū)間不再細分,對不滿足要求的還要繼續(xù)上述過程對不滿足要求的還要繼續(xù)上述過程.最后還要應用龍貝格法則求出相應區(qū)間的積分近似值最后還要應用龍貝格法則求出相應區(qū)間的積分近似值.例例7 計算積分計算積分,112 . 02dxx 若用復合辛普森法若用復合辛普森法(3.5)式計算結果見表式計算結果見表此處此處hn即為公式中的即為公式中的h,積分精確值為積分精確值為4).nhnSn|Sn-Sn-1|10.84.9481480.7611120.44.1870370.16281930.24.0242180.02205440.14.0021640.00201050.054.0

38、00154計算到計算到|Sn-Sn-1|0.02為止為止,此時此時I(f)的近似值的近似值S50.2,1=4.000154,若再用龍貝格法則得到若再用龍貝格法則得到00002. 4151 , 2 . 0455 SSSRS整個計算將區(qū)間整個計算將區(qū)間32等分等分,計算計算33個個f(x)的值的值.現(xiàn)在若用自適應積分法現(xiàn)在若用自適應積分法,當當h2=0.4時有時有S20.2,0.6=3.51851852,S2 0.6,1=0.66851852,由于由于S2= S20.2,1= S20.2,0.6+ S20.6,1=4.187037,|S1-S2|=0.761111大于允許誤差大于允許誤差0.02,

39、故對故對0.2,0.6及及 0.6,1兩區(qū)間再用兩區(qū)間再用h3=h2/2做積分做積分.先計算先計算0.6,1的積分的積分S3 0.6,0.8=0.41678477, S3 0.8,1=0.25002572.由于由于S20.6,1-( S30.6,0.8+ S30.8,1)=0.66851852-0.66681049=0.001708 小于允許誤差小于允許誤差0.01,故在故在0.6,1區(qū)間的積分值為區(qū)間的積分值為66669662. 0)66851852. 066681049. 0 0 1 , 6 . 0 RS再計算再計算0.2,0.6的積分的積分, S2 0.2,0.

40、6=3.51851852,而對而對h3=h2/2得得S3 0.2,0.4=2.52314815, S3 0.4,0.6=0.83425926.由于由于S20.2,0.6-( S30.2,0.4+ S30.4,0.6)=0. 161111大于允許大于允許誤差誤差0.01,因此還要分別計算因此還要分別計算0.2,0.4及及0.4,0.6 積分積分.當當h4=h3/2時可求得時可求得S4 0.4,0.5=0.50005144, S4 0.5,0.6=0.33334864.而而S30.4,0.6-( S40.4,0.5+ S40.5,0.6)=0.000859 小于允許小于允許誤差誤差0.01,故在故

41、在0.4,0.6區(qū)間的積分值為區(qū)間的積分值為.8333428. 06 . 0 , 4 . 0 RS而對而對0.2,0.4的積分的積分, S3 0.2,0.4- S4 0.2,0.4不小于不小于0.005,故故還要分別計算還要分別計算 0.2,0.3及及0.3,0.4的積分的積分,其中其中 S40.3,0.4=0.83356954,當當h5=h4/2可求得可求得S5 0.3,0.35=0.47620166, S5 0.35,0.4=0.35714758,且且S40.3,0.4-(S5 0.3,0.35+ S5 0.35,0.4)=0.000220小于允小于允許誤差許誤差0.0025,故有故有.8

42、3333492. 04 . 0 , 3 . 0 RS最后子區(qū)間最后子區(qū)間0.2,0.3的積分可檢驗出它的誤差小于的積分可檢驗出它的誤差小于0.0025,且可得且可得.666686. 13 . 0 , 2 . 0 RS將以上各區(qū)間的積分近似值相加可得到將以上各區(qū)間的積分近似值相加可得到.00005957. 41 , 6 . 06 . 0 , 4 . 04 . 0 , 3 . 03 . 0 , 2 . 0)( RSRSRSRSfI它一共只需要計算它一共只需要計算17個個f(x)的值的值. 在構造在構造Newton-Cotes公式公式時，限定用積分區(qū)間的時，限定用積分區(qū)間的等分點等分點作為求積節(jié)點作

43、為求積節(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 的條件下，能否通過的條件下，能否通過適當選取求積適當選取求積節(jié)點節(jié)點xk的位置以及相應的求積系數(shù)的位置以及相應的求積系數(shù)Ak，使求積公式，使

44、求積公式具有盡可能高具有盡可能高(最高最高)的代數(shù)精度？的代數(shù)精度？ bankkkxfAxxf0)(d)(這類求積公式稱為這類求積公式稱為高斯高斯(Gauss)求積公式求積公式。4.5 高斯求積公式高斯求積公式 將節(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)()( 為使問題更具一般性為使問題更具一

45、般性,我們研究帶權積分我們研究帶權積分 (x)為權函數(shù)為權函數(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)()()( 構造具有構造具有2n+1次代數(shù)精度的求積公式次代數(shù)精度的求積公式(5.1) 從例中可看到求解非線性方程組從例中可看到求解非線性方程組(5.2)較復雜，通常較復雜，通常n2就很難求解故一般不通過解方程就很難求解故一般不通過解方程

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

47、的構造問題由插值型公式構由插值型公式構造知造知,關鍵求關鍵求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，則代入公式應則代入公式應精確成立精確成立： nkkkmkbamxwxPAdxxwxPx0)()()()()( 0= 0 “” 要證明要證明 x0 xn 為為 Gauss 點，即要證公式對任意次點，即要證

48、公式對任意次數(shù)數(shù)不大于不大于2n+1 的多項式的多項式 Pm(x) 精確成立，即證明：精確成立，即證明： nkkmkbamxPAdxxPx0)()()( 設設)()()()(xrxqxwxPm bababamdxxrxdxxqxwxdxxPx)()()()()()()( 0 nkkkxrA0)( nkkmkxPA0)( x0 xn 為為 Gauss 點點 與任意次數(shù)與任意次數(shù)不大于不大于n 的多項式的多項式 P(x) （帶權）正交（帶權）正交。 nkkxxxw0)()(定理定理求求 Gauss 點點 求求w(x)的的零點零點一、高斯點的基本特性一、高斯點的基本特性 Gauss 公式的余項：公式

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

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

51、大有性質(zhì)：任意次數(shù)不大于于n 的多項式的多項式 P(x) 必與必與 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：

52、求求 2 = 0 的的 2 個根，即為個根，即為 Gauss 點點 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 注：注：構造正交多項式也可以利用最小二乘數(shù)據(jù)擬合中介構造正交多項式也可以

53、利用最小二乘數(shù)據(jù)擬合中介紹過的遞推式進行。紹過的遞推式進行。 特殊正交多項式族：特殊正交多項式族： Legendre 多項式族：多項式族：1)( x 定義在定義在 1, 1上，上，以以 Pn+1 的根為節(jié)點的求積公式稱為的根為節(jié)點的求積公式稱為Gauss-Legendre 公式公式。注：注：一般一般a,b上的積分可化為上的積分可化為-1,1上特殊高斯公式進行計算。上特殊高斯公式進行計算。例例 用用3點點Gauss公式計算積分公式計算積分 解解 查表得查表得x1=-0.7745966692,x2=0,x3=0.7745966692, A1=A3=0.5555555556,A2=0.8888888

54、889, 所以有所以有 Gauss-Legendre求積公式的余項為求積公式的余項為 ) 1 , 1(,)() 12()!2() !(2)2(3412nnfnnnfR11.cos xdxI68300355. 1coscoscos332211xAxAxAI誤差為 5347103492. 6)cos(7) ! 6(62fR實際上實際上,I=2sin1=1.68294197, 誤差為誤差為|R |=6.158 10-5 . 用用Simpson公式公式,則有則有I 1.69353487, 誤差為誤差為|R |=1.06 10-2 . 由于由于因此因此,a,b上權函數(shù)上權函數(shù) (x)=1的的Gauss型

55、求積公式為型求積公式為batabbaxdttabbafabdxxf)2)()()22(2)(11baniiixabbafAabdxxf1)22(2)(求積誤差可表示為求積誤差可表示為),(,)() 12()!2() !()()2(3412bafnnnabfRnn例例 用用3點點Gauss公式計算積分公式計算積分結果遠比結果遠比Simpson公式的結果精確公式的結果精確.102.14dxxI 解解 這里這里Gauss點和積分系數(shù)與上例相同點和積分系數(shù)與上例相同,所以所以 312112102141068. 32/ )1(1421)2121(142114iiixAdttdxxI Chebyshev

56、多項式族：多項式族：211)(xx 定義在定義在 1, 1上，上，) arccos( cos)(xkxTk Tn+1 的根為的根為 2212cosnkxkk = 0, , n以此為節(jié)點構造公式以此為節(jié)點構造公式 1102)()(11nkkkxfAdxxfx稱為稱為 Gauss-Chebyshev 公式公式。注意到積分端點注意到積分端點 1 可能是積分可能是積分的的奇點奇點，用普通，用普通Newton-Cotes公公式在端點會出問題。而式在端點會出問題。而Gauss公公式可能避免此問題的發(fā)生。式可能避免此問題的發(fā)生。Gauss-Laguerre求積公式為求積公式為 求積公式的誤差為求積公式的誤差

57、為 由于由于 所以所以,對對0, + )上權函數(shù)上權函數(shù) (x)=1的積分的積分,也可以構造類似的也可以構造類似的Gauss-Laguerre求積公式求積公式:01)()(niiixxfAdxxfe), 0(,)()!2() !()2(2nfnnfR00)()(dxxfeedxxfxx01)()(niixixfeAdxxfi4.6 4.6 數(shù)值微分數(shù)值微分 數(shù)值微分的數(shù)值微分的概念概念 數(shù)值微分的數(shù)值微分的計算方法計算方法 原始概念近似原始概念近似: :中點法及外推法中點法及外推法 函數(shù)近似函數(shù)近似: :插值型的求導公式插值型的求導公式 函數(shù)相互關系轉化函數(shù)相互關系轉化: :利用數(shù)值積分求導利

58、用數(shù)值積分求導 數(shù)值微分的數(shù)值微分的誤差分析誤差分析 泰勒展式估計泰勒展式估計 事后誤差估計事后誤差估計 基本關系轉化基本關系轉化 數(shù)值微分數(shù)值微分就是就是用函數(shù)值的線性組合近似函數(shù)在某點用函數(shù)值的線性組合近似函數(shù)在某點的導數(shù)值的導數(shù)值一、中點法和外推法一、中點法和外推法 按導數(shù)定義按導數(shù)定義 , 是差商是差商 當當 時的極限時的極限取取差商差商作為作為導數(shù)導數(shù)的近似值的近似值，建立簡單的數(shù)值微分方法，建立簡單的數(shù)值微分方法:)(0 xf hxfhxf)()(00 0h hxfhxfxf000 (8.1)向后差商近似導數(shù)向后差商近似導數(shù)(8.2)(8.3)中心差商近似導數(shù)中心差商近似導數(shù) hh

59、xfxfxf 000 hhxfhxfxf2000 hhxfhxfhD200 容易看出，就精度而言，以（容易看出，就精度而言，以（8.3）式更為可取，稱）式更為可取，稱(8.4)為為 的的中點公式中點公式, 其中其中h為一增量，稱為為一增量，稱為步長步長 這種數(shù)值這種數(shù)值微分方法稱為微分方法稱為中點方法中點方法, 它是前兩種方法的算術平均它是前兩種方法的算術平均)(0 xf 分別將分別將在在 x=a 處做處做Taylor展開有展開有)(haf )(! 5)(! 4)(! 3)(! 2)()()() 5(5) 4(432afhafhafhafhafhafhaf代入代入D(h)得得 )(! 5)(!

60、 3)()()5(42afhafhafhD,6)()(2MhhDaf 其中其中)(maxxfMhax 現(xiàn)在來考慮中點公式現(xiàn)在來考慮中點公式 的截斷誤差，的截斷誤差，hhafhafhD2)()()( (8.5)所以所以截斷誤差截斷誤差 )(! 5)(! 3)()() 5(42afhafhafhD(8.6)從截斷誤差的角度來看，步長從截斷誤差的角度來看，步長h越小，計算結果越準確。且越小，計算結果越準確。且 所以所以, 在在實際計算時實際計算時，通常，通常采用采用二分步長二分步長及及誤差事后估誤差事后估計法計法, 在變步長的過程中實現(xiàn)步長的自動選擇，在保證截斷在變步長的過程中實現(xiàn)步長的自動選擇，在