1、第5章 微積分的數(shù)值解法 在現(xiàn)代工程領(lǐng)域或者科研過(guò)程中所遇到的許多實(shí)際問(wèn)題的求解，常常歸結(jié)為定積分的計(jì)算。例如，力學(xué)和電氣中功和功率的計(jì)算、電流電壓平均值及有效值的計(jì)算，幾何學(xué)中面積和何種及其重心的計(jì)算，化學(xué)化工中的熱容、熱焓、轉(zhuǎn)化率、反應(yīng)時(shí)間和反應(yīng)器體積等等的計(jì)算，都與某些定積分的計(jì)算有關(guān)。但很多情況下，其積分函數(shù)f(x)的關(guān)系難以明確表示，即使能用解析明確表示，有時(shí)表達(dá)式也很復(fù)雜，不能用于實(shí)際計(jì)算。本章介紹計(jì)算積分和導(dǎo)數(shù)的實(shí)用數(shù)值方法。 5.1 數(shù)值積分的基本思想 一、定積分的數(shù)學(xué)含義 由定積分的定義： 其中 為第i個(gè)積分區(qū)間 內(nèi)任一點(diǎn)。當(dāng)n足夠大時(shí)有以下近似表達(dá)式：將被積函數(shù)在積分區(qū)間上

2、足夠多的離散點(diǎn)的函數(shù)值與其所在小區(qū)間長(zhǎng)度乘積之和作為定積分的近似數(shù)值。 1( )lim( )bniiniaIf x dxfx(5-1)iix1( )( )bniiiaIf x dxfx(5-2)5.1 數(shù)值積分的基本思想二、定積分的幾何含義 （5-2）式告訴我們，定積分的值就是由下面四條曲線所轉(zhuǎn)面的面積，即： ,0,( )xa xb yyf xab( )f x( )f圖5-1 定積分的幾何意義 5.1 數(shù)值積分的基本思想如圖5.1中所示的陰影部分的面積。這就是定積分的幾何意義。 當(dāng)求積區(qū)間只有一個(gè)時(shí)當(dāng)求積區(qū)間只有一個(gè)時(shí)1. 由積分中值定理，在（a, b）中必存在一點(diǎn)，使得 2.用f(a)和f(

3、b)的平均值近似代替f()，則：3.用中點(diǎn)的函數(shù)值近似代替f()，則： I( )() ( )baf x dxba f( )( )I( )()2baf af bf x dxbaI( )()2baabf x dxba f5.1 數(shù)值積分的基本思想推廣到更一般情況，當(dāng)把求積區(qū)間(a,b)(a,b)分成n個(gè)小求積區(qū)間時(shí)，則根據(jù)定積分的幾何定義，我們有其中 為小區(qū)間長(zhǎng)度分點(diǎn)為 ，它們所對(duì)應(yīng)的縱坐標(biāo)（函數(shù)）分別為 。如圖5-2，取其中任一個(gè)小區(qū)間 進(jìn)行分析。該區(qū)間為一曲邊梯形，y用以下方法求其面積。1I( )( )bniiiaf x dxfxibaxxn ix0121,kknxa x xxxxb121(

4、),( ),(),(),(),( )kkf af xf xf xf xf b1(,)(0,1,2,1)kkxxkn5.1 數(shù)值積分的基本思想1. 用小矩形或小梯形近似表示時(shí)，然后求和得到I , 這就是高等數(shù)學(xué)中的矩形積分法和梯形積分法。它們都是用直線（一次函數(shù)）代替曲線（f(x)函數(shù)）求得的近似值，精度較低，一般要求N取得足夠大才能保證精度。2. 當(dāng)用曲線去逼近f(x)時(shí)，只要曲線近似程度足夠高，一般都能滿足計(jì)算精度要求。 本章按照上述思想，主要介紹梯形積分法、拋物線積分法（Simpson法)、Newton-Cotes積分法，龍貝格積分法以及高斯積分法5.2 梯形積分法 一 定步長(zhǎng)梯形積分法1

5、. 算法分析 設(shè)在區(qū)間(a, b)上有可積函數(shù)f(x),求積分值將區(qū)間(a, b)分成n個(gè)相等的小區(qū)間，每個(gè)小區(qū)間之長(zhǎng)為： h=(b-a)/n各分點(diǎn)： x0=a, x1, x2, , xi, xi+1, , xn=b其中 xi=a+i*h (i=1,2,3, ,n) 用p(x)=c*x+d直線近似代替f(x)(參見圖5-2） T( )baf x dx5.2 梯形積分法用插值的方法，我們可求得將其代入積分公式有11(,()iixf x( ,( )iixf x( )p xcxd1111( )()( )iiiiiiiixxxxp xf xf xxxxx5.2 梯形積分法第i個(gè)區(qū)間： Ti=(f(xi

6、)+f(xi+1)*h/2 第i-1個(gè)區(qū)間： Ti-1=(f(xi-1)+f(xi)*h/2 求和得到： 定步長(zhǎng)梯形積分法的截?cái)嗾`差為: 1101T( ( )( )/2( ) ,1,2,3, -1nniiiihf af bf xxaxxhin(5-3)(1)02( ) ( )( )d()d(1)!()( )12nnbbnkaakfR ff xpxxxxxnbah f5.2 梯形積分法 該誤差按h的2次方的速度下降,即定步長(zhǎng)梯形積分法的誤差階為2，也就是說(shuō)，定步長(zhǎng)梯形積分法具有一階代數(shù)精度（因?yàn)榇鷶?shù)精度等于誤差階減1，所以2-1=1）。2. 定步長(zhǎng)梯形積分法的程序框圖與通用程序設(shè)計(jì) 定步長(zhǎng)梯形積

7、分法的通用程序框圖如圖5-3所示，共分為三個(gè)程序框圖，即在圖中（1）為調(diào)用子程序的主控程序框圖，（2）為函數(shù)子程序框圖，（3）為積分通用子程序框圖。5.2 梯形積分法（1）主控程序框圖（2）函數(shù)子程序框圖STARTINPUT a,b,nCALL sub TOUTPUT solENDFUNCTION FF=exprEND FUNCTION圖5-3 梯形積分法的程序框圖5.2 梯形積分法subroutineh=(b-a)/nComp f(a),f(b)T=(f(a)+f(b)/2DO i=1,n-1x=a+i*hComp f(x)T=T+f(x)END DO iT=T*hEND subroutin

8、e（3）梯形積分子程序框圖5.2 梯形積分法10C 主控程序20 PROGRAM main30 read(5,*) a,b,n40 call SUBRPUTINE50 &txjf(a,b,n,T) 60 write(6,*) T70 END80C 函數(shù)子程序90 FUNCTION f(x)100 f=expr(x)110 END FUNCTION f120C 梯形積分子程序130 SUBRPUTINEtxjf(a,b,n,T) 140 h=(b-a)/n150 T=(f(a)+f(b)/2160 x=a170 DO i=1,n-1180 x=a+i*h(x=x+h)190 T=T+f(

9、x)200 END DO210 T=T*h220 END SUBRPUTINE txjf 11( ) niif x5.2 梯形積分法二、變步長(zhǎng)梯形積分法 前面介紹了定步長(zhǎng)梯形積分法。采用定步長(zhǎng)梯形積分法求定積分時(shí)，步長(zhǎng)的選擇應(yīng)該是適當(dāng)?shù)?。如果采用定步長(zhǎng)梯形積分法，只有改變N值，比較幾組計(jì)算結(jié)果才能判是否達(dá)到精度要求，這樣程序的執(zhí)行就受到人工干擾，占用了機(jī)時(shí)，降低了計(jì)算機(jī)的效率。為此，需要尋求新的改善計(jì)算精度途徑。截?cái)嗾`差分析表明，只要步長(zhǎng)h充分小，必可滿足精度要求。因此，通常采取逐步縮小步長(zhǎng)h的辦法。自動(dòng)變步長(zhǎng)梯形積分法是計(jì)算機(jī)根據(jù)精度要求進(jìn)行自動(dòng)判斷，逐步逐步縮小步長(zhǎng)h。1. 算法分析 其基

10、本思想為：逐步變更步長(zhǎng)，用二分法使步長(zhǎng)逐次變小， 直到滿足精度。 5.2 梯形積分法N=1:N=1: h1=b-a T1=f(a)+f(b)*h1/2 N=2:N=2: h2=h1/2 T2=f(a)+f(b)*h2/2+f(x1)*h2 =T1/2+f(x1)*h2 N=4: N=4: h4=h2/2T4=f(a)+f(b)*h4/2+f(x1)+f(x2)+f(x3)*h4 =T2/2+f(x1)+f(x3)*h4 abab1xab1x2x3x5.2 梯形積分法N=2n:N=2n:注：注：xi=a+(2*i-1)*h2n：上一次積分區(qū)間的中點(diǎn)。 在計(jì)算 的過(guò)程中，每當(dāng)算出一新的近似值時(shí)，便

11、檢查判斷下列條件是否成立： 如果條件成立，則 就是符合精度要求的積分近似值，否則縮小步長(zhǎng)h2n再行計(jì)算，直到滿足條件為止。 22221 /2 (2 -1) TT /2( ) nninnnnniihhxaihhf x12482,nnT T T TT T2nnTTE （E為容許誤差） 2nT(5-4)5.2 梯形積分法變步長(zhǎng)梯形積分法的截?cái)嗾`差為:2. 變步長(zhǎng)梯形積分法的程序框圖與通用程序設(shè)計(jì) 變步長(zhǎng)梯形積分法的通用程序框圖如圖5-4所示，也共分為三個(gè)程序框圖，即在圖中（1）為調(diào)用子程序的主控程序框圖，（2）為函數(shù)子程序框圖，（3）為積分通用子程序框圖。（1）和（2）與上一節(jié)雷同，因此這里只畫出梯

12、形變步長(zhǎng)積分法通用子程序框圖（3）。2() ( )12 4nbaR fh f 5.2 梯形積分法subroutineh=(b-a)Comp f(a),f(b)T1=(f(a)+f(b)/2DO i=1,nh=h/2;T2=T1/2Comp f(x)T2=T2+f(x)END DO in=2*n;T1=T2END subroutine圖5-4 梯形積分法子程序框圖n=1x=a+(2*i-1)*habs(T2-T1)e) then100 n=2*n;T1=T2110 goto 45120 endif130 END SUBRPUTINE btxjf5.2 梯形積分法3.應(yīng)用實(shí)例分析例例5.15.1

13、已知圓周率的近似值可以由下列定積分求出試求出圓周率的似近值。解：解： 因?yàn)樗郧蠼獯藛?wèn)題的函數(shù)程序程序?yàn)?0C 函數(shù)子程序90 FUNCTION f(x)100 f=4.0/(1+x*x)110 END FUNCTION f 12041dxx24( )1f xx5.2 梯形積分法n n101020204040100100200200300300T T33139926331399263.1411763.1411763.1414893.1414893.1415763.1415763.1415893.1415893 3141591141591表5.1 梯形積分法輸入不同N地的計(jì)算機(jī)輸出結(jié)果 由上表可

14、以看出，定步長(zhǎng)梯形積分法的收斂速度是很慢的，由上表可以看出，定步長(zhǎng)梯形積分法的收斂速度是很慢的，對(duì)于本例，對(duì)于本例，N=300N=300時(shí)，即把積分區(qū)間時(shí)，即把積分區(qū)間00，11劃分為劃分為300300個(gè)小區(qū)個(gè)小區(qū)間時(shí)才使計(jì)算結(jié)果精確到間時(shí)才使計(jì)算結(jié)果精確到 0.000010.00001。 5.3 Simpson積分法 梯形積分法是用直線段（一次二項(xiàng)式）去逼近曲線f(x),而Simpson 是用拋物線段（二次三項(xiàng)式）去逼近曲線f(x). 一、定步長(zhǎng)Simpson積分法1.算法分析 設(shè)將（a, b）分為n個(gè)（偶數(shù)）相等的小區(qū)間，則 求積節(jié)點(diǎn)坐標(biāo) 對(duì)應(yīng)求積點(diǎn)函數(shù)值( )baSf x dx( - )

15、/hb an012, , , nxa xxxb011( ), ( ),( )nyf ayf xyf b5.3 Simpson積分法被積曲線用 y=f(x) 表達(dá)，而逼近曲線則采用 y=Ax2+Bx+C 0 xa1x2xix1ix2ixnxb( )yf x2yAxBxC5.3 Simpson積分法第一段：x0,x1,x2 有同理，第二段：x2,x3,x4 有第末段： xn-2,xn-1,xn 有200021112222yAxBxCyAxBxCyAxBxC求出系數(shù)求出系數(shù)A, B, C 2021012 ()d(4)/3xxSAxBxCxyyyh2234(4)/3Syyyh/221(4)/3nnnn

16、Syyyh5.3 Simpson積分法共有n/2小塊面積，對(duì)它們求和即可得Simpson法求積公式或這里01234-1013-124-2(42424)/3()4()2()/3nnnnnSyyyyyyyhyyyyyyyyh(5-5)/2/2 121211/2/221211 ( )( )4()2() /3 ( )-( )4()2() /3nnijijnnijijSf af bf xf xhf af bf xf xh(5-6)212(21),2iixaihxaih5.3 Simpson積分法同理我們可以推出辛普森積分法的截?cái)嗾`差為：由上式可知辛普森積分法的誤差階和代數(shù)精度分別為4和3，在積分小區(qū)間總

17、數(shù)與梯形積分法相同，辛普森積分法的精度較高。2.通用程序框圖與通用程序設(shè)計(jì) 這里只畫出Simpson定步長(zhǎng)積分法通用子程序框圖（3）。 4(4)() ( )180baR fh f (5-7)5.3 Simpson積分法subroutineh=(b-a)/nComp f(a),f(b)S=f(a)-f(b)DO i=1,mComp f(x1), f(x2)S=S+4*f(x1)+ 2*f(x2)END DO ix2=a+(2*i)*hEND subroutine圖5-5 Simpson法子程序框圖m=n/2x1=a+(2*i-1)*hS=S*h/35.3 Simpson積分法10 SUBRPUT

18、INE Simpson(a,b,n)10 SUBRPUTINE Simpson(a,b,n)20 h=(b-a)/n20 h=(b-a)/n30 m=n/230 m=n/240 S=f(a)-f(b)40 S=f(a)-f(b)50 DO i=1,m50 DO i=1,m60 x1=a+(260 x1=a+(2* *i-1)i-1)* *h h70 x2=a+270 x2=a+2* *i i* *h h80 S=S+4.080 S=S+4.0* *f(x1)+2.0f(x1)+2.0* *f(x2)f(x2)90 END DO90 END DO100 S=S100 S=S* *h/3.0h/3

19、.0110 END SUBROUTINE Simpson110 END SUBROUTINE Simpson5.3 Simpson積分法二、變步長(zhǎng)Simpson積分法1.算法分析 步長(zhǎng)縮小和誤差控制方法均與變步長(zhǎng)積分雷同，即步長(zhǎng)折半、 ?；舅枷肴缦拢篘=1N=1： 求積點(diǎn) x0=a,x1,x2=b H1=(b-a)/2 S1=f(a)+4*f(x1)+f(b)*H1/3N=2N=2： 求積點(diǎn) x0=a,x1,x2,x3,x4=b H2=H1/2 S2=f(a)+4f(x2)+2f(x1)+4f(x3)+f(b)*H2/3 =f(a) +f(b)+2*f(x2) +4*f(x1)+f(x3)*

20、H2/3 2 |-|nnSSab1xab1x2x3x5.3 Simpson積分法N=2n:N=2n:關(guān)于變步長(zhǎng)Simpson積分法的遞推計(jì)算公式以及程序框圖和通用FORTRAN程序?qū)⒃诤竺鎺追N積分比較部分提示。這里給出一個(gè)通用的FORTRAN程序清單： 12212211 ( )( )4()2 () /3nnniiniiSf af bf xf xh(5-8)5.3 Simpson積分法10 SUBRPUTINE Simpb(a,b,E)10 SUBRPUTINE Simpb(a,b,E)20 h=(b-a)/220 h=(b-a)/230 s2=0;n=130 s2=0;n=140 s0=f(a

21、)+f(b)40 s0=f(a)+f(b)45 s1=f(a+h)45 s1=f(a+h)50 s=(s0+4.050 s=(s0+4.0* *s1)s1)* *h/3.0 h/3.0 6060 n=2 n=2* *n;h=h/2.0n;h=h/2.070 s2=s2+s170 s2=s2+s180 s1=0.080 s1=0.090 x=a+h90 x=a+h100 DO i=1,n100 DO i=1,n110 s1=s1+f(x)110 s1=s1+f(x)120 x=x+h+h130 END DO140 s2n=(s0+2.0*s2+4.0*s1)150&*h/3.0160 i

22、f(abs(s2n-s)E) then170 s=s2n180 goto 60190 endif200 END SUBROUTINE Simpb 5.4 Newton-Cotes積分法 在梯形積分法中，是用一次二項(xiàng)式（直線）去逼近曲線f(x);在Simpson積分法中，是用二次三項(xiàng)式（拋物線）去逼近曲線f(x)。而在N-C積分法中，用的是m 次多項(xiàng)式去逼近曲線f(x)。 設(shè)在被積區(qū)間a, b選取m+1個(gè)插值點(diǎn)，用m 次多項(xiàng)式pm去逼近曲線f(x),則：利用拉格郎日插值多項(xiàng)式 ( )( )bbmaaf x dxpx dx()0( )()(),0,1,2, bmmkkkaf x dxbaCf xk

23、m(5-9)5.4 Newton-Cotes積分法其中:Ck (m)為Cotes系數(shù)，可以具體求出起其值：m=1: C0(1) =1/2; C1(1)=1/2m=2: C0(2)=1/6; C1(2)=4/6; C2(2)=1/6見下表 注：在實(shí)際應(yīng)用中，一般m9，否則會(huì)產(chǎn)生計(jì)算不穩(wěn)定。 上面所談到的問(wèn)題，是在整個(gè)區(qū)間a, b上用Pm去逼近f(x)。下面看看用N-C公式即(5-1)式用在小區(qū)間時(shí)情況，將a, b分成n個(gè)小區(qū)間。有： ()0( 1)(1)(1)(1)()! ()!mm kmkCt ttktktm dtm kmk5.4 Newton-Cotes積分法 K m 0 0 1 1 2 2

24、 3 3 4 4 5 5 6 6 7 7 8 8 11/21/21/21/2 21/61/64/64/61/61/6 31/61/63/83/83/83/81/61/6 47/907/9032/9032/9012/9012/9032/9032/907/907/90 519/28819/28825/9625/9625/14425/14425/14425/14425/9625/9619/28819/288 641/84041/8409/359/359/2809/28034/10534/1059/2809/2809/359/3541/84041/840 7751/17751/172802803577

25、/173577/172802801323/171323/172802802989/172298917229891721323173577172751/1728080 8989/28989/283503505838/285838/28350350- -928/283928/283505010496/2810496/28350350- -4540/2834540/283505010496/2810496/28350350- -928/283928/28350505838/285838/28350

26、350989/283989/28350505.4 Newton-Cotes積分法 x0=a, x1, x2, , xn-1, xn=b, h=(b-a)/n當(dāng)項(xiàng)次m=1時(shí)： x0, x1 x1, x2xn-1,xn 1(1)(1)10001110( )()()() ( )() /2xxf x dxxxCf xCf xf af xh2121( )( ()() /2xxf x dxf xf xh11( )( ()( ) /2nnxnxf x dxf xf b h5.4 Newton-Cotes積分法求和得到 *復(fù)化梯形求積公式同理，當(dāng)m=2時(shí)，在每一個(gè)小區(qū)間xk, xk+1上增加一個(gè)分點(diǎn)（中點(diǎn)），

27、xk+1/2=（xk+ xk+1）/2= xk+ h/2，此時(shí)區(qū)間a,b上共有2n+1個(gè)等距求積節(jié)點(diǎn)，即： x0, x0+1/2, x1 , x1+1/2, x2, , xn+1/2, xn 兩節(jié)點(diǎn)之間的距離（新的步長(zhǎng)）h=h/2。因此， m=2時(shí)的求積公式為1111( ) ( )( )2( ) /2 ( )( )2( )()/2bniianiif x dxf af bf xhf af bf xban(5-10)5.4 Newton-Cotes積分法*復(fù)化Simpson求積公式求和當(dāng)項(xiàng)次m=4時(shí): 在xk, xk+1上增加三個(gè)分點(diǎn)xk+1/4,xk+2/4, xk+3/4 使xk+1/4 =

28、xk+ h/4， xk+2/4 = xk+ 2h/4， xk+3/4 = xk+ 3h/4 ，此時(shí)區(qū)間a,b上共有4n+1個(gè)等距求積節(jié)點(diǎn)，兩節(jié)點(diǎn)之間的距離（新的步長(zhǎng)）h=h/4。同理，求和得111/410112/43/400( )7 ( ) ( ) 14( )32()+12()32()()/90bnnikikannkkkkf x dxf af bf xf xf xf xban(5-12)111/201121211( ) ( )( )4()2( ) /6 ( )( )4()2() ()/6bnnkikianniiiif x dxf af bf xf xhf af bf xf xban(5-11)

29、5.5 RombergRomberg積分法*復(fù)化Newton-Cotes求積公式一、一、 RombergRomberg積分積分算法分析算法分析 RombergRomberg積分法積分法是是用計(jì)算精度較低的二分前后兩次積分近似值，用計(jì)算精度較低的二分前后兩次積分近似值，通過(guò)誤差補(bǔ)償法來(lái)獲得收斂速度較快的一種積分方法。通過(guò)誤差補(bǔ)償法來(lái)獲得收斂速度較快的一種積分方法。1 1、 梯形積分法的誤差估計(jì)（見式梯形積分法的誤差估計(jì)（見式3-33-3）( (線性組合線性組合) ) T T* *-T-T2n2n(T(T2n2n-T-Tn n)/3)/3 其中其中T T* *為精確值為精確值 T T* * T T

30、2n2n +(T(T2n2n-T-Tn n)/3)/3 4 4* *T T2n2n/3-T/3-Tn n/3 /3 =S =Sn n (5-13) (5-13)式說(shuō)明式說(shuō)明 T T2n2n的誤差約等于的誤差約等于（T T2n2n-T-Tn n）/3/3, ,將這個(gè)誤差作為將這個(gè)誤差作為一種補(bǔ)償加到一種補(bǔ)償加到T T2n2n上，可期望得到更精確的上，可期望得到更精確的SimpsonSimpson的積分公式的積分公式 它是它是梯形積分法二分前后兩次積分近似值進(jìn)行線性組合。梯形積分法二分前后兩次積分近似值進(jìn)行線性組合。 (5-13)(5-14)5.5 RombergRomberg積分法2 2、Sim

31、psonSimpson積分法的誤差估計(jì)積分法的誤差估計(jì)( (線性組合線性組合) ) S S* *-S-S2n2n(S(S2n2n-S-Sn n)/15)/15 其中其中S S* *為精確值為精確值 S S* *S S2n2n +(S(S2n2n-S-Sn n)/15 )/15 16 16* *S S2n2n/15-S/15-Sn n/15 /15 =C =Cn nSimpsonSimpson積分法二分前后兩次積分近似值進(jìn)行線性組合，得到積分法二分前后兩次積分近似值進(jìn)行線性組合，得到Newton-CotesNewton-Cotes的積分公式（精度更高）。的積分公式（精度更高）。3 3、Cotes

32、Cotes積分公式進(jìn)行線性加速積分公式進(jìn)行線性加速C C* *-C-C2n2n(C(C2n2n-C-Cn n)/63)/63其中其中C C* *為精確值。為精確值。同理同理 C C* * C C2n2n +(C(C2n2n-C-Cn n)/63 )/63 64 64* *C C2n2n/63-C/63-Cn n/63 /63 =R =Rn n(5-15)(5-16)(5-17)(5-18)5.5 RombergRomberg積分法 Cotes Cotes積分法二分前后兩次積分近似值進(jìn)行線性組合，得到積分法二分前后兩次積分近似值進(jìn)行線性組合，得到RombergRomberg的積分公式（精度更高）

33、。注意，此時(shí)不能再對(duì)的積分公式（精度更高）。注意，此時(shí)不能再對(duì)RombergRomberg的積分進(jìn)行線性加速。因?yàn)榇藭r(shí)插值多項(xiàng)式已經(jīng)為八的積分進(jìn)行線性加速。因?yàn)榇藭r(shí)插值多項(xiàng)式已經(jīng)為八次，它具有九次插值多項(xiàng)式的精度。依據(jù)插值多項(xiàng)式的理論，次，它具有九次插值多項(xiàng)式的精度。依據(jù)插值多項(xiàng)式的理論，當(dāng)插值多項(xiàng)式的次數(shù)超過(guò)九次時(shí)，將出現(xiàn)不穩(wěn)定情況。因此，當(dāng)插值多項(xiàng)式的次數(shù)超過(guò)九次時(shí)，將出現(xiàn)不穩(wěn)定情況。因此，到此為止。到此為止。 RombergRomberg積分法不僅精度高，而且它從最簡(jiǎn)單的梯形積分開始，積分法不僅精度高，而且它從最簡(jiǎn)單的梯形積分開始，逐步加工成收斂速度快的逐步加工成收斂速度快的CotesC

34、otes積分法。積分法。二、二、 RombergRomberg積分法積分法通用程序框圖與通用程序設(shè)計(jì)通用程序框圖與通用程序設(shè)計(jì)1.1.算法過(guò)程算法過(guò)程5.5 RombergRomberg積分法T1T2T4T8T16S1S2S4S8C1C2C3R1R2圖5-6 Romberg法子積分圖end subroutineh=b-a;t1=(f(a)+f(b)*h/2;K=1s=0; x=a+h/2s=s+f(x); x=x+hxb?t2=t1/2+s*h/2s2=t2+(t2-t1)/3k=1?k=k+1;h=h/2;t1=t2;s1=s2c2=s2+(s2-s1)/15k=2?r2=c2+(c2-c1

35、)/63c1=c2k=3?r1=r2abs(r2-r1)E?YNNNNYNYYsubroutine圖5-7 Romberg法子程序框圖5.5 RombergRomberg積分法10 subroutine romberg(a,b,E)10 subroutine romberg(a,b,E)20 h=b-a;t1=(f(a)+f(b)20 h=b-a;t1=(f(a)+f(b)* *h/2h/230 k=130 k=140 s=0;x=a+h/240 s=0;x=a+h/250 if(xb) then 50 if(xE) then220 if(abs(r2-r1)E) then230 goto 1

36、90 230 goto 190 240 endif240 endif250 end subroutine romberg250 end subroutine romberg3.3.程序語(yǔ)言程序語(yǔ)言5.5 RombergRomberg積分法 讀懂本程序的關(guān)鍵在于看懂控制變量k(步長(zhǎng)折減次數(shù)）。當(dāng)k=1時(shí)，程序剛開始，程序只求出t2和s2(實(shí)際上s1)；當(dāng)k=2時(shí)，步長(zhǎng)折二分，t1=t2, s1=s2,且程序又求出新的t2(實(shí)際上為t4)和s2，同時(shí)求出c2，并且令c1=c2；當(dāng)k=3時(shí)，步長(zhǎng)再折二分，t1=t2, s1=s2, 同時(shí)又求出新的t2(實(shí)際上為t8)和s2(實(shí)際上為s4)。求出新的c

37、2,且執(zhí)行c1=c2,求出r2,且執(zhí)行r1=r2, c1=c2；當(dāng)k=4時(shí)，步長(zhǎng)再折二分， t1=t2, s1=s2, 同時(shí)又求出新的t2(實(shí)際上為t16)和s2(實(shí)際上為s8)以及c2(實(shí)際上為c4)，新的r2。此后，只要當(dāng)k3，程序每次都執(zhí)行判斷abs(r2-r1)是否滿足精度要求，若滿足，則結(jié)束，否則執(zhí)行r2=r1和c1=c2，且重復(fù)k=4（只是k值不同）的過(guò)程，直至滿足精度要求。5.6 GaussGauss積分法 前面介紹的幾種積分法都是在等分區(qū)間情況下進(jìn)行討論的，而Gauss積分法則是適用于不等分區(qū)間的一種非常有效的數(shù)值積分方法。一、高斯點(diǎn)及高斯公式簡(jiǎn)介 首先介紹一下代數(shù)精度的概念。

38、一個(gè)求積公式，如果對(duì)于任何一個(gè)不超過(guò)m次多項(xiàng)式fm(x)均能使其余項(xiàng)積分為零，即Rf=0，但對(duì)于m+1次多項(xiàng)式不能能使其余項(xiàng)積分為零，即Rf0,則稱此求積公式具有m次代數(shù)精度。 內(nèi)插求積公式可以表示為0( )d( )nbiiaif xxA f x(5-19)5.6 GaussGauss積分法 在等距節(jié)點(diǎn)的情況下，我們可以驗(yàn)證它具有n次代數(shù)精度。假定若節(jié)點(diǎn)數(shù)目一定，我們能否通過(guò)適當(dāng)?shù)剡x擇節(jié)點(diǎn)的位置，使積分公式具有更高的代數(shù)精度呢？Gauss首先指出，若我們選擇的接點(diǎn)x0,x1,x2, xn所構(gòu)成的多項(xiàng)式 與任意次數(shù)不超過(guò)n次的多項(xiàng)式p(x)都正交，即 則因 為次數(shù)不超過(guò)2n+1次的多項(xiàng)式，故求積

39、公式(5-19)式具有2n+1次代數(shù)精度。 考慮積分012( )()()()()nxxxxxxxxx( ) ( )d0bax p xx( ) ( )x p x( )dbaf xx5.6 GaussGauss積分法 的上下限a,b是變動(dòng)的，通過(guò)積分變換可將上下限變?yōu)楣潭ǖ闹?1和1。事實(shí)上，令 ，可以將積分區(qū)間a,b變換為-1,1，即因此，不失一般性，可考察區(qū)間-1,1上的Gauss積分公式 () /2()/2xba tba11( )d()d222babababaf xxftt110( )d( )niiif xxA f x(5-20)5.6 GaussGauss積分法常用的Gauss積分公式如下

40、：（1）Gauss-Legendre積分公式 （4）Gauss-Hermite積分公式（2）Gauss-Chebyshev積分公式（3）Gauss-Laguerre積分公式110( )d( )niiif xxA f x11021( )d(cos)122niiif xxA fnn00( )d( )nxiiief xxA f x20( )d( )nxiiief xxA f x5.6 GaussGauss積分法nxiAinxiAi00230.86113630.33998100.34785480.652145211208/95/9400.90617980.53846930.56888890.23692

41、690.47862873/53/3表Gauss-Legendre求積公式的xi與Ai5.6 GaussGauss積分法nxiAinxiAi01130.32254771.74576114.53662039.35907090.60315410.35741870.03888790.000539310.58578643.41421360.85355340.146446620.41577462.29428046.28994510.71109300.27851770.010389340.26356031.41340313.59642587.085810012.64080080.52175560.39866

42、680.07594240.00361180.0000234表2 Gauss-Laguerre求積公式的xi與Ai5.6 GaussGauss積分法nxiAinxiAi001.772453930.52464761.65068010.80491410.081312810.70710680.8862269201.22474491.18163590.2954110400.95857262.02018290.94530870.39361630.0199532表3 Gauss-Hermite求積公式的xi與Ai5.6 GaussGauss積分法例用階數(shù)n=2的Gauss求積公式計(jì)算以下積分的近似值解：由于

43、n=2，在表中查出Gauss求積點(diǎn)坐標(biāo)和權(quán)系數(shù)1204Id1xx01013/3,3/3,1xxAA 001122I()()22222144112131131112322323.14754086bababababaA fxA fx 5.6 GaussGauss積分法subroutinexm=(b+a)/2; xr= (b-a)/2; s=0DO i=1,ndx=xr*x(i)Comp f(xm+dx), f(xm-dx)s=s+w(i)*(f(xm+dx)+f(xm-dx)END DO is=xr*sEND subroutineINPUT x(i), w(i)圖5-8 Gauss積分子程序框圖5

44、.6 GaussGauss積分法10 subroutine gausint(a,b,s)10 subroutine gausint(a,b,s)20 dimension x(i),w(i)20 dimension x(i),w(i)30 data w/0.2955242,0.269266730 data w/0.2955242,0.269266740 &,0.2190864,0.149451340 &,0.2190864,0.149451350 &,0.066671350 &,0.066671360 data x/0.1488743,0.433395460 da

45、ta x/0.1488743,0.433395470 &,0.6794096,0.865063470 &,0.6794096,0.865063480 &,0.973906580 &,0.973906590 xm=(b+a)/2;xr=(b-a)/290 xm=(b+a)/2;xr=(b-a)/2100 s=0100 s=0110 do i=1,5110 do i=1,5120 dx=xr120 dx=xr* *x(i)x(i)130 s=s+w(i)130 s=s+w(i)* *(f(xm+dx)+(f(xm+dx)+140 & f(xm+dx) 140

46、 & f(xm+dx) 150 end do150 end do160 s=xr160 s=xr* *s s170 end subroutine gausint170 end subroutine gausint5.7 數(shù)值微分?jǐn)?shù)值微分法 在現(xiàn)代工程技術(shù)領(lǐng)域和科研過(guò)程中，數(shù)值微分法雖然不像數(shù)值積分法那樣用途廣泛，但還是經(jīng)常碰到數(shù)值微分問(wèn)題。嚴(yán)格地說(shuō)，數(shù)值微分法應(yīng)叫做數(shù)值導(dǎo)數(shù)法（微分是差，而導(dǎo)數(shù)才是差商的極限形式），但由于歷史的原因，我們這里仍延用數(shù)值微分這一概念。一、差商法求數(shù)值微分1、算法分析 考慮一個(gè)在點(diǎn)附近能夠解析的函數(shù)f(x)，在 點(diǎn)它可以展開為Taylor級(jí)數(shù)，并去前兩項(xiàng)，有

47、 ()ixh2()( )( )()iiif xhf xhfxO h(5-21)5.7 數(shù)值微分?jǐn)?shù)值微分法其中 為 的同階無(wú)窮小。由（5.21）式求得函數(shù)的一階導(dǎo)數(shù)表達(dá)示為我們稱（5.22）式為f(x)在 點(diǎn)的誤差階為h的向前差分式，即同樣使用泰勒級(jí)數(shù)將f(x)在 點(diǎn)展開并整理可得我們稱（5.24）式為f(x)在 點(diǎn)的誤差階為h的向后差分式2()O h2h()( )( )( )iiif xhf xfxO hh(5-22)()( )( )iiif xhf xfxh(5-23)()ixh( )()( )iiif xf xhfxh(5-24)ixix5.7 數(shù)值微分?jǐn)?shù)值微分法?。?.23）式和（5.2

48、4）式的算術(shù)平均值，有（5.25）式稱為f(x)在 點(diǎn)一階中心差商近似公式。2、程序框圖與程序設(shè)計(jì) ()()( )2iiif xhf xhfxh(5-25)ix5.7 數(shù)值微分?jǐn)?shù)值微分法subroutinex1=x-h;x2=x+hDO x=a,b,hcomp f(x1), f(x2),f(x)END DO xEND subroutineddf1=(f(x2)-f(x)/hddf2=(f(x)-f(x1)/hddf3=(f(x2)-f(x1)/2/h圖5-9 差商數(shù)值微分子程序框圖5.7 數(shù)值微分?jǐn)?shù)值微分法10 subroutine diffder(a,b,h)20 do x=a,b,h30 x1=x-h;x2=x+h40 ddf1=(f(x2)-f(x)/h50 ddf2=(f(x)-f(x1)/h60 ddf3=(f(x2)-f(x1)/h70 end do80 End subroutine diffder5.7 數(shù)值微分?jǐn)?shù)值微分法二、插值法求數(shù)值微分1、算法分析 用插值公式求數(shù)值微分有兩點(diǎn)公式、三點(diǎn)公式、四點(diǎn)公式和五點(diǎn)公式等等，下面介紹常用的三點(diǎn)公式。 已知f(x)在等距節(jié)點(diǎn) 的函數(shù)值 ，且函數(shù)f(x)的二階導(dǎo)數(shù)存在，作三點(diǎn)二次插值函數(shù)F(x)并對(duì)F(x)求導(dǎo)得三