FORTRAN數(shù)值方法及其在物理學(xué)中應(yīng)用

1、FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用1 第三章 物理學(xué)中定積分的數(shù)值計算方法物理學(xué)中定積分的數(shù)值計算方法 3.1 定積分基本數(shù)值算法及其應(yīng)用定積分基本數(shù)值算法及其應(yīng)用 3.2 龍貝格法及其應(yīng)用龍貝格法及其應(yīng)用 3.3 高斯求積法高斯求積法 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用2 3.1 定積分基本數(shù)值算法及其應(yīng)用定積分基本數(shù)值算法及其應(yīng)用 一、矩形法、梯形法和拋物線法一、矩形法、梯形法和拋物線法(辛普森法辛普森法) 1. 矩形法矩形法 矩形法計算定積分示意圖矩形法計算定積分示意圖 b a N i i xxfdxxfI 1 0 )(

2、)( N ab x )( 考慮定積分考慮定積分 ： b a dxxfI)( FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用3 2. 梯形法梯形法 N ab x )( xxfxfS iii )()( 2 1 1 曲線下所有梯形面積和為：曲線下所有梯形面積和為： xxfxfxxfxfxxfxfSS NNi )()( 2 1 )()( 2 1 )()( 2 1 12110 , 2 1 0 N CC1 121 N CCC其中：其中： b a N i ii xxfCdxxfI 0 )()( FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用4 3. 拋物線法（辛普

3、森法）拋物線法（辛普森法） v把區(qū)間把區(qū)間 分成分成 （偶數(shù)）（偶數(shù)） ,baN N ab x )( v整個曲線整個曲線 用用 個以拋物個以拋物 )(xf2N v設(shè)設(shè) , ,對于第一個曲邊四邊形有對于第一個曲邊四邊形有 CBxAxy 2 CBxAxxf CBxAxxf CBxAxxf 2 2 22 1 2 11 0 2 00 )( )( )( 個均等的小區(qū)間個均等的小區(qū)間 線為邊界的四邊形來替代線為邊界的四邊形來替代 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用5 該四邊形的面積為該四邊形的面積為: 2 0 2 0 )()( 2 0 x x x x dxCBxAxdxx

4、fS 顯然有顯然有 ， ，xxx2 02 102 2xxx )()(4)( 3 )( 2100 2 0 xfxfxf x dxxfS x x )()(4)( 3 )( 4321 4 2 xfxfxf x dxxfS x x )()(4)( 3 )( 6542 6 4 xfxfxf x dxxfS x x )()(4)( 3 )( 122 2 NNN x x N xfxfxf x dxxfS N N 因此因此 4)(2 )()()( 6 1 20 2 202 2 20 2 002 CxxB xxACBxAxCBxAxxx FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用6 積

5、分值可近似可以表示為：積分值可近似可以表示為： )()(4)(2 )(2)(4)(2)(4)( 3 12 43210 NNN xfxfxf xfxfxfxfxf x N i ii xxfC 0 )( ,4,2 3 2 ,3 , 1 3 4 ,0 3 1 i i Ni C i 需要說明的是被積函數(shù)有時并未寫需要說明的是被積函數(shù)有時并未寫 成成 的形式，給出的被積函數(shù)的形式，給出的被積函數(shù) 是一組離散的數(shù)據(jù)點是一組離散的數(shù)據(jù)點 ，此時三，此時三 種計算方法仍舊適用，只需將公式種計算方法仍舊適用，只需將公式 中中 的換成的換成 即可。即可。 )(xfy ),( ii yx )( i xf i y b

6、 a x x N dxxfdxxfI 0 )()( FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用7 例例1: 將區(qū)間將區(qū)間 兩等分，用矩形法、梯形法、兩等分，用矩形法、梯形法、 ) 2 , 0( 矩形法：矩形法： xxfxxfxxfdxxfG b a i i )()()()( 10 1 0 1 梯形法：梯形法： 兩等分：兩等分： 4 x xxfxxfxxfxxfCdxxfG b a i ii )( 2 1 )()( 2 1 )()( 210 2 0 2 xx) 4 cos(0cos 4 ) 2 cos 2 1 4 cos10cos 2 1 ( 精確值：精確值：1G 2

7、0 cos xdxG拋物線法計算積分拋物線法計算積分 3408. 1 4 70711. 1 42 2 4 1 94806. 0 4 20711. 1)0 2 1 2 2 11 2 1 ( 4 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用8 拋物線法：拋物線法： xxfxxfxxfxxfCdxxfG b a i ii )( 3 1 )( 3 4 )( 3 1 )()( 210 2 0 3 0.00.20.40.60.81.01.21.41.6 0.0 0.2 0.4 0.6 0.8 1.0 2 1/2/2=0.707 /4 y=cos(x) X Y 若把區(qū)間若把區(qū)間10等

8、分，有等分，有 0764. 1 1 G 9974. 0 2 G 000003. 1 3 G 4 ) 2 cos 3 1 4 cos 3 4 0cos 3 1 ( 00228. 1 4 27615. 1)0 3 1 2 2 3 4 1 3 1 ( 4 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用9 open(1,file=int.dat) write(*,*)input a,b,N=? read(*,*)a,b,N C method 1：矩形法：矩形法 y1=0.0 do 10 j=0,N-1 x1=a x1=x1+float(j)*(b-a)/float(N) 10 y

9、1=y1+f(x1)*(b-a)/float(N) write(1,*)N,y1 write(*,*)N,y1 計算程序計算程序 b a N i i xxfdxxfI 1 0 )()( FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用10 C method 2：梯形法：梯形法 y2=0.0 do 20 j=0,N x2=a x2=x2+float(j)*(b-a)/float(N) If(j.eq.0.or.j.eq.N) then y2=y2+0.5*f(x2)*(b-a)/float(N) else y2=y2+f(x2)*(b-a)/float(N) end if 2

10、0 continue write(1,*)N,y2 write(*,*)N,y2 計算程序計算程序 b a N i ii xxfCdxxfI 0 )()( , 2 1 0 N CC 1 121 N CCC FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用11 計算程序計算程序 C method 3：拋物線法：拋物線法 N: 偶數(shù)偶數(shù) y3=0.0 do 30 j=0,N x3=a x3=x3+float(j)*(b-a)/float(N) If(j.eq.0.or.j.eq.N) then y3=y3+1./3.*f(x3)*(b-a)/float(N) else k=j-

11、2*int(j/2) if(k.eq.0) then y3=y3+2./3.*f(x3)*(b-a)/float(N) else y3=y3+4./3.*f(x3)*(b-a)/float(N) end if end if 30 continue write(1,*)N,y3 write(*,*)N,y3 end function f(x) f=cos(x) end b a N i ii xxfCdxxfI 0 )()( 1 0, 3 4 1,3, 3 2 2,4, 3 i iN Ci i FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用12 nEX3-1：用矩形法、梯形法

12、、拋物線法編程：用矩形法、梯形法、拋物線法編程 計算定積分計算定積分( ) 0 2 ) 1 1 cos(dx x 1000N 作業(yè)作業(yè) FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用13 二、二、電磁學(xué)中數(shù)值積分的應(yīng)用電磁學(xué)中數(shù)值積分的應(yīng)用 點電荷的電勢：點電荷的電勢： r Q k r Q u 0 4 1 （庫侖常數(shù)米（庫侖常數(shù)米.伏伏/庫；真空中介電常數(shù)為庫；真空中介電常數(shù)為 法法/米）米） 9 109k 12 0 1085. 8 1 1電勢的計算電勢的計算 Q電量為電量為 的點電荷，的點電荷， 處產(chǎn)生的電勢為：處產(chǎn)生的電勢為：r在距離在距離 ),(yxPr Q * *

13、 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用14 例例2：設(shè)有一長直導(dǎo)線均勻帶電，線電荷密度為：設(shè)有一長直導(dǎo)線均勻帶電，線電荷密度為 ， 解：在導(dǎo)線上取一小段解：在導(dǎo)線上取一小段 ，視為點電荷，其電量為，視為點電荷，其電量為 ，dxdx 212 0 2 000 )(4 1 4 1 yxx dx r dx du l l yxx dx duu 212 0 2 00 )(4 l2P長度為長度為 ，求空間任一點，求空間任一點 的電勢。的電勢。 P它在它在 點產(chǎn)生的電勢為：點產(chǎn)生的電勢為： FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用15 解析法可得結(jié)果

14、：解析法可得結(jié)果： )()( )()( ln 4 0 212 0 2 0 0 212 0 2 0 0 lxylx lxylx u 數(shù)值法：數(shù)值法： 212 0 2 00 )( 1 4 )( yxx xf 取取 llba, 9 00 , ,5 10 ,0.125,0, =100 x yN 075. 0l 將將 代入并用三種數(shù)值方法即可得到數(shù)值解。代入并用三種數(shù)值方法即可得到數(shù)值解。)(xf dxxfu b a )( v解析解：解析解： 62.326450 u v數(shù)值解：數(shù)值解： 326450.62 329600.62 823810.61 u u u 拋物線法： 梯形法： 矩形法： FORTRAN

15、數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用16 x 例例3：求帶電圓環(huán)（半徑為：求帶電圓環(huán)（半徑為 ）在軸線上點）在軸線上點 的電勢的電勢 （線電荷密度（線電荷密度 ） r l ds du 0 4 1 2/122 00 2 0 0 )(4 2 4 1 4rx Q l r ds l duu r X l x Y r -1.0-0.50.00.51.0 0 u x r=0.08 r=0.16 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用17 例例4：帶電圓盤（半徑為：帶電圓盤（半徑為 ）在軸線上）在軸線上 點的電勢點的電勢 Rx 2/122 0 2/122 0 )

16、( 2 4 1 )(4 1 rx rdr rx dQ du 環(huán)帶電量環(huán)帶電量 ，在，在 處的電勢為：處的電勢為：xrdrdQ2 解析解：解析解： RR rx rd rx rdr 0 2/122 2 0 2/122 )( )( 2 1 )( |)(|)( 2/122 0 2/122 xRxrx R = = |)( 2 2/122 0 xRxu 數(shù)值解：數(shù)值解： R drrfu 0 )( 2/122 0 )(2 )( rx r rf FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用18 -0.10-0.08-0.06-0.04-0.020.000.020.040.060.080

17、.10 0.00E+000 u x R=0.08 R=0.16 圓盤軸線上電勢分布圓盤軸線上電勢分布 練習(xí)：練習(xí)：(1)(1)5 . 0,08. 0,06. 0Rx 用三種數(shù)值積分方法計算例用三種數(shù)值積分方法計算例4中的電勢。中的電勢。 (2) ( (2) ( 隨隨 變化的曲線變化的曲線) )08. 008. 0 xux FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用19 2電場強度的計算電場強度的計算 電量為電量為 的點電荷，在距離的點電荷，在距離 處產(chǎn)生的電場強度為處產(chǎn)生的電場強度為Qr 點電荷的電場強度：點電荷的電場強度： rE 2 0 4 1 r Q )( r r

18、 r FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用20 解：圓環(huán)上任一點電荷解：圓環(huán)上任一點電荷 在軸線上點在軸線上點 的電場強度的電場強度 dsxEd 例例5：求帶電圓環(huán)（半徑為：求帶電圓環(huán)（半徑為 ）在軸線上點）在軸線上點 的電場強度的電場強度 rx )2(rQ 2 0 4 1 l ds dE 3 0 4 1 cos l xds dEdE x )(cos l x r x rx Qx l xds E 2 0 2/322 0 3 0 )( 1 44 1 （線電荷密度（線電荷密度 ，如下圖所示），如下圖所示） 大小為大小為 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在

19、物理學(xué) 中應(yīng)用中應(yīng)用21 -1.0-0.50.00.51.0 0 Ex x r=0.08 r=0.16 圓環(huán)軸線上電場強度分布圓環(huán)軸線上電場強度分布 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用22 例例6：已知圓環(huán)帶電量為：已知圓環(huán)帶電量為q q ，桿的線密度為，桿的線密度為 ，長為，長為L L。 求：桿對圓環(huán)的作用力。求：桿對圓環(huán)的作用力。 xqdd 解：解： 圓環(huán)在圓環(huán)在 dq 處產(chǎn)生的電場處產(chǎn)生的電場 2/322 0 )(4 1 xR qx Ex x E R q L qd x o 在桿上位置在桿上位置x處取一個電量為處取一個電量為dq的微元的微元 FORTRAN

20、數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用23 x E R q L qd x o L xR xxq F 0 2322 0 )(4 d xEqEF xx ddd ) 11 ( 4 22 0LR R q dq所受電場力為：所受電場力為： 取取q=1， =1，R=1，L=1 其中其中 F/m1082187854. 8 12 0 由上式可得由上式可得F=2.6323932E+09 a=0，b=1，N=100 矩形法矩形法 F=2.6164170E+09 梯形法梯形法 F=2.6323049E+09 辛普森法辛普森法 F=2.6323930E+09 帶電圓環(huán)與桿間作用力帶電圓環(huán)與桿間作用力F

21、: FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用24 0246810 R=1 R=2 F a 當(dāng)桿沿當(dāng)桿沿x軸正向平移距離軸正向平移距離a，桿對環(huán)的作用力又是多少？，桿對環(huán)的作用力又是多少？ FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用25 取一半徑為取一半徑為 ，寬為，寬為 的圓環(huán)帶，其電量的圓環(huán)帶，其電量rdrrdrdQ2 2/322 0 2/322 0 )( 2 4 1 )( 1 4xr rxdr xr xdQ dE R xr rdrx dEE 0 2/322 0 )(2 2/322 0 )(2 )( xr rx rf 其中其中數(shù)值解：數(shù)值

22、解： b a drrfE)(Rba,已知已知 例例7: 求均勻帶電圓盤在軸線上的電場強度。求均勻帶電圓盤在軸線上的電場強度。 為電荷為電荷 E整個圓盤的電場整個圓盤的電場 解析解：解析解： )( 1 2 2/122 0 xR x E R面密度，為圓盤半徑。面密度，為圓盤半徑。 dE它在軸線上的電場強度為它在軸線上的電場強度為 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用26 -0.10-0.08-0.06-0.04-0.020.000.020.040.060.080.10 Ex x R=0.08 R=0.16 圓盤軸線上電場強度分布圓盤軸線上電場強度分布 FORTRAN

23、數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用27 練習(xí)：練習(xí)：(1)(1)5 . 0,08. 0,06. 0Rx (2) ( (2) ( 隨隨 變化的曲線變化的曲線) )08. 008. 0 x E x （用三種方法，注意（用三種方法，注意 的選?。┑倪x?。?N 用三種數(shù)值積分方法計算例用三種數(shù)值積分方法計算例7中的場強中的場強 。E FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用28 EX3-2：對圓盤軸線上一點：對圓盤軸線上一點E 2 1 2/322 0 )(2 R R rx rdrx E 采用解析法和數(shù)值法（三種）計算采用解析法和數(shù)值法（三種）計算 01.

24、 0 , 5 . 0 x其中 m10 ,01. 0/, 21 RR （1 1） （2 2） 觀察觀察 隨隨 的變化。的變化。 m1 ,01. 0/, 21 RR E 2 R 作業(yè)作業(yè) FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用29 例例8：求載流圓線圈軸線上一點：求載流圓線圈軸線上一點 P 的磁感應(yīng)強度的磁感應(yīng)強度 2 00 d 4 d r rlI B 2 0 d 4 d r lI B 畢薩定律：畢薩定律： 圖中電流元在圖中電流元在P點產(chǎn)生磁感應(yīng)強度點產(chǎn)生磁感應(yīng)強度dB： )( d 4 22 0 xR lI I lI d B d B d P x R o r x 3磁場

25、的計算磁場的計算 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用30 根據(jù)對稱性根據(jù)對稱性0 B I lI d B d B d P x R o r x P x B 2/122 )( cos xR R r R 2/322 2 0 )(2xR IR B cos d 4 cosd 2 0 r lI B x BBd 方向滿足右手螺旋法則方向滿足右手螺旋法則 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用31 載流圓線圈軸線上磁感應(yīng)強度分布載流圓線圈軸線上磁感應(yīng)強度分布 0246810 -5.00E-008 0.00E+000 5.00E-008 1.00E-

26、007 1.50E-007 2.00E-007 2.50E-007 3.00E-007 3.50E-007 4.00E-007 4.50E-007 5.00E-007 5.50E-007 6.00E-007 6.50E-007 7.00E-007 B x R=2, I=1 R=1, I=1 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用32 例例9：求繞軸旋轉(zhuǎn)的帶電圓盤軸線上的磁場：求繞軸旋轉(zhuǎn)的帶電圓盤軸線上的磁場 x q r O R B d P 解：解： 2 / Rq 帶電圓盤電帶電圓盤電荷面密度為：荷面密度為： rrqd2d所取環(huán)帶上帶電量為：所取環(huán)帶上帶電量為： t

27、 q I d d rr rr d 2 d2 等效電流為：等效電流為： FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用33 q r O R B d P 2/322 3 0 2/322 2 0 )(2 d )(2 d d xr rr xr Ir B 圓環(huán)圓環(huán)dq在在P點產(chǎn)生的磁感應(yīng)強度點產(chǎn)生的磁感應(yīng)強度dB： x Rx xR dBB R 2 2 2 22 22 0 0 對上式積分：對上式積分： 取取=1，=600，R=1，x=1 其中其中 27 0 AN104 由上式可得由上式可得B=4.536746E-05 a=0，b=1，N=100 矩形法矩形法 B=4.5071927E

28、-05 梯形法梯形法 B=4.5738361E-05 辛普森法辛普森法 B=4.5736691E-05 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用34 繞軸旋轉(zhuǎn)的帶電圓盤軸線上的磁場分布繞軸旋轉(zhuǎn)的帶電圓盤軸線上的磁場分布 0246810 -0.00005 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 0.00040 B x R FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用35 三、三、分子物理中數(shù)值積分的應(yīng)用分子物理中數(shù)值積分的應(yīng)用 理想氣體在平衡態(tài)下分子的速率分

29、布函數(shù)理想氣體在平衡態(tài)下分子的速率分布函數(shù): kTv ev kT vf 2/22/3 2 ) 2 (4)( 是分子質(zhì)量，是分子質(zhì)量， 是氣體溫度，是氣體溫度， 為波耳茲曼常數(shù)。為波耳茲曼常數(shù)。 TJ/K 1038. 1 23 k 麥克斯韋速率分布律麥克斯韋速率分布律 dvev kT dvvf N dN kTv2/22/3 2 ) 2 (4)( dvvv氣體中速率在氣體中速率在 到到 間的分子數(shù)的比率為間的分子數(shù)的比率為: FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用36 例例10：利用氣體分子麥克斯韋速率分布律，求：利用氣體分子麥克斯韋速率分布律，求0度，度，127度度

30、 下氮分子、氧分子運動的速率分布曲線，并求當(dāng)速下氮分子、氧分子運動的速率分布曲線，并求當(dāng)速 率在率在300500m/s范圍內(nèi)的分子所占的比例，討論溫范圍內(nèi)的分子所占的比例，討論溫 度度 及分子量及分子量 對速率分布曲線的影響。對速率分布曲線的影響。T 解：根據(jù)麥克斯韋速率分布律可知，任一速率間隔解：根據(jù)麥克斯韋速率分布律可知，任一速率間隔 到到 中中 1 v 2 v 2 1 2 2 1 2/22/3 ) 2 (4)( v v kTv v v dvev kT dvvf N N 對于氮氣，分子質(zhì)量分別為對于氮氣，分子質(zhì)量分別為 kg/ kg/ ，而氧氣分子，而氧氣分子028. 0 0 N 的分子數(shù)

31、所占的比率可用積分法求出，即的分子數(shù)所占的比率可用積分法求出，即 0 N 0 N032. 0質(zhì)量質(zhì)量 kg/ ，其中，其中 =6.022為阿伏伽德羅常數(shù)。為阿伏伽德羅常數(shù)。 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用37 020040060080010001200 0.0 5.0 x10 -4 1.0 x10 -3 1.5x10 -3 2.0 x10 -3 2.5x10 -3 v (m/s) f(v) 氧氣， T=273K 氮氣， T=273K 氧氣， T=400K 氮氣， T=400K 氣體分子運動麥克斯韋速率分布曲線氣體分子運動麥克斯韋速率分布曲線 當(dāng)溫度升高和分

32、子量減小時，速率小的分子數(shù)減小而速率大的分子當(dāng)溫度升高和分子量減小時，速率小的分子數(shù)減小而速率大的分子 隨著溫度升高和分子量減小，分布曲線將變得越來越平坦。隨著溫度升高和分子量減小，分布曲線將變得越來越平坦。 氮氣分子所占比例為氮氣分子所占比例為0.3957 (T=273K)，0.3085(T=400K)； 氧氣分子所占比例為氧氣分子所占比例為0.4189 (T=273K)，0.3410(T =400K)。 數(shù)增多，分布曲線的極大值隨著溫度升高和分子量減小而向右移動。數(shù)增多，分布曲線的極大值隨著溫度升高和分子量減小而向右移動。 利用辛普森法可得，當(dāng)速率在利用辛普森法可得，當(dāng)速率在300500m

33、/s范圍內(nèi)：范圍內(nèi)： FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用38 v EX3-5: 結(jié)合理想氣體平衡態(tài)下的麥克斯韋速率分布律，利結(jié)合理想氣體平衡態(tài)下的麥克斯韋速率分布律，利 用辛普森法計算用辛普森法計算0 0C C下氮氣、氧氣分子運動的平均速率下氮氣、氧氣分子運動的平均速率 和和 方均根速率方均根速率 ，并與解析解，并與解析解 ， 做比較，其中做比較，其中 是氣體的摩爾質(zhì)量，是氣體的摩爾質(zhì)量， 為摩爾氣體常數(shù)。為摩爾氣體常數(shù)。 2 v M RT v59. 1 M M RT v73. 1 2 31. 8R EX3-3 p.67 EX3-4 p.67 作業(yè)作業(yè) FOR

34、TRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用39 dimension y(0:47) coe=2*(4.561*1.0e-5)*9.81/0.0956 data y/0.10135,0.20064,0.27303,0.31095,0.33094,0.33991,0.34474, * 0.35577,0.31508,0.29578,0.27717,0.26131,0.24545,0.23097,0.21718, * 0.20339,0.19167,0.17995,0.16823,0.15789,0.14824,0.13927,0.13288, * 0.12548,0.11859

35、,0.11238,0.10687,0.10204,0.09215,0.09308,0.08894, * 0.08480,0.08067,0.07722,0.07377,0.07032,0.06757,0.06481,0.06205, * 0.05929,0.05654,0.05378,0.05102,0.04826,0.04550,0.04274,0.04067, * 0.03861/ C method 3：辛普森法：辛普森法 y3=0.0 do 30 j=0,47 If(j.eq.0.or.j.eq.47) then y3=y3+1./3.*y(j)*0.0127 else k=j-2*in

36、t(j/2) if(k.eq.0) then y3=y3+2./3.*y(j)*0.0127 else y3=y3+4./3.*y(j)*0.0127 end if end if 30 continue write(*,*) sqrt(coe*y3) end FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用40 小結(jié)：小結(jié)： b a N i i xxfdxxfI 1 0 )()(1矩形法：矩形法： b a N i ii xxfCdxxfI 0 )()(2梯形法：梯形法： , 2 1 0 N CC1 121 N CCC其中：其中： b a N i ii xxfCdxxfI 0

37、 )()(3拋物線法拋物線法(辛普森法辛普森法)： , 4 , 2 3 2 , 3 , 1 3 4 , 0 3 1 i i Ni Ci其中：其中： 問題：問題： 1.二重積分問題？二重積分問題？ 2.不等分間隔問題不等分間隔問題 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用41 3.2 龍貝格法及其應(yīng)用龍貝格法及其應(yīng)用 三種簡單定積分的計算方法缺點是要在求積之前給三種簡單定積分的計算方法缺點是要在求積之前給 在實際計算中，經(jīng)常采用變步長的計在實際計算中，經(jīng)常采用變步長的計 算過程，即在步長逐步分半的過程中，算過程，即在步長逐步分半的過程中， 反復(fù)利用求積公式進行計算，直

38、到求反復(fù)利用求積公式進行計算，直到求 得的積分值滿足精度要求為止得的積分值滿足精度要求為止! 出步長，如何選擇最合適的步長來統(tǒng)一精度和計算出步長，如何選擇最合適的步長來統(tǒng)一精度和計算 量的矛盾往往有困難！量的矛盾往往有困難！ FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用42 一、變步長的梯形法一、變步長的梯形法 共有共有 個分點個分點1n,khaxk, n ab h nk, 2 , 1 , 0 )()( 2 11 kk xfxf h T )()(2)( 4 1 2 12 k k k xfxfxf h T n T用用 表示簡單梯形法所求得的積分值。表示簡單梯形法所求得的積

39、分值。 n,ba在簡單的梯形法中，將求積區(qū)間在簡單的梯形法中，將求積區(qū)間 分成分成 等份，則等份，則 ),( 2 1 1 2 1 kk k xxx其中點為其中點為, 1kk xx 1、考察一個子段、考察一個子段 在該子段上二分前和二分后的兩個積分值分別為在該子段上二分前和二分后的兩個積分值分別為 1 T 2 T 和和 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用43 2、利用這個、利用這個關(guān)系式當(dāng)關(guān)系式當(dāng) 從從0到到 累加累加時，可導(dǎo)出下面時，可導(dǎo)出下面 k1n 1 0 2 12 )( 22 1 n k k nn xf h TT 3、注意到：、注意到： , ) 2 1

40、( )( 22 1 2 1 1 0 2 12 hkax xf h TT k n k k nn 1 0 2 ) 2 1 ( 22 1 n k nn hkaf h TT 二分步長以后的積分值可由二分步長以前的積分值來計算二分步長以后的積分值可由二分步長以前的積分值來計算! 遞推公式：遞推公式： )()(2)( 4 1 2 12 k k k xfxfxf h T )()( 2 11 kk xfxf h T )( 22 1 2 112 k xf h TT FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用44 二、變步長的辛普森求積法二、變步長的辛普森求積法 用梯形法二分前后的兩個積

41、分值用梯形法二分前后的兩個積分值 和和 作如下線性組合，作如下線性組合， n T n T2 不進行證明！不進行證明！ 說明說明變步長梯形法算法簡單，但精度差，收斂慢變步長梯形法算法簡單，但精度差，收斂慢! | 2nn SS 1| 2 n S | 2 2 n nn S SS 1| 2 n S 其中其中 為允許的誤差限。為允許的誤差限。 )( 3 1 3 1 3 4 222nnnnnn TTTTTS ?。喝。?結(jié)果為辛普森積分值結(jié)果為辛普森積分值 h h 2 nn2重復(fù)上述積分過程，將求積區(qū)間逐步折半重復(fù)上述積分過程，將求積區(qū)間逐步折半( , ), n S2 n S 直到相鄰兩次的積分值直到相鄰兩

42、次的積分值 與與滿足下列關(guān)系：滿足下列關(guān)系： FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用45 例例11：我國第一顆人造地球衛(wèi)星軌道是一個橢圓，橢圓周長：我國第一顆人造地球衛(wèi)星軌道是一個橢圓，橢圓周長 2 /2 2 0 41sin c sad a 其中其中 是橢圓的半長軸，是地球的中心與軌道中心（橢圓中心）是橢圓的半長軸，是地球的中心與軌道中心（橢圓中心） 的距離，地球半徑的距離，地球半徑 =6371km，近地點距離，近地點距離 =439km， 遠地點距離遠地點距離 =2384km, , ac Rh H)2( 2 1 hHRa)( 2 1 hHc 試?yán)米儾介L辛普森法求

43、解定積分計算人造地球衛(wèi)星的試?yán)米儾介L辛普森法求解定積分計算人造地球衛(wèi)星的 的計算公式為：的計算公式為： 軌道周長。軌道周長。 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用46 subroutine simp(a,b,ep,s2,f) h=b-a t1=h/2.*(f(a)+f(b) n=1 5 s=0. do k=0,n-1 s=s+f(a+(k+0.5)*h) end do t2=t1/2.+h/2.*s s2=t2+(t2-t1)/3. if(n/=1) goto 20 15 n=n+n h=h/2. t1=t2 s1=s2 goto 5 20 if(abs(s2

44、)1.) d=abs(s2-s1)/s2) if(d=ep) goto 15 return end program main external f call simp(0.,1.5708,1e-5,s2,f) write(*,(20 x,s=,f10.4)s2 end function f(x) f=4*7782.5*sqrt(1.-(972.5/7782.5)*2 *sin(x)*2) end 計算程序計算程序 程序運行結(jié)果：程序運行結(jié)果： S=48707.5519 km )( 3 1 3 1 3 4 222nnnnnn TTTTTS 1 0 2 12 )( 22 1 n k k nn xf

45、 h TT | 2nn SS1| 2 n S | 2 2 n nn S SS 1| 2 n S FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用47 三、龍貝格求積法三、龍貝格求積法 不進行證明！不進行證明！ )( 15 1 22nnnn SSSC 重復(fù)上述步驟，便可導(dǎo)出下列龍貝格公式重復(fù)上述步驟，便可導(dǎo)出下列龍貝格公式: )( 63 1 22nnnn CCCR | 2nn RR1| 2 n R | 2 2 n nn R RR 1| 2 n R 其中其中 為允許的誤差限。為允許的誤差限。 為了使求積結(jié)果的精度更高，可以利用龍貝格法進行計算。為了使求積結(jié)果的精度更高，可以利用

46、龍貝格法進行計算。 n S2 n S用變步長辛普森法二分前后的兩個積分值用變步長辛普森法二分前后的兩個積分值 和和 作如下的作如下的 n C線性組合，得到的是柯特斯的積分值線性組合，得到的是柯特斯的積分值 重復(fù)積分過程，直到相鄰兩次的積分值重復(fù)積分過程，直到相鄰兩次的積分值 與與 滿足下列關(guān)系：滿足下列關(guān)系： n R2 n R FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用48 例例12：利用龍貝格法計算單擺的振蕩周期利用龍貝格法計算單擺的振蕩周期 T /2 1 0 22 0 2 4 1 sin ()sin 2 ld T g 2 m/s8 . 9g其中其中 ， ， ，m5

47、l/20 0 5 10 計算程序：計算程序： program main external f call romb(0.,1.5708,1e-5,r2,f) write(*,(3x,B=,e12.6) r2 end function f(x) f=4*sqrt(5/9.8)/(1-sin(10/3.14159)*2*sin(x)*2)*(1/2) end FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用49 subroutine romb(a,b,ep,r2,f) h=b-a t1=h/2.*(f(a)+f(b) n=1 5 s=0.0 do k=0,n-1 s=s+f(a+

48、(k+0.5)*h) end do t2=t1/2.+h/2.*s s2=t2+(t2-t1)/3. if(n/=1)goto 20 15 n=n+n h=h/2 t1=t2 s1=s2 goto 5 20 c2=s2+(s2-s1)/15. if(n/=2)goto 40 30 c1=c2 goto 15 40 r2=c2+(c2-c1)/63. if(n/=4) goto 60 50 r1=r2 goto 30 60 if(abs(r2)-1.) 70,70,80 70 if(abs(r2-r1)=ep) goto 50 return 80 if(abs(r2-r1)/r2)=ep) go

49、to 50 return end 程序運行結(jié)果如下：程序運行結(jié)果如下： B=.448800E+01 )( 15 1 22nnnn SSSC )( 63 1 22nnnn CCCR FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用50 1 2 2 2 2 b y a x dxxbxaL 0 2222 cossin 1 100400 22 yx EX3-6: 若橢圓方程為若橢圓方程為 ，則橢圓周長的計算公式為，則橢圓周長的計算公式為 利用辛普森法和變步長辛普森法計算橢圓利用辛普森法和變步長辛普森法計算橢圓 的周長。的周長。 EX3-7: 應(yīng)用變步長辛普森法編程計算積分應(yīng)用變步長

50、辛普森法編程計算積分 此積分值稱為此積分值稱為Catalan常數(shù)，常數(shù)， 的真值為的真值為0.915965。要求誤差。要求誤差 不超過不超過10-5。 1 0 arctan dx x x G G EX3-8:應(yīng)用龍貝格法計算積分應(yīng)用龍貝格法計算積分: 要求誤差不超過要求誤差不超過10-5。 6 . 0 1 . 0 2 . 18 . 0 )1 () 1449. 1 ( )2(02792. 0 dx xxx x I 作業(yè)作業(yè) FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用51 3.3 高斯求積法高斯求積法 一、代數(shù)精度、代數(shù)精度 點，用不同的點，用不同的m次多項式次多項式 0

51、 2 2 1 1 )(axaxaxaxy m m m m m m 矩形法：矩形法：0m 常數(shù)近似常數(shù)近似 2m辛普森（拋物線）法：辛普森（拋物線）法：拋物線近似拋物線近似 1m梯形法：梯形法： 線性近似線性近似 2、 可否任意大？可否任意大？mm越大，方法精度越高？越大，方法精度越高？1、是否、是否 前面介紹的三種求積法，都是限定用等分點作為求積節(jié)前面介紹的三種求積法，都是限定用等分點作為求積節(jié) 去近似被積函數(shù)去近似被積函數(shù) )(xf。 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用52 12 nm可能達到的最高代數(shù)精度為：可能達到的最高代數(shù)精度為： 1對有對有個積分點的

52、求積公式：個積分點的求積公式：n b a n k kk xfAdxxfG 0 )()( 給定的條件下，可找到一組節(jié)點的給定的條件下，可找到一組節(jié)點的n2. 在積分點的數(shù)目在積分點的數(shù)目 和對應(yīng)的系數(shù)和對應(yīng)的系數(shù) k x k A，使以上求積公式達到 ，使以上求積公式達到 坐標(biāo)坐標(biāo) 12 nm的代數(shù)精確度。的代數(shù)精確度。 結(jié)結(jié) 論論 n xxx , , , 10 適當(dāng)選取節(jié)點適當(dāng)選取節(jié)點的位置，可以使求積的位置，可以使求積公式精度公式精度 盡可能高。當(dāng)節(jié)點位置不同時，積分結(jié)果可能是不一樣的。盡可能高。當(dāng)節(jié)點位置不同時，積分結(jié)果可能是不一樣的。 說明說明 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及

53、其在物理學(xué) 中應(yīng)用中應(yīng)用53 v如上圖所示，當(dāng)節(jié)點位置不同時，積分結(jié)果是不一樣的。如上圖所示，當(dāng)節(jié)點位置不同時，積分結(jié)果是不一樣的。 圖圖 (a) 一般而言，代數(shù)精度越高，數(shù)值求積公一般而言，代數(shù)精度越高，數(shù)值求積公 式越準(zhǔn)確。要使求積公式具有式越準(zhǔn)確。要使求積公式具有 次代數(shù)精度，次代數(shù)精度， 只要令它對于只要令它對于 m xxxxf, 1)( 2 都能準(zhǔn)確成立即可。都能準(zhǔn)確成立即可。 m 圖圖 (b) a b FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用54 例例13：求證梯形公式：求證梯形公式: b a bfaf ab dxxf)()( 2 )( 具有一階代數(shù)精度

54、。具有一階代數(shù)精度。 證明：首先驗證證明：首先驗證xxf, 1)(時，梯形公式準(zhǔn)確成立。時，梯形公式準(zhǔn)確成立。 b a bfaf abab abdx)()( 2 11 2 1 b a bfaf ab ba ab abxdx)()( 2 2 )( 2 1 22 再驗證，當(dāng)再驗證，當(dāng) 2 )(xxf時，梯形公式能否準(zhǔn)確成立時，梯形公式能否準(zhǔn)確成立 b a bfaf ab ba abab dxx)()( 2 23 22 33 2 另外可以證明辛普森計算公式具有另外可以證明辛普森計算公式具有3階代數(shù)精度。階代數(shù)精度。 的形式，的形式，而而xcc 10 xxf, 1)( 由于任意一個一次多項式可表示成由

55、于任意一個一次多項式可表示成 ，梯形公式準(zhǔn)確成立，所以梯形公式具有，梯形公式準(zhǔn)確成立，所以梯形公式具有 一階代數(shù)精度。一階代數(shù)精度。 對于對于 梯形公式梯形公式 不成立！不成立！ FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用55 二、高斯型代數(shù)求積公式二、高斯型代數(shù)求積公式 先考慮積分區(qū)間先考慮積分區(qū)間-1-1，11上的求積公式上的求積公式 1 1 0 )()( n k kk xfAdxxf n xxx , , , 10 1n如果節(jié)點如果節(jié)點 是是次多項式次多項式 )()()()( 210n xxxxxxxxx )(xn的根，并且的根，并且與任意一個次數(shù)不超過與任意一個

56、次數(shù)不超過的多項式的多項式)(xq 正交，即：正交，即： 1 1 0)()(dxxqx 則求積公式對一切次數(shù)不超過則求積公式對一切次數(shù)不超過12 n 的多項式都準(zhǔn)確成立，的多項式都準(zhǔn)確成立， 求積系數(shù)：求積系數(shù)： 1 1 )()( )( dx xxx x A kk k 定定 理理 k x該定理未給出節(jié)點該定理未給出節(jié)點 的具體取法。的具體取法。 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用56 由特殊函數(shù)理論知，勒讓德（由特殊函數(shù)理論知，勒讓德（Legendre）多項式）多項式 ) 1( !2 1 )( 2n n n n n x dx d n xp 在在-1，1上是正交

57、的，即上是正交的，即 1 1 1 0)()(dxxpxp nn )( 1 xpn的首項系數(shù)為：的首項系數(shù)為： 21 )!1(2 )!1(2 n n n ，故常取，故常取 ) 1( )!1(2 )!1( )( )!1(2 )!1(2 )( 12 1 1 1 21 n n n n n x dx d n n xp n n x k x 的選取方法的選取方法 k A及求積系數(shù)及求積系數(shù)節(jié)點節(jié)點 1 1 0 )()( n k kk xfAdxxf 的的)( 1 xpn1n此時此時個零點就是求積式個零點就是求積式: : n xxx , , , 10 的節(jié)點的節(jié)點 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及

58、其在物理學(xué) 中應(yīng)用中應(yīng)用57 求積系數(shù)：求積系數(shù)： 1 1 )()( )( dx xxx x A kk k 經(jīng)過計算得：經(jīng)過計算得： 2 1 2 )()1 ( 2 knk k xpx A -1,1 ),( )!22( )!1( 32 2 )( )22( 3 432 n n f n n n fR 可以估得高斯求積公式的截斷誤差為：可以估得高斯求積公式的截斷誤差為： k x 的選取方法的選取方法 k A及求積系數(shù)及求積系數(shù)節(jié)點節(jié)點 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用58 幾個低階高斯求積公式及其余項：幾個低階高斯求積公式及其余項： (1) 點）(1 0n ,) 1

59、( 2 1 )( 2 1 xx dx d xp 1)( 1 x p 2 0 A可得：可得： 截斷誤差：截斷誤差： )( 3 1 )(ffR 0 0 x節(jié)點節(jié)點，由節(jié)點系數(shù)公式 由節(jié)點系數(shù)公式 2 1 2 )()1 ( 2 knk k xpx A 1 1 0 )()( n k kk xfAdxxf 1 1 )0(2)(fdxxf ) 1( !2 1 )( 2n n n n n x dx d n xp FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用59 （2） 點）(2 1n ) 13( 2 1 ) 1( 8 1 )( 222 2 2 2 xx dx d xpxxp3)( 2

60、 1 10 AA 2點求積公點求積公式：式： 1 1 1 0 11 ( )()()() 33 kk k f x dxA f xff 截斷誤差：截斷誤差： )( 135 1 )( )4( ffR , 3 1 0 x 3 1 1 x 1 1 ) 5 15 ( 9 5 )0( 9 8 ) 5 15 ( 9 5 )(fffdxxf 3點點 高斯求積公式為（不做推導(dǎo)）高斯求積公式為（不做推導(dǎo)） )2( n 2 1 2 )()1 ( 2 knk k xpx A ) 1( !2 1 )( 2n n n n n x dx d n xp FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用60