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

1、FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用1 第三章 物理學(xué)中定積分的數(shù)值計(jì)算方法物理學(xué)中定積分的數(shù)值計(jì)算方法 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. 矩形法矩形法 矩形法計(jì)算定積分示意圖矩形法計(jì)算定積分示意圖 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整個(gè)曲線整個(gè)曲線 用用 個(gè)以拋物個(gè)以拋物 )(xf2N v設(shè)設(shè) , ,對(duì)于第一個(gè)曲邊四邊形有對(duì)于第一個(gè)曲邊四邊形有 CBxAxy 2 CBxAxxf CBxAxxf CBxAxxf 2 2 22 1 2 11 0 2 00 )( )( )( 個(gè)均等的小區(qū)間個(gè)均等的小區(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ù)的形式，給出的被積函數(shù) 是一組離散的數(shù)據(jù)點(diǎn)是一組離散的數(shù)據(jù)點(diǎn) ，此時(shí)三，此時(shí)三 種計(jì)算方法仍舊適用，只需將公式種計(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拋物線法計(jì)算積分拋物線法計(jì)算積分 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 計(jì)算程序計(jì)算程序 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 計(jì)算程序計(jì)算程序 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 計(jì)算程序計(jì)算程序 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、、拋物線法編程：用矩形法、梯形法、拋物線法編程 計(jì)算定積分計(jì)算定積分( ) 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)用 點(diǎn)電荷的電勢(shì)：點(diǎn)電荷的電勢(shì)： r Q k r Q u 0 4 1 （庫侖常數(shù)米（庫侖常數(shù)米.伏伏/庫；真空中介電常數(shù)為庫；真空中介電常數(shù)為 法法/米）米） 9 109k 12 0 1085. 8 1 1電勢(shì)的計(jì)算電勢(shì)的計(jì)算 Q電量為電量為 的點(diǎn)電荷，的點(diǎn)電荷， 處產(chǎn)生的電勢(shì)為：處產(chǎn)生的電勢(shì)為：r在距離在距離 ),(yxPr Q * *

13、 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用14 例例2：設(shè)有一長(zhǎng)直導(dǎo)線均勻帶電，線電荷密度為：設(shè)有一長(zhǎng)直導(dǎo)線均勻帶電，線電荷密度為 ， 解：在導(dǎo)線上取一小段解：在導(dǎo)線上取一小段 ，視為點(diǎn)電荷，其電量為，視為點(diǎn)電荷，其電量為 ，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長(zhǎng)度為長(zhǎng)度為 ，求空間任一點(diǎn)，求空間任一點(diǎn) 的電勢(shì)。的電勢(shì)。 P它在它在 點(diǎn)產(chǎn)生的電勢(shì)為：點(diǎn)產(chǎn)生的電勢(shì)為： 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)（半徑為 ）在軸線上點(diǎn)）在軸線上點(diǎn) 的電勢(shì)的電勢(shì) （線電荷密度（線電荷密度 ） 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：帶電圓盤（半徑為：帶電圓盤（半徑為 ）在軸線上）在軸線上 點(diǎn)的電勢(shì)點(diǎn)的電勢(shì) Rx 2/122 0 2/122 0 )

16、( 2 4 1 )(4 1 rx rdr rx dQ du 環(huán)帶電量環(huán)帶電量 ，在，在 處的電勢(shì)為：處的電勢(shì)為：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 圓盤軸線上電勢(shì)分布圓盤軸線上電勢(shì)分布 練習(xí)：練習(xí)：(1)(1)5 . 0,08. 0,06. 0Rx 用三種數(shù)值積分方法計(jì)算例用三種數(shù)值積分方法計(jì)算例4中的電勢(shì)。中的電勢(shì)。 (2) ( (2) ( 隨隨 變化的曲線變化的曲線) )08. 008. 0 xux FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用19 2電場(chǎng)強(qiáng)度的計(jì)算電場(chǎng)強(qiáng)度的計(jì)算 電量為電量為 的點(diǎn)電荷，在距離的點(diǎn)電荷，在距離 處產(chǎn)生的電場(chǎng)強(qiáng)度為處產(chǎn)生的電場(chǎng)強(qiáng)度為Qr 點(diǎn)電荷的電場(chǎng)強(qiáng)度：點(diǎn)電荷的電場(chǎng)強(qiáng)度： rE 2 0 4 1 r Q )( r r

18、 r FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用20 解：圓環(huán)上任一點(diǎn)電荷解：圓環(huán)上任一點(diǎn)電荷 在軸線上點(diǎn)在軸線上點(diǎn) 的電場(chǎng)強(qiáng)度的電場(chǎng)強(qiáng)度 dsxEd 例例5：求帶電圓環(huán)（半徑為：求帶電圓環(huán)（半徑為 ）在軸線上點(diǎn)）在軸線上點(diǎn) 的電場(chǎng)強(qiáng)度的電場(chǎng)強(qiáng)度 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)軸線上電場(chǎng)強(qiáng)度分布圓環(huán)軸線上電場(chǎng)強(qiáng)度分布 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用22 例例6：已知圓環(huán)帶電量為：已知圓環(huán)帶電量為q q ，桿的線密度為，桿的線密度為 ，長(zhǎng)為，長(zhǎng)為L(zhǎng) L。 求：桿對(duì)圓環(huán)的作用力。求：桿對(duì)圓環(huán)的作用力。 xqdd 解：解： 圓環(huán)在圓環(huán)在 dq 處產(chǎn)生的電場(chǎng)處產(chǎn)生的電場(chǎng) 2/322 0 )(4 1 xR qx Ex x E R q L qd x o 在桿上位置在桿上位置x處取一個(gè)電量為處取一個(gè)電量為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所受電場(chǎng)力為：所受電場(chǎng)力為： 取取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，桿對(duì)環(huán)的作用力又是多少？，桿對(duì)環(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: 求均勻帶電圓盤在軸線上的電場(chǎng)強(qiáng)度。求均勻帶電圓盤在軸線上的電場(chǎng)強(qiáng)度。 為電荷為電荷 E整個(gè)圓盤的電場(chǎng)整個(gè)圓盤的電場(chǎng) 解析解：解析解： )( 1 2 2/122 0 xR x E R面密度，為圓盤半徑。面密度，為圓盤半徑。 dE它在軸線上的電場(chǎng)強(qiáng)度為它在軸線上的電場(chǎng)強(qiáng)度為 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 圓盤軸線上電場(chǎng)強(qiáng)度分布圓盤軸線上電場(chǎng)強(qiáng)度分布 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 （用三種方法，注意（用三種方法，注意 的選取）的選?。?N 用三種數(shù)值積分方法計(jì)算例用三種數(shù)值積分方法計(jì)算例7中的場(chǎng)強(qiáng)中的場(chǎng)強(qiáng) 。E FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用28 EX3-2：對(duì)圓盤軸線上一點(diǎn)：對(duì)圓盤軸線上一點(diǎn)E 2 1 2/322 0 )(2 R R rx rdrx E 采用解析法和數(shù)值法（三種）計(jì)算采用解析法和數(shù)值法（三種）計(jì)算 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：求載流圓線圈軸線上一點(diǎn)：求載流圓線圈軸線上一點(diǎn) P 的磁感應(yīng)強(qiáng)度的磁感應(yīng)強(qiáng)度 2 00 d 4 d r rlI B 2 0 d 4 d r lI B 畢薩定律：畢薩定律： 圖中電流元在圖中電流元在P點(diǎn)產(chǎn)生磁感應(yīng)強(qiáng)度點(diǎn)產(chǎn)生磁感應(yīng)強(qiáng)度dB： )( d 4 22 0 xR lI I lI d B d B d P x R o r x 3磁場(chǎng)

25、的計(jì)算磁場(chǎng)的計(jì)算 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用30 根據(jù)對(duì)稱性根據(jù)對(duì)稱性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)強(qiáng)度分布載流圓線圈軸線上磁感應(yīng)強(qiá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)的帶電圓盤軸線上的磁場(chǎng)：求繞軸旋轉(zhuǎn)的帶電圓盤軸線上的磁場(chǎng) 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點(diǎn)產(chǎn)生的磁感應(yīng)強(qiáng)度點(diǎn)產(chǎn)生的磁感應(yīng)強(qiáng)度dB： x Rx xR dBB R 2 2 2 22 22 0 0 對(duì)上式積分：對(duì)上式積分： 取取=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)的帶電圓盤軸線上的磁場(chǎng)分布繞軸旋轉(zhuǎn)的帶電圓盤軸線上的磁場(chǎng)分布 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、 下氮分子、氧分子運(yùn)動(dòng)的速率分布曲線，并求當(dāng)速下氮分子、氧分子運(yùn)動(dòng)的速率分布曲線，并求當(dāng)速 率在率在300500m/s范圍內(nèi)的分子所占的比例，討論溫范圍內(nèi)的分子所占的比例，討論溫 度度 及分子量及分子量 對(duì)速率分布曲線的影響。對(duì)速率分布曲線的影響。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 對(duì)于氮?dú)猓肿淤|(zhì)量分別為對(duì)于氮?dú)?，分子質(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 氮?dú)猓?T=273K 氧氣， T=400K 氮?dú)猓?T=400K 氣體分子運(yùn)動(dòng)麥克斯韋速率分布曲線氣體分子運(yùn)動(dòng)麥克斯韋速率分布曲線 當(dāng)溫度升高和分

32、子量減小時(shí)，速率小的分子數(shù)減小而速率大的分子當(dāng)溫度升高和分子量減小時(shí)，速率小的分子數(shù)減小而速率大的分子 隨著溫度升高和分子量減小，分布曲線將變得越來越平坦。隨著溫度升高和分子量減小，分布曲線將變得越來越平坦。 氮?dú)夥肿铀急壤秊榈獨(dú)夥肿铀急壤秊?.3957 (T=273K)，0.3085(T=400K)； 氧氣分子所占比例為氧氣分子所占比例為0.4189 (T=273K)，0.3410(T =400K)。 數(shù)增多，分布曲線的極大值隨著溫度升高和分子量減小而向右移動(dòng)。數(shù)增多，分布曲線的極大值隨著溫度升高和分子量減小而向右移動(dòng)。 利用辛普森法可得，當(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)下的麥克斯韋速率分布律，利 用辛普森法計(jì)算用辛普森法計(jì)算0 0C C下氮?dú)?、氧氣分子運(yùn)動(dòng)的平均速率下氮?dú)?、氧氣分子運(yùn)動(dòng)的平均速率 和和 方均根速率方均根速率 ，并與解析解，并與解析解 ， 做比較，其中做比較，其中 是氣體的摩爾質(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其中：其中： 問題：?jiǎn)栴}： 1.二重積分問題？二重積分問題？ 2.不等分間隔問題不等分間隔問題 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用41 3.2 龍貝格法及其應(yīng)用龍貝格法及其應(yīng)用 三種簡(jiǎn)單定積分的計(jì)算方法缺點(diǎn)是要在求積之前給三種簡(jiǎn)單定積分的計(jì)算方法缺點(diǎn)是要在求積之前給 在實(shí)際計(jì)算中，經(jīng)常采用變步長(zhǎng)的計(jì)在實(shí)際計(jì)算中，經(jīng)常采用變步長(zhǎng)的計(jì) 算過程，即在步長(zhǎng)逐步分半的過程中，算過程，即在步長(zhǎng)逐步分半的過程中， 反復(fù)利用求積公式進(jìn)行計(jì)算，直

38、到求反復(fù)利用求積公式進(jìn)行計(jì)算，直到求 得的積分值滿足精度要求為止得的積分值滿足精度要求為止! 出步長(zhǎng)，如何選擇最合適的步長(zhǎng)來統(tǒng)一精度和計(jì)算出步長(zhǎng)，如何選擇最合適的步長(zhǎng)來統(tǒng)一精度和計(jì)算 量的矛盾往往有困難！量的矛盾往往有困難！ FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用42 一、變步長(zhǎng)的梯形法一、變步長(zhǎng)的梯形法 共有共有 個(gè)分點(diǎn)個(gè)分點(diǎn)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用用 表示簡(jiǎn)單梯形法所求得的積分值。表示簡(jiǎn)單梯形法所求得的積

39、分值。 n,ba在簡(jiǎn)單的梯形法中，將求積區(qū)間在簡(jiǎn)單的梯形法中，將求積區(qū)間 分成分成 等份，則等份，則 ),( 2 1 1 2 1 kk k xxx其中點(diǎn)為其中點(diǎn)為, 1kk xx 1、考察一個(gè)子段、考察一個(gè)子段 在該子段上二分前和二分后的兩個(gè)積分值分別為在該子段上二分前和二分后的兩個(gè)積分值分別為 1 T 2 T 和和 FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用43 2、利用這個(gè)、利用這個(gè)關(guān)系式當(dāng)關(guān)系式當(dāng) 從從0到到 累加累加時(shí)，可導(dǎo)出下面時(shí)，可導(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 二分步長(zhǎng)以后的積分值可由二分步長(zhǎng)以前的積分值來計(jì)算二分步長(zhǎng)以后的積分值可由二分步長(zhǎng)以前的積分值來計(jì)算! 遞推公式：遞推公式： )()(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 二、變步長(zhǎng)的辛普森求積法二、變步長(zhǎng)的辛普森求積法 用梯形法二分前后的兩個(gè)積

41、分值用梯形法二分前后的兩個(gè)積分值 和和 作如下線性組合，作如下線性組合， n T n T2 不進(jìn)行證明！不進(jìn)行證明！ 說明說明變步長(zhǎng)梯形法算法簡(jiǎn)單，但精度差，收斂慢變步長(zhǎng)梯形法算法簡(jiǎn)單，但精度差，收斂慢! | 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：我國(guó)第一顆人造地球衛(wèi)星軌道是一個(gè)橢圓，橢圓周長(zhǎng)：我國(guó)第一顆人造地球衛(wèi)星軌道是一個(gè)橢圓，橢圓周長(zhǎng) 2 /2 2 0 41sin c sad a 其中其中 是橢圓的半長(zhǎng)軸，是地球的中心與軌道中心（橢圓中心）是橢圓的半長(zhǎng)軸，是地球的中心與軌道中心（橢圓中心） 的距離，地球半徑的距離，地球半徑 =6371km，近地點(diǎn)距離，近地點(diǎn)距離 =439km， 遠(yuǎn)地點(diǎn)距離遠(yuǎn)地點(diǎn)距離 =2384km, , ac Rh H)2( 2 1 hHRa)( 2 1 hHc 試?yán)米儾介L(zhǎng)辛普森法求

43、解定積分計(jì)算人造地球衛(wèi)星的試?yán)米儾介L(zhǎng)辛普森法求解定積分計(jì)算人造地球衛(wèi)星的 的計(jì)算公式為：的計(jì)算公式為： 軌道周長(zhǎng)。軌道周長(zhǎng)。 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 計(jì)算程序計(jì)算程序 程序運(yùn)行結(jié)果：程序運(yùn)行結(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 三、龍貝格求積法三、龍貝格求積法 不進(jìn)行證明！不進(jìn)行證明！ )( 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é)果的精度更高，可以利用龍貝格法進(jìn)行計(jì)算。為了使求積結(jié)果的精度更高，可以利用

46、龍貝格法進(jìn)行計(jì)算。 n S2 n S用變步長(zhǎng)辛普森法二分前后的兩個(gè)積分值用變步長(zhǎng)辛普森法二分前后的兩個(gè)積分值 和和 作如下的作如下的 n C線性組合，得到的是柯特斯的積分值線性組合，得到的是柯特斯的積分值 重復(fù)積分過程，直到相鄰兩次的積分值重復(fù)積分過程，直到相鄰兩次的積分值 與與 滿足下列關(guān)系：滿足下列關(guān)系： n R2 n R FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用48 例例12：利用龍貝格法計(jì)算單擺的振蕩周期利用龍貝格法計(jì)算單擺的振蕩周期 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 計(jì)算程序：計(jì)算程序： 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 程序運(yùn)行結(jié)果如下：程序運(yùn)行結(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: 若橢圓方程為若橢圓方程為 ，則橢圓周長(zhǎng)的計(jì)算公式為，則橢圓周長(zhǎng)的計(jì)算公式為 利用辛普森法和變步長(zhǎng)辛普森法計(jì)算橢圓利用辛普森法和變步長(zhǎng)辛普森法計(jì)算橢圓 的周長(zhǎng)。的周長(zhǎng)。 EX3-7: 應(yīng)用變步長(zhǎng)辛普森法編程計(jì)算積分應(yīng)用變步長(zhǎng)

50、辛普森法編程計(jì)算積分 此積分值稱為此積分值稱為Catalan常數(shù)，常數(shù)， 的真值為的真值為0.915965。要求誤差。要求誤差 不超過不超過10-5。 1 0 arctan dx x x G G EX3-8:應(yīng)用龍貝格法計(jì)算積分應(yīng)用龍貝格法計(jì)算積分: 要求誤差不超過要求誤差不超過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ù)精度 點(diǎn)，用不同的點(diǎn)，用不同的m次多項(xiàng)式次多項(xiàng)式 0

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

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

53、其在物理學(xué) 中應(yīng)用中應(yīng)用53 v如上圖所示，當(dāng)節(jié)點(diǎn)位置不同時(shí)，積分結(jié)果是不一樣的。如上圖所示，當(dāng)節(jié)點(diǎn)位置不同時(shí)，積分結(jié)果是不一樣的。 圖圖 (a) 一般而言，代數(shù)精度越高，數(shù)值求積公一般而言，代數(shù)精度越高，數(shù)值求積公 式越準(zhǔn)確。要使求積公式具有式越準(zhǔn)確。要使求積公式具有 次代數(shù)精度，次代數(shù)精度， 只要令它對(duì)于只要令它對(duì)于 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ù)精度。 證明：首先驗(yàn)證證明：首先驗(yàn)證xxf, 1)(時(shí)，梯形公式準(zhǔn)確成立。時(shí)，梯形公式準(zhǔn)確成立。 b a bfaf abab abdx)()( 2 11 2 1 b a bfaf ab ba ab abxdx)()( 2 2 )( 2 1 22 再驗(yàn)證，當(dāng)再驗(yàn)證，當(dāng) 2 )(xxf時(shí)，梯形公式能否準(zhǔn)確成立時(shí)，梯形公式能否準(zhǔn)確成立 b a bfaf ab ba abab dxx)()( 2 23 22 33 2 另外可以證明辛普森計(jì)算公式具有另外可以證明辛普森計(jì)算公式具有3階代數(shù)精度。階代數(shù)精度。 的形式，的形式，而而xcc 10 xxf, 1)( 由于任意一個(gè)一次多項(xiàng)式可表示成由

55、于任意一個(gè)一次多項(xiàng)式可表示成 ，梯形公式準(zhǔn)確成立，所以梯形公式具有，梯形公式準(zhǔn)確成立，所以梯形公式具有 一階代數(shù)精度。一階代數(shù)精度。 對(duì)于對(duì)于 梯形公式梯形公式 不成立！不成立！ 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é)點(diǎn)如果節(jié)點(diǎn) 是是次多項(xiàng)式次多項(xiàng)式 )()()()( 210n xxxxxxxxx )(xn的根，并且的根，并且與任意一個(gè)次數(shù)不超過與任意一個(gè)

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

57、的，即上是正交的，即 1 1 1 0)()(dxxpxp nn )( 1 xpn的首項(xiàng)系數(shù)為：的首項(xiàng)系數(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é)點(diǎn)節(jié)點(diǎn) 1 1 0 )()( n k kk xfAdxxf 的的)( 1 xpn1n此時(shí)此時(shí)個(gè)零點(diǎn)就是求積式個(gè)零點(diǎn)就是求積式: : n xxx , , , 10 的節(jié)點(diǎn)的節(jié)點(diǎn) FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及

58、其在物理學(xué) 中應(yīng)用中應(yīng)用57 求積系數(shù)：求積系數(shù)： 1 1 )()( )( dx xxx x A kk k 經(jīng)過計(jì)算得：經(jīng)過計(jì)算得： 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 可以估得高斯求積公式的截?cái)嗾`差為：可以估得高斯求積公式的截?cái)嗾`差為： k x 的選取方法的選取方法 k A及求積系數(shù)及求積系數(shù)節(jié)點(diǎn)節(jié)點(diǎn) FORTRAN數(shù)值方法及其在物理學(xué)數(shù)值方法及其在物理學(xué) 中應(yīng)用中應(yīng)用58 幾個(gè)低階高斯求積公式及其余項(xiàng)：幾個(gè)低階高斯求積公式及其余項(xiàng)： (1) 點(diǎn)）(1 0n ,) 1

59、( 2 1 )( 2 1 xx dx d xp 1)( 1 x p 2 0 A可得：可得： 截?cái)嗾`差：截?cái)嗾`差： )( 3 1 )(ffR 0 0 x節(jié)點(diǎn)節(jié)點(diǎn)，由節(jié)點(diǎn)系數(shù)公式 由節(jié)點(diǎn)系數(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） 點(diǎn)）(2 1n ) 13( 2 1 ) 1( 8 1 )( 222 2 2 2 xx dx d xpxxp3)( 2

60、 1 10 AA 2點(diǎn)求積公點(diǎn)求積公式：式： 1 1 1 0 11 ( )()()() 33 kk k f x dxA f xff 截?cái)嗾`差：截?cái)嗾`差： )( 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點(diǎn)點(diǎn) 高斯求積公式為（不做推導(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