連帶Legendre函數(shù)的高精度計(jì)算

1、第 16 卷第 1 期1999 年 1 月計(jì)算物理CHINESE JOURNAL OF COM PUTATIONAL PHYSICSVol .16 , No.1 Jan .,1999連帶 Legendre 函數(shù)的高精度計(jì)算朱世德 向安平 楊 進(jìn)(成都?xì)庀髮W(xué)院基礎(chǔ)科學(xué)系, 610041)摘要mm給出了計(jì)算復(fù)數(shù)次整數(shù)階連帶 Legendre 函數(shù) P (cos)和 Q (cos)的高精度算法。 計(jì)算結(jié)果表明, 在 486 以上的微機(jī)上, 如采用雙精度, 至少可精確到小數(shù)點(diǎn)后 9 位有效數(shù)字。關(guān)鍵詞 連帶 Legendre 函數(shù) 超幾何級數(shù) Gauss 公式 G auss 連續(xù)分?jǐn)?shù)法中圖分類號O17

2、4.60 引言在空氣動力學(xué)等重要的力學(xué)應(yīng)用領(lǐng)域中, 有時(shí)需要計(jì)算連帶 Legendre 函數(shù)。對中等 v值 ,第一類連帶 Legendre 函數(shù) Pm (cos)可通過超幾何級數(shù)2F 1(a , b ;c;z)和遞推關(guān)系計(jì)算 1 ;對較大的 值 ,級數(shù)算法既收斂緩慢又有嚴(yán)重的舍入誤差。雖然可使用遞推關(guān)系改善其收斂性 ,但是 ,由此將導(dǎo)致嚴(yán)重的上溢或下溢發(fā)生 2 。第二類連帶 Legendre 函數(shù) Qm(cos) 并沒有很好收斂性的級數(shù)算法 ,對小的 和 更是如此 。在本文中, 對參變量 090, 0 m 12 和 Re()-1/2 的變化范圍, 建立了計(jì)算Pm(cos)和 Qm(cos)的高

3、精度算法。對中等 值 ,我們用超幾何級數(shù)及其遞推關(guān)系法計(jì)算Pm (cos),而用 Gauss 連續(xù)分?jǐn)?shù)和 Wronskian 關(guān)系法計(jì)算 Qm (cos)。對小角 , 用對數(shù)表達(dá)式計(jì)算 Qm(cos)。對大的 值, 用漸近展開法計(jì)算 Pm(cos)和 Qm(cos)。1 第一類連帶 Legendre 函數(shù)的算法對中等 值 ,第一類 Legendre 函數(shù)可用超幾何級數(shù)精確計(jì)算1, 3 ,即-m1m2P(cos)=m !tan(/2)2F 1(-, +1 ;m +1;sin (/2)(1)mm(+m +1) -mP(cos)=(-1)(-m +1)P(cos)(2)式(1)和(2)中, 為復(fù)數(shù)

4、, m 為整數(shù)。式(1)中超幾何級數(shù)的定義為a(a +1)b(b +1) zab z22 F1(a ,b ;c;z)=1 +1 ! +2 ! +cc(c +1)a(a +1)(a +k -1)b(b +1)(b +k -1)zk+c(c +1)(c +k-1)k ! + (3)收稿日期:1997-03-19 ;修回日期:1998-07-31朱世德男 52 副教授60計(jì)算物理第 16 卷對小的 ( 2)值, Pm(cos)可直接采用式(1)計(jì)算 3 , 對較大的 值, 式(1)中的超幾何級數(shù)收斂緩慢, 且有嚴(yán)重的舍入誤差。改善超幾何級數(shù)(3)的收斂性, 一個(gè)可能的方法是利用 Legendre 遞

5、推公式使 m 增加到較大值。但是 ,用 Legendre 遞推公式計(jì)算時(shí), 不可避免地要出現(xiàn)嚴(yán)重的上溢或下溢 2 。遞推公式 4Pm+1(cos)+2m cotPm(cos)+(+m )(-m +1)Pm-1(cos)=0(4)中包含了線性和平方系數(shù) ,平方系數(shù)將導(dǎo)致遞推計(jì)算中 Legendre 函數(shù)的快速變化 。為了回避這一困難, 我們使用超幾何級數(shù)的 Gauss 公式(遞推關(guān)系)計(jì)算 ,即 1c(c -1)(z -1)2 F1(a , b ;c -1 ;z)+c c -1 -(2c -a -b -1)z 2 F1(a , b ;c;z)+(c -a)(c -b)z 2 F1(a , b ;

6、c +1 ;z)=0(5)式(5)中, 超幾何級數(shù)的系數(shù)對 c 是同級的 ,因此, 在遞推中,超幾何級數(shù)隨系數(shù)的改變是漸變的。利用基本的方法 ,可以證明超幾何連續(xù)公式的穩(wěn)定性是非常好的 。而且,通過使用式(5),Pm (cos)Qm (cos)能夠?qū)崿F(xiàn)導(dǎo)數(shù) 和 的靈活計(jì)算。數(shù)值計(jì)算表明, 如果 不是太大的話 ,使用式(5)是非常有效的。對大的 , 此法可能會遇到精度丟失的困難 ,采用漸近展開法既能獲得高的精度又能減少耗時(shí)。2 第二類連帶 Legendre 函數(shù)的算法2.1Gauss 連續(xù)分?jǐn)?shù)算法對小的 和, 第二類 Legendre 函數(shù) Qm(cos)并不具有很好收斂性的級數(shù)算法。我們采用

7、Elder 4 計(jì)算第二類 Bessel 函數(shù)的方法, 即 Gauss 連續(xù)分?jǐn)?shù)算法和 Wronskian 關(guān)系。利用Q(z )在割線 -1 x 1 上之值, 能夠證明 1, 3Qm(x)=1e-m i e-m i/ 2 Qm(x +i0)+e mi/2 Qm(x -i0)(6)2是當(dāng)復(fù)變量 z 趨于割線時(shí)的極限。在復(fù)平面上,我們有 1式(6)中, x =cos, x +i0 和 x -i0e-m i Qm(cos+i0)(+3/2)=1/2 (+m +1)emi/2(2sin)me-(+m +1)i 2 F 1(m +1/2 , +m +1 ;+3/2 ;e-2i)(7)用 -1 替換式(7

8、)中的 , 并將結(jié)果除式(7)有2 F1(m +1/2 , +m +1 ;+3/2 ;e-2i)mQm(x +i0)+m-iR(x +i0)=Qmx i0=eFm 1 2m1 2 e-2i(8)+1/22, +-1( + )1( + /;+ / ; )令 a = m +1/2 , b =+m , c =+1/2 , z =e-2i ,則mbe-i2 F1(a , b +1 ;c +1 ;z)R(x +i0)=c2 F 1(a , b ;c;z)(9)式(9)可用 Gauss 連續(xù)分?jǐn)?shù)算法計(jì)算 6 ,即em i/2第 1 期朱世德等:連帶 Legendre 函數(shù)的高精度計(jì)算612 F1(a ,

9、b +1 ;c +1 ;z)=1a(c -b)z2 F 1(a , b ;c;z)-1c(c +1)1 -(b +1)(c -a +1)z(c +1)(c +2)1 -(a +1)(c -b +1)z(c +2)(c +3)1 -(b +2)(c -a +2)z(c +3)(c +4)1 -(10)在復(fù)平面上 ,Gauss 連續(xù)分?jǐn)?shù)處處收斂, 除了在實(shí)軸(+1 x )和2 F1(a , b ;c ;z)的孤立零點(diǎn)外。m 6為了由Gauss 連續(xù)分?jǐn)?shù)獲得 Q(x +i0),使用 W ronskian 關(guān)系dP-m (z )-mme mi-md Qm(z )m式W P(z ), Q(z ) =(1

10、 -z2)=P(z)dz-Q(z)dz(115)和導(dǎo)數(shù)關(guān)系 12dP-m(z)-m-m(z-1)dz=zP(z)-(-m )P-1(z)(12)2dQm(z)mm(z-1)dz=zQ(z )-(+m )Q-1(z)(13)容易求得(+m)P-m (z)Qm-1(z)-(-m)P-m1(z)Qm(z)=em i(14)m式(14)兩邊除以 Q-1(z),并用 cos+i0 代替 z ,我們得到(-m)P-m (cos+i0)-(-m)P-m1(cos+i0)Rm(cos+i0)=emi(15)Qm-1(cos+i0)用 +1 代替式(15)中的 ,并求解之 ,有Qm(cos+i0)em i=(+

11、m +1)P+-m1(cos+i0)-(-m +1)P-m (cos+i0)Rm+1(cos+i0) (16)由于 1P-m (cos+i0)=e mi/2 P-m (cos)(17)所以 ,代式(17)入式(16)中,有Qm cosi0( + )=-m-mm(+m +1)P+1(cos)-(-m +1)P (cos)R+1(cos+i0)(18)62計(jì)算物理第 16 卷由式(18),并使用關(guān)系Qm(cos)=e-3m i/2 Qm(cos+i0)+(i/2)Pm(cos)(19)我們最終得到mQm(cos)=(-1)-m-mm(+m+1)P+1(cos)-(-m +1)P(cos)R+1(c

12、os+i0)+iPm(cos).(20)2雖然隨著 的增加Gauss 連續(xù)分?jǐn)?shù)似乎收斂得更快 ,但是, 上述算法對大的 , 可能會遇到 Pm(cos)精度丟失的困難 。因而,對大的 , 應(yīng)使用其漸近展開式計(jì)算 Pm(cos)和 Qm(cos)。出現(xiàn)困難的另一區(qū)域是在小的 (1)值。此時(shí) ,復(fù)數(shù) Gauss 連續(xù)分?jǐn)?shù)算法仍然是精確的 ,然而卻收斂得很慢,當(dāng)宗量趨于割線 +1 時(shí)更是如此 。采用對數(shù)表達(dá)式來計(jì)算 Qm(cos) 即能回避這一困難。采用實(shí)分析 ,可以建立計(jì)算 Qm(cos)的類似算法 。由于其結(jié)果對整數(shù) 沒有定義, 而且必須計(jì)算收斂性較差的兩個(gè)連續(xù)分?jǐn)?shù), 因此我們僅把其作為檢驗(yàn)上述復(fù)

13、方法收斂性 、精度和速度的一種算法。2.2 對數(shù)表達(dá)式算法當(dāng) 1時(shí), 采用對數(shù)表達(dá)式不僅能獲得同 Gauss 連續(xù)分?jǐn)?shù)算法相同的精度 ,而且耗時(shí)更少。由文 1 ,3 ,可導(dǎo)出如下公式m2 s -1Qm(cos)=1P-m (cos)lnctg2(/2) -2 -2 (+1)+2s =1 (+s)(+1 -s)(-1)mm(m -r -1)!m22222+2ctg(/2)(-)(1 -)(4 -)(r -1) -)(+r)r !r =0m(/2)2222r2rsin(-)(1 -)(m+l -1) -)(m +l )(-1)sin(/2)+2l !(m +l )!l =12 l(-1)m 1 +

14、1/2 +1/3 + +1/ l sin(/2)+2(-m +1)(-m +2)(+m)2222m(-)(1-)(r -l) -)(+r)tan(/2)!(r +m )!1 +1/2 +1/3r=0+ +1/(m +r) sin2r(/2)(21)式(21)中 , 是 Euler 常數(shù), ()是 Psi 函數(shù) 。當(dāng) m =0 , 不是整數(shù)時(shí), 有0102Q(cos)=2P)(cos)lnctg(/P2)-2 -2 (+e1)122222l(-)(1-)(l -1) -)(l +)21 +1/2 +1/3 +1/ l sin(/2) (22)l 1(l !)=當(dāng) m =0 , 是整數(shù) n 時(shí),

15、有y uhrQ0n(cos)=12 P0n(cos)lnctg2(/2) -2 1 +1/2 +1/3 + +1/ nnl1 +1/2 +1/3 + +1/ l (-1)(n +l)!2l(/2) (23)(l2sinl =1!)(n -1)!第 1 期朱世德等:連帶 Legendre 函數(shù)的高精度計(jì)算633 算例按上述算法 ,我們用 Microsoft FORTRAN V5.10 和 VAX FORTRAN V5.0 編制了程序 , 程序中的所有實(shí)變量和復(fù)變量均采用雙精度 。為了驗(yàn)證算法和程序的有效性, 我們在 486 、奔騰微機(jī)和 VAX/VM S 750 計(jì)算機(jī)上進(jìn)行了大量的計(jì)算.在此僅

16、給出在微機(jī)上計(jì)算得到的幾種典型情況的結(jié)果 ,VAX 機(jī)上的計(jì)算結(jié)果用來評價(jià)微機(jī)上計(jì)算結(jié)果的精度。為了進(jìn)一步證實(shí)本算法的精度,用 Maple V Release 4 編制了下述算例的 Maple 程序代碼 , 運(yùn)算結(jié)果證實(shí)了算例能獲得的精度 。本文 FORTRAN 程序中的 值就是取自 M aple 的計(jì)算結(jié)果, 用 Maple 代碼“ ealf (gamma , 64 );” , 立即得到有 64 位有效數(shù)字的 值 , 即0.5772156649015328606065120900824024310421593359399235988057672349 。之所以在此用 M aple 作為驗(yàn)證精

17、度的工具 ,是因?yàn)?Maple 可以實(shí)現(xiàn)任意精度的浮點(diǎn)運(yùn)算。首先計(jì)算 Q 20.0(cos),由式(23)容易導(dǎo)出Q20(cos)=ln ctg2(/2)(3cos2 -1)/4 -(3/2)cos m(24)在 0.0180范圍內(nèi) ,用Gauss 連續(xù)分?jǐn)?shù)算法和式(24)計(jì)算的結(jié)果列于表 1 。在表 1 中,列出了小數(shù)點(diǎn)后 15 位數(shù)字 ,對不同的算法和計(jì)算機(jī)它們是完全相同的 。為了檢驗(yàn)對數(shù)表達(dá)式算法的精度和速度 ,用 Gauss 連續(xù)分?jǐn)?shù)算法和對數(shù)表達(dá)式算法計(jì)算Q10.1 cos 和 Q 110.5 cos(0.005010.000),其結(jié)果列于表 2 。在表2 中 ,列出了小數(shù)點(diǎn)后 14

18、 位數(shù)字,對不同的算法和計(jì)算機(jī)它們也是完全相同的。計(jì)算表明 ,利用對數(shù)表達(dá)式算法能獲得同樣精度的結(jié)果, 但除了小角度外收斂較慢 。表 35101列出了 為純實(shí)數(shù)時(shí) Q 0.5 和 Q 0.5 的計(jì)算值, 表 4列出了錐函數(shù) P -0.5+i y 0.5 、P10-0.5 +i y 0.5 和 Q1-0.5+i y 0.5 、Q10-0.5+i y 0.5 的計(jì)算值 。在表 3 和 4 中僅列出了小數(shù)點(diǎn)后 10 位數(shù)字。表 1 Q20.0 cos 計(jì)算結(jié)果的比較Table 1The comparisons of the computation results for Q20.0 cos0Com

19、plex GAUSS ContinuedEQ(24)Fration Method0.010.784654392482865D +010.784654392482865D +010.100.554392908372049D +010.554392908372049D +010.500.393395131812442D +010.393395131812442D +011.000.323941099146171D +010.323941099146171D +0110.00.848841713232357D +000.848841713232357D +0020.00.21368715888554

20、0D -010.213687158885540D -0130.0-0.475939420098648D +00-0.475939420098648D +0040.0-0.764768397072655D +00-0.764768397072655D +0050.0-0.872812404617566D +00-0.872812404617566D +0060.0-0.818663268041757D +00-0.818663268041757D +0070.0-0.628686918700802D +00-0.628686918700802D +0080.0-0.340250577301736

21、D +00-0.340250577301736D +0064計(jì)算物理第 16 卷表 2 用 Gauss 連續(xù)分?jǐn)?shù)和對數(shù)表達(dá)式算法計(jì)算 Q10.1 cos 和 Q110.5 cos 的結(jié)果 Table 2 Calculated results of Q10.1 cos and Q110.5 cos by using Gauss fraction and logarithmic expression algorithmQ 1 cos()Q 1 cos()0.110.50.005-0.11459155967011D +05-0.11459155967011D +050.010-0.572957807

22、34440D +04-0.57295780734440D +040.050-0.11459161236929D +04-0.11459161236929D +040.100-0.57295879545620D +03-0.57295879545620D +030.250-0.22918539897471D +03-0.22918539897471D +030.500-0.11459578820329D +03-0.11459578820359D +030.750-0.76400424582516D +02-0.76400424582516D +021.000-0.57303572639239D

23、 +02-0.57303572639239D +025.000-0.11490416444659D +02-0.11490416444659D +0210.00-0.578555544318911D +01-0.578555544318911D +01表 3 用 Gauss 連續(xù)分?jǐn)?shù)和對數(shù)表達(dá)式算法計(jì)算 Q5 0.5 和 Q10 0.5 的結(jié)果 Table 3Calculated results of Q5 0.5and Q10 0.5 by using Gauss fractionand logarithmic expression algorithmQ50.5Q10 0.50.100.18

24、90080188D +030.4422019441D +080.500.1954423849D +030.4499152949D +081.00.2093856976D +030.4654023111D +081.500.2338751664D +030.4880508920D +082.000.2709697263D +030.5189632000D +082.500.3210500676D +030.5597232144D +083.000.3941377838D +030.6124885333D +085.000.1850243061D +040.1024546133D +0910.00

25、0.2965453943D +040.1932505391D +1015.000.2840985316D +060.1433781264D +1220.000.6037537129D +060.5535387943D +1225.000.1358738794D +070.2316748892D +1430.000.6349051171D +070.9582603845D +1440.000.1307875271D +080.2309968419D +1650.000.2432863913D +080.1297283270D +17表 4錐函數(shù) Pm 0.5 Qm 0.5 的計(jì)算結(jié)果-0.5+i

26、y-0.5 +iyTable 4Calculated results of Pm-0.5 +iy 0.5 and Qm-0.5 +i y 0.5yP 1-0.5 +iy 0.5P10-0.5+iy 0.5Q1-0.5+iy 0.5Q 10-0.5 +iy 0.50.10.1555200727 D +000.4891358513D +030.7762014449 D +000.4378198940D +080.50.3088993114 D +000.1147841024D +040.8309762935 D +000.4349983194D +081.00.8524933343 D +000.

27、5007038619D +040.1439505070 D +010.4263082576D +081.50.1997326275 D +010.2185782502D +050.3159850134 D +010.4122435524D +082.00.4187708476 D +010.9207521121D +050.6582936938 D +010.3933997825D +082.50.8246410482 D +010.3739295323D +060.1295450012 D +020.3705817184D +083.00.1565137984 D +020.14652184

28、10D +070.2458536301 D +020.3453079918D +085.00.1714349478 D +030.2469554929D +090.2692893868 D +030.3885844271D +0910.00.4678102294 D +050.1372346802D +140.7348345900 D +050.2155677315D +1415.00.1085296078 D +080.1288283463D +180.1704779092 D +080.2023630932D +1820.00.2364063884 D +100.4336908303D +

29、210.3713462863 D +100.6812399633D +2125.00.4978136689 D +120.7790426050D +240.7819418950 D +120.1223719727D +2530.00.1026296081 D +150.9252636097D +270.1081704866 D +150.1844475912D +2840.00.4192810395 D +190.5908936634D +330.6585579927 D +190.9281735972D +3350.00.1657040143 D +240.193827364D +390.2

30、602872570 D +240.3044633118D +39第 1 期朱世德等:連帶 Legendre 函數(shù)的高精度計(jì)算654 結(jié)論與討論我們導(dǎo)出了計(jì)算復(fù)數(shù)次整數(shù)階連帶 Legendre 函數(shù)精確而有效的算法。通常, 在微機(jī)上計(jì)算 Pm(cos)和 Qm(cos)的最高精度近似為采用雙精度運(yùn)算所能獲得的精度(15 位有效數(shù)字)。大量計(jì)算表明 ,對參變量 090,0 m 12 和 Re()-1/2 的變化范圍, 本文提出的算法精度在微機(jī)上至少精確到小數(shù)點(diǎn)后 9 位有效數(shù)字 。在參變量 090,0 m 12 和 Re()-1/2 的變化范圍之外, 算法的精度將降低 。本文的 M aple 數(shù)值計(jì)算表明,M aple 在高精度數(shù)值計(jì)算方面遠(yuǎn)遠(yuǎn)勝過 FORTRAN 。其最大缺點(diǎn)是計(jì)算效率低 ,通常 Maple 程序比相應(yīng)的 FORT RAN 程序慢 1000 倍。參考文獻(xiàn)1Erdelyi A , et al.Higher T ranscendental Functions.Vol.1, Bateman M anuscript Project, M cGraw Hill BookCo ., Inc ., New York , 1953.2Smith J M , et al.Legendre Polynomials.A CM T ransactions on M ath