Chapt5數值微積分的數值解法ppt課件

《Chapt5數值微積分的數值解法ppt課件》由會員分享，可在線閱讀，更多相關《Chapt5數值微積分的數值解法ppt課件（69頁珍藏版）》請在裝配圖網上搜索。

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

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

3、和f(b)的平均值近似代替的平均值近似代替f()，那么：，那么：3.用中點的函數值近似代替用中點的函數值近似代替f()，那么：，那么：I()()()baf x dxba f()()I()()2baf af bf x dxbaI()()2baabf x dxba f5.1 數值積分的基本思想推廣到更一般情況，當把求積區間(a,b)分成n個小求積區間時，則根據定積分的幾何定義，我們有其中 為小區間長度分點為 ，它們所對應的縱坐標函數分別為 。如圖5-2，取其中任一個小區間 進行分析。該區間為一曲邊梯形，y用以下方法求其面積。1I()()bniiiaf x dxfxibaxxn ix0121,kkn

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

5、分法1.算法分析 設在區間(a,b)上有可積函數f(x),求積分值將區間(a,b)分成n個相等的小區間，每個小區間之長為：h=(b-a)/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個區間：Ti=(f(xi)+f(xi+1)*h/2 第i-1個區間：

6、Ti-1=(f(xi-1)+f(xi)*h/2 求和得到：定步長梯形積分法的截斷誤差為:1101T()()/2(),1,2,3,-1nniiiihf af bf xxaxxhin(5-3)(1)02()()()d()d(1)!()()12nnbbnkaakfR ff xpxxxxxnbah f5.2 梯形積分法 該誤差按h的2次方的速度下降,即定步長梯形積分法的誤差階為2，也就是說，定步長梯形積分法具有一階代數精度因為代數精度等于誤差階減1，所以2-1=1）。2.定步長梯形積分法的程序框圖與通用程序設計 定步長梯形積分法的通用程序框圖如圖5-3所示，共分為三個程序框圖，即在圖中1為調用子程序的

7、主控程序框圖，（2為函數子程序框圖，（3為積分通用子程序框圖。5.2 梯形積分法（1主控程序框圖（2函數子程序框圖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 subroutine（3梯形積分子程序框圖5.2 梯形積分法10C 主控程序20 PROGRAM main3

8、0 read(5,*)a,b,n40 call SUBRPUTINE50&txjf(a,b,n,T)60 write(6,*)T70 END80C 函數子程序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(x)200 END DO210 T=T*h220 END SUBRPUTINE txjf 11()niif x

9、5.2 梯形積分法二、變步長梯形積分法 前面介紹了定步長梯形積分法。采用定步長梯形積分法求定積分時，步長的選擇應該是適當的。如果采用定步長梯形積分法，只有改變N值，比較幾組計算結果才能判是否達到精度要求，這樣程序的執行就受到人工干擾，占用了機時，降低了計算機的效率。為此，需要尋求新的改善計算精度途徑。截斷誤差分析表明，只要步長h充分小，必可滿足精度要求。因而，通常采取逐步縮小步長h的辦法。自動變步長梯形積分法是計算機根據精度要求進行自動判斷，逐步逐步縮小步長h。1.算法分析 其基本思想為：逐步變更步長，用二分法使步長逐次變小，直到滿足精度。5.2 梯形積分法N=1:N=1:h1=b-a h1=

10、b-a T1=f(a)+f(b)T1=f(a)+f(b)*h1/2 h1/2 N=2:N=2:h2=h1/2 h2=h1/2 T2=f(a)+f(b)T2=f(a)+f(b)*h2/2+f(x1)h2/2+f(x1)*h2h2 =T1/2+f(x1)=T1/2+f(x1)*h2 h2 N=4:N=4:h4=h2/2h4=h2/2T4=f(a)+f(b)T4=f(a)+f(b)*h4/2+f(x1)+f(x2)+f(x3)h4/2+f(x1)+f(x2)+f(x3)*h4h4 =T2/2+f(x1)+f(x3)=T2/2+f(x1)+f(x3)*h4 h4 abab1xab1x2x3x5.2 梯

11、形積分法N=2n:N=2n:注：注：xi=a+(2xi=a+(2*i-1)i-1)*h2nh2n：上一次積分區間的中點。：上一次積分區間的中點。在計算在計算 的過程中，每當算出一新的近似值的過程中，每當算出一新的近似值時，便檢查判斷下列條件是否成立：時，便檢查判斷下列條件是否成立：如果條件成立，那么如果條件成立，那么 就是符合精度要求的積分近似值，否就是符合精度要求的積分近似值，否則縮小步長則縮小步長h2nh2n再行計算，直到滿足條件為止。再行計算，直到滿足條件為止。22221/2 (2-1)TT/2()nninnnnniihhxaihhf x12482,nnT T T TT T2nnTTE（

12、E為容許誤差）2nT(5-4)5.2 梯形積分法變步長梯形積分法的截斷誤差為:2.變步長梯形積分法的程序框圖與通用程序設計 變步長梯形積分法的通用程序框圖如圖5-4所示，也共分為三個程序框圖，即在圖中1為調用子程序的主控程序框圖，（2為函數子程序框圖，（3為積分通用子程序框圖。（1和2與上一節雷同，因此這里只畫出梯形變步長積分法通用子程序框圖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

13、*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.應用實例分析例5.1 已知圓周率的近似值可以由下列定積分求出試求出圓周率的似近值。解：因為所以求解此問題的函數程序程序為80C 函數子程序90 FUNCTION f(x)100 f=4.0/(1+x*x)110 END FUNCTION f 12041dxx24()1f xx5.2 梯形積分法n n10102020404

14、0100100200200300300T T33139926331399263.1411763.1411763.1414893.1414893.1415763.1415763.1415893.1415893 3141591141591表5.1 梯形積分法輸入不同N地的計算機輸出結果 由上表可以看出，定步長梯形積分法的收斂速度是很慢的，由上表可以看出，定步長梯形積分法的收斂速度是很慢的，對于本例，對于本例，N=300N=300時，即把積分區間時，即把積分區間00，11劃分為劃分為300300個小區個小區間時才使計算結果精確到間時才使計算結果精確到 0.000010.00001。5.3 Simps

15、on積分法 梯形積分法是用直線段一次二項式去逼近曲線f(x),而Simpson 是用拋物線段二次三項式去逼近曲線f(x).一、定步長Simpson積分法1.算法分析 設將a,b分為n個偶數相等的小區間，那么 求積節點坐標 對應求積點函數值()baSf x dx(-)/hb an012,nxa xxxb011(),(),()nyf ayf xyf b5.3 Simpson積分法被積曲線用 y=f(x)表達，而逼近曲線則采用 y=Ax2+Bx+C 0 xa1x2xix1ix2ixnxb()yf x2yAxBxC5.3 Simpson積分法第一段：x0,x1,x2 有同理，第二段：x2,x3,x4

16、有第末段：xn-2,xn-1,xn 有200021112222yAxBxCyAxBxCyAxBxC求出系數求出系數A,B,A,B,C C 2021012 ()d(4)/3xxSAxBxCxyyyh2234(4)/3Syyyh/221(4)/3nnnnSyyyh5.3 Simpson積分法共有n/2小塊面積，對它們求和即可得Simpson法求積公式或這里01234-1013-124-2(42424)/3()4()2()/3nnnnnSyyyyyyyhyyyyyyyyh(5-5)/2/2 121211/2/221211()()4()2()/3()-()4()2()/3nnijijnnijijSf

17、af bf xf xhf af bf xf xh(5-6)212(21),2iixaihxaih5.3 Simpson積分法同理我們可以推出辛普森積分法的截斷誤差為：由上式可知辛普森積分法的誤差階和代數精度分別為4和3，在積分小區間總數與梯形積分法相同，辛普森積分法的精度較高。2.通用程序框圖與通用程序設計 這里只畫出Simpson定步長積分法通用子程序框圖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)

18、+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 SUBRPUTINE 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

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

20、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)*H2/3 2|-|nnSSab1xab1x2x3x5.3 Simpson積分法N=2n:N=2n:關于變步長關于變步長SimpsonSimpson積分法的遞推計算公式以及程序框圖和通用積分法的遞推計算公式以及程序框圖和通用FORTRANFORTRAN程序將在后面幾種積分比較部分提示。這里給出一個程序將在后面幾種積分比較部分提示。這里給出一個通用的通用的FORTRANFORTRAN程序清單：程序清單：12212211()()4(

21、)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)+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 60 n=260 n=2*n;h=h/2.0n;h=h/2.070 s2=s2+s170 s2=s2+s180 s

22、1=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 if(abs(s2n-s)E)then170 s=s2n180 goto 60190 endif200 END SUBROUTINE Simpb 5.4 Newton-Cotes積分法 在梯形積分法中，是用一次二項式直線去逼近曲線f(x);在Simpson積分法中，是用二次三項式拋物線去逼近曲線f(x)。而在N-

23、C積分法中，用的是m 次多項式去逼近曲線f(x)。設在被積區間a,b選取m+1個插值點，用m 次多項式pm去逼近曲線f(x),那么：利用拉格郎日插值多項式 ()()bbmaaf x dxpx dx()0()()(),0,1,2,bmmkkkaf x dxbaCf xkm(5-9)5.4 Newton-Cotes積分法其中:Ck(m)為Cotes系數，可以具體求出起其值：m=1:C0(1)=1/2;C1(1)=1/2m=2:C0(2)=1/6;C1(2)=4/6;C2(2)=1/6見下表 注：在實際應用中，一般m9，否則會產生計算不穩定。上面所談到的問題，是在整個區間a,b上用Pm去逼近f(x)

24、。下面看看用N-C公式即(5-1)式用在小區間時情況，將a,b分成n個小區間。有：()0(1)(1)(1)(1)()!()!mm kmkCt ttktktm dtm kmk5.4 Newton-Cotes積分法 K m 0 0 1 1 2 2 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/14

25、425/9625/9619/28819/288 641/84041/8409/359/359/2809/28034/10534/1059/2809/2809/359/3541/84041/840 7751/17751/172802803577/173577/172802801323/171323/172802802989/1722989/17280802989/1722989/17280801323/1721323/17280803577/173577/17280280751/172751/1728080 8989/28989/283503505838/285838/28350350-928/

26、283928/283505010496/2810496/28350350-4540/2834540/283505010496/2810496/28350350-928/283928/28350505838/285838/28350350989/283989/28350505.4 Newton-Cotes積分法 x0=a,x1,x2,xn-1,xn=b,h=(b-a)/n當項次m=1時：x0,x1 x1,x2xn-1,xn 1(1)(1)10001110()()()()()()/2xxf x dxxxCf xCf xf af xh2121()()()/2xxf x dxf xf xh11()()

27、()/2nnxnxf x dxf xf b h5.4 Newton-Cotes積分法求和得到*復化梯形求積公式同理，當m=2時，在每一個小區間xk,xk+1上增加一個分點中點），xk+1/2=（xk+xk+1）/2=xk+h/2，此時區間a,b上共有2n+1個等距求積節點，即：x0,x0+1/2,x1,x1+1/2,x2,xn+1/2,xn 兩節點之間的距離新的步長h=h/2。因而，m=2時的求積公式為1111()()()2()/2()()2()()/2bniianiif x dxf af bf xhf af bf xban(5-10)5.4 Newton-Cotes積分法*復化Simpson

28、求積公式求和當項次m=4時:在xk,xk+1上增加三個分點xk+1/4,xk+2/4,xk+3/4 使xk+1/4=xk+h/4，xk+2/4=xk+2h/4，xk+3/4=xk+3h/4，此時區間a,b上共有4n+1個等距求積節點，兩節點之間的距離新的步長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

29、bf xf xhf af bf xf xban(5-11)5.5 Romberg積分法*復化Newton-Cotes求積公式一、Romberg積分算法分析 Romberg積分法是用計算精度較低的二分前后兩次積分近似值，通過誤差補償法來獲得收斂速度較快的一種積分方法。1、梯形積分法的誤差估計見式3-3）(線性組合)T*-T2n(T2n-Tn)/3 其中T*為精確值 T*T2n+(T2n-Tn)/3 4*T2n/3-Tn/3 =Sn(5-13)式說明 T2n的誤差約等于T2n-Tn）/3,將這個誤差作為一種補償加到T2n上，可期望得到更精確的Simpson的積分公式 它是梯形積分法二分前后兩次積分

30、近似值進行線性組合。(5-13)(5-14)5.5 Romberg積分法2 2、SimpsonSimpson積分法的誤差估計積分法的誤差估計(線性組合線性組合)S S*-S2n(S2n-Sn)/15-S2n(S2n-Sn)/15 其中其中S S*為精確值為精確值 S S*S2n+(S2n-Sn)/15 16S2n+(S2n-Sn)/15 16*S2n/15-Sn/15 S2n/15-Sn/15 =Cn =CnSimpsonSimpson積分法二分前后兩次積分近似值進行線性組合，得到積分法二分前后兩次積分近似值進行線性組合，得到Newton-CotesNewton-Cotes的積分公式精度更高）

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

32、意，此時不能再對的積分公式精度更高）。注意，此時不能再對RombergRomberg的積分進行線性加速。因為此時插值多項式已經為八的積分進行線性加速。因為此時插值多項式已經為八次，它具有九次插值多項式的精度。依據插值多項式的理論，次，它具有九次插值多項式的精度。依據插值多項式的理論，當插值多項式的次數超過九次時，將出現不穩定情況。因而，當插值多項式的次數超過九次時，將出現不穩定情況。因而，到此為止。到此為止。RombergRomberg積分法不僅精度高，而且它從最簡單的梯形積分開始，積分法不僅精度高，而且它從最簡單的梯形積分開始，逐步加工成收斂速度快的逐步加工成收斂速度快的CotesCotes

33、積分法。積分法。二、二、RombergRomberg積分法通用程序框圖與通用程序設計積分法通用程序框圖與通用程序設計1.1.算法過程算法過程5.5 Romberg積分法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)/63c1=c2k=3?r

34、1=r2abs(r2-r1)E?YNNNNYNYYsubroutine圖5-7 Romberg法子程序框圖5.5 Romberg積分法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 190 230 goto 190 240 endif2

35、40 endif250 end subroutine romberg250 end subroutine romberg3.3.程序語言程序語言5.5 Romberg積分法 讀懂本程序的關鍵在于看懂控制變量k(步長折減次數）。當k=1時，程序剛開始，程序只求出t2和s2(實際上s1)；當k=2時，步長折二分，t1=t2,s1=s2,且程序又求出新的t2(實際上為t4)和s2，同時求出c2，并且令c1=c2；當k=3時，步長再折二分，t1=t2,s1=s2,同時又求出新的t2(實際上為t8)和s2(實際上為s4)。求出新的c2,且執行c1=c2,求出r2,且執行r1=r2,c1=c2；當k=4時

36、，步長再折二分，t1=t2,s1=s2,同時又求出新的t2(實際上為t16)和s2(實際上為s8)以及c2(實際上為c4)，新的r2。以后，只要當k3，程序每次都執行判斷abs(r2-r1)是否滿足精度要求，若滿足，則結束，否則執行r2=r1和c1=c2，且重復k=4只是k值不同的過程，直至滿足精度要求。5.6 Gauss積分法 前面介紹的幾種積分法都是在等分區間情況下進行討論的，而Gauss積分法則是適用于不等分區間的一種非常有效的數值積分方法。一、高斯點及高斯公式簡介 首先介紹一下代數精度的概念。一個求積公式，如果對于任何一個不超過m次多項式fm(x)均能使其余項積分為零，即Rf=0，但對

37、于m+1次多項式不能能使其余項積分為零，即Rf0,則稱此求積公式具有m次代數精度。內插求積公式可以表示為0()d()nbiiaif xxA f x(5-19)5.6 Gauss積分法 在等距節點的情況下，我們可以驗證它具有n次代數精度。假定若節點數目一定，我們能否通過適當地選擇節點的位置，使積分公式具有更高的代數精度呢？Gauss首先指出，若我們選擇的接點x0,x1,x2,xn所構成的多項式 與任意次數不超過n次的多項式p(x)都正交，即 則因 為次數不超過2n+1次的多項式，故求積公式(5-19)式具有2n+1次代數精度?？紤]積分012()()()()()nxxxxxxxxx()()d0ba

38、x p xx()()x p x()dbaf xx5.6 Gauss積分法 的上下限a,b是變動的，通過積分變換可將上下限變為固定的值 1和1。事實上，令 ，可以將積分區間a,b變換為-1,1，即因而，不失一般性，可考察區間-1,1上的Gauss積分公式()/2()/2xba tba11()d()d222babababaf xxftt110()d()niiif xxA f x(5-20)5.6 Gauss積分法常用的Gauss積分公式如下：（1Gauss-Legendre積分公式（4Gauss-Hermite積分公式（2Gauss-Chebyshev積分公式（3Gauss-Laguerre積分公

39、式110()d()niiif xxA f x11021()d(cos)122niiif xxA fnn00()d()nxiiief xxA f x20()d()nxiiief xxA f x5.6 Gauss積分法nxiAinxiAi00230.86113630.33998100.34785480.652145211208/95/9400.90617980.53846930.56888890.23692690.47862873/53/3表Gauss-Legendre求積公式的xi與Ai5.6 Gauss積分法nxiAinxiAi01130.32254771.74576114.53662039.

40、35907090.60315410.35741870.03888790.000539310.58578643.41421360.85355340.146446620.41577462.29428046.28994510.71109300.27851770.010389340.26356031.41340313.59642587.085810012.64080080.52175560.39866680.07594240.00361180.0000234表2 Gauss-Laguerre求積公式的xi與Ai5.6 Gauss積分法nxiAinxiAi001.772453930.52464761.6

41、5068010.80491410.081312810.70710680.8862269201.22474491.18163590.2954110400.95857262.02018290.94530870.39361630.0201932表3 Gauss-Hermite求積公式的xi與Ai5.6 Gauss積分法例用階數n=2的Gauss求積公式計算以下積分的近似值解：由于n=2，在表中查出Gauss求積點坐標和權系數1204Id1xx01013/3,3/3,1xxAA 001122I()()22222144112131131112322323.14754086bababababaA fxA

42、fx 5.6 Gauss積分法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.6 Gauss積分法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

43、.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 data 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

44、 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&f(xm+dx)150 end do150 end do160 s=xr160 s=xr*s s170 end subroutine gausint170 end subroutine gausint5.7 數值微分法 在現代工程技術領域和科研過程中，數值微分法雖然不像數值積分法那樣用途廣泛，但還是經常碰到數值微分問題。嚴格地說，數值微分法應叫做數值導

45、數法微分是差，而導數才是差商的極限形式），但由于歷史的原因，我們這里仍延用數值微分這一概念。一、差商法求數值微分1、算法分析 考慮一個在點附近能夠解析的函數f(x)，在 點它可以展開為Taylor級數，并去前兩項，有 ()ixh2()()()()iiif xhf xhfxO h(5-21)5.7 數值微分法其中 為 的同階無窮小。由5.21式求得函數的一階導數表達示為我們稱5.22式為f(x)在 點的誤差階為h的向前差分式，即同樣使用泰勒級數將f(x)在 點展開并整理可得我們稱5.24式為f(x)在 點的誤差階為h的向后差分式2()O h2h()()()()iiif xhf xfxO hh(5

46、-22)()()()iiif xhf xfxh(5-23)()ixh()()()iiif xf xhfxh(5-24)ixix5.7 數值微分法取5.23式和5.24式的算術平均值，有（5.25式稱為f(x)在 點一階中心差商近似公式。2、程序框圖與程序設計()()()2iiif xhf xhfxh(5-25)ix5.7 數值微分法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

47、圖5-9 差商數值微分子程序框圖5.7 數值微分法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 數值微分法二、插值法求數值微分1、算法分析 用插值公式求數值微分有兩點公式、三點公式、四點公式和五點公式等等，下面介紹常用的三點公式。已知f(x)在等距節點 的函數值 ，且函數f(x)的二階導數存在，作三點二次插值函數F(x)并對F

48、(x)求導得三點數值導數公式如下：ix(),1,2,if xik kk2123()4()()()()2kkkkf xf xf xfxO hh221()()()()2kkkf xf xfxO hh(5-26)(5-27)5.7 數值微分法即用這種方法求數值微分時，節點數必須是3的倍數。要特別強調的是，當插值多項式pn(x)收斂于f(x)時，不一定能保正也收斂于 ，而且當節點間距h縮小時，雖然截斷誤差逐步減小，但舍入誤差卻可能會增大，因此減小步長h不一定能提高計算精度。2、程序框圖與程序設計 212()4()3()()()2kkkkf xf xf xfxO hh(5-28)0,3,6,9,3(0,

49、1,2,3,)ki in()npx()fx5.7 數值微分法subroutinex0=x;x1=x+h;x2=x+2hdo x=a,b,3hcomp f(x0),f(x1),f(x2)END DO xEND subroutineidf1=(4f(x1)-3f(x0)-f(x2)/2hidf2=(f(x2)-f(x0)/2h圖5-10 插值數值微分子程序框圖idf3=(f(x0)-4f(x1)+3f(x2)/2h5.7 數值微分法10 subroutine intpder(a,b,h)20 do x=a,b,3*h30 x0=x;x1=x+h;x2=x+2h40 idf1=(4*f(x1)-3*

50、f(x0)-f(x2)/h/250 idf2=(f(x2)-f(x)/h/260 idf3=(f(x)+3*f(x2)-4*f(x1)/h/270 end do80 end subroutine intpder5.7 數值微分法三、實用的通用求數值微分 前兩節分別介紹了連續函數和離散函數表格函數的數值微分計算方法，已經能夠解決現代科學或工程技術領域或科研工作中所遇見的絕大部分微分問題。但是有些數值計算過程要求的精度較高，用前面介紹的方法求出的數值導數滿足不了要求，所以本章節再介紹一種精度較高的連續函數通用數值微分方法。1、算法分析已知連續函數f=f(x)，求任意點 處的導數 ，可用下述方法：此

51、時ix()fx0.5iixx1,imm5.7 數值微分法 此時 此時則有這里m是210之間的正整數，一般取m=10。2、程序框圖與通用程序設計 通用的數值微分法的程序框圖見圖5-11，下面是根據框圖編寫的通用程序：1()()iif xf xdxx1im()()iif xf xdxxim1()2ifxdd5.7 數值微分法subroutined=0do x1=a,b,hcomp f(x),f(x1)END DO x1END subroutineX=x1+0.5*iD=(f(x)-f(x1)/(x-x1)圖5-11 實用數值微分子程序框圖d1=ddo i=9,10END DO i5.7 數值微分法

52、10 subroutine popder(a,b,h)20 do x1=a,b,h30 d=040 do i=9,1050 d1=d60 x=x1+0.5*i70 d=(f(x)-f(x1)/(x-x1)80 end do90 end do100 end subroutine popder練習題5-1 利用9節點定步長復化梯形與Simpson求積公式計算下列積分5-2 已知函數y=f(x)的數值如下，利用差商或插值數值微分公式計算f(1)。x 0.6 0.8 0.9 1.0 1.1 1.2 1.4 y 0.736 0.8365 0.9095 1.0000 1.1105 1.2446 1.6017/2100sinln(1)(1)d,(2)dxxxxxx上機實習題5-3 用變步長Simpson求積公式與Romberg求積公式計算5-1題中積分。