馮佳計(jì),賈曉偉,沈建琪

(上海理工大學(xué)理學(xué)院,上海200093)

?

小變量情況下第一類整數(shù)階Bessel函數(shù)的計(jì)算

馮佳計(jì),賈曉偉,沈建琪

(上海理工大學(xué)理學(xué)院,上海200093)

在計(jì)算第一類整數(shù)階Bessel函數(shù)時(shí),后向遞推算法穩(wěn)定高效.然而,起始點(diǎn)的選取必須有足夠高的階數(shù),并且需要進(jìn)行歸一化處理.本文對(duì)Taylor級(jí)數(shù)展開(kāi)算法進(jìn)行研究,并對(duì)其級(jí)數(shù)展開(kāi)規(guī)律、計(jì)算精度,以及求和項(xiàng)與參數(shù)間的關(guān)系進(jìn)行了討論.此外,本文利用指數(shù)形式,極大擴(kuò)展了該算法的可計(jì)算范圍.與du Toit算法、MATLAB和Mathematica應(yīng)用軟件的計(jì)算結(jié)果比較顯示,本文的算法具有較高的準(zhǔn)確性和穩(wěn)定性.

Bessel函數(shù);Taylor級(jí)數(shù)展開(kāi);指數(shù)擴(kuò)展

1 引言

第一類整數(shù)階Bessel函數(shù)Jn(z)在物理學(xué)、力學(xué)、電磁學(xué)、環(huán)境科學(xué)等領(lǐng)域都有著重要的用途[1,2].由于應(yīng)用廣泛,大量文獻(xiàn)致力于討論它的性質(zhì)并給出了多種不同的算法[3～16].其中后向遞推算法選取兩個(gè)連續(xù)的高階函數(shù)作為起始點(diǎn)賦任意值并后向遞推計(jì)算,最后做歸一化處理[9～11].這種算法穩(wěn)定且高效.然而,為了確保計(jì)算的順利進(jìn)行,后向遞推算法起始點(diǎn)的階數(shù)必須足夠高,需要更多的計(jì)算機(jī)資源和更大的計(jì)算量.在采用后向遞推算法計(jì)算Bessel函數(shù)時(shí),可考慮采用Taylor級(jí)數(shù)展開(kāi)法精確計(jì)算出起始階數(shù)的函數(shù)值以提高計(jì)算效率.原則上,Bessel函數(shù)的Taylor級(jí)數(shù)展開(kāi)對(duì)所有變量均收斂.然而,其收斂速度卻依賴于階數(shù)n和變量z的大小.在階數(shù)n很小、參數(shù)z很大的情況下,需要計(jì)算的級(jí)數(shù)展開(kāi)項(xiàng)增多,且在計(jì)算過(guò)程中有可能會(huì)丟失有效數(shù)字,導(dǎo)致計(jì)算效率降低、計(jì)算結(jié)果出現(xiàn)偏差.

2 遞推算法

整數(shù)階Bessel函數(shù)滿足前向遞推關(guān)系

(1)

或者后向遞推關(guān)系

(2)

對(duì)于給定參數(shù)z,前向遞推從兩個(gè)連續(xù)低階級(jí)數(shù)開(kāi)始計(jì)算,而后向遞推從兩個(gè)連續(xù)高階級(jí)數(shù)開(kāi)始計(jì)算.像大多數(shù)遞推過(guò)程一樣,上述公式易受傳播誤差的影響[7].

當(dāng)z為實(shí)數(shù)或者虛部很小(|lmz|≤1.0)的復(fù)數(shù)時(shí),Bessel函數(shù)的絕對(duì)值|Jn(z)|在n<0.9|z|范圍變化不大,而在n≥0.9|z|范圍隨階數(shù)的增大快速衰減.對(duì)于n<0.9|z|的范圍,前向遞推或者后向遞推都能用來(lái)計(jì)算Jn(z);但當(dāng)n>0.9|z|時(shí),則只能采用后向遞推法[5,9].如果z的虛部較大,Bessel函數(shù)的絕對(duì)值|Jn(z)|在n≥0范圍內(nèi)隨階數(shù)的增大迅速衰減,只能采用后向遞推法計(jì)算.

用后向遞推方法計(jì)算Bessel函數(shù)Jn(z)需從兩個(gè)連續(xù)的高階級(jí)數(shù)Bq(z)和Bq+1(z)開(kāi)始,Bq(z)和Bq+1(z)的值可以任意選取,du.Toit建議Bq(z)=1和Bq+1(z)=0[9].然而,階數(shù)q必須同時(shí)滿足env|Jq+1(z)|≤10-p/2env|JN(z)|和env|2Jq(z)|≤10-penv|J0(z)|以確保收斂,其中p為所需的有效位數(shù)、N是所需計(jì)算的最高階數(shù),env|x|表示函數(shù)的包絡(luò)線.如果N≤|z|,q的最小值需滿足

(3)

若N>|z|,q值的選取則相對(duì)繁瑣[5,11].完成Bn(z)的計(jì)算后,需進(jìn)行歸一化處理得到相應(yīng)的函數(shù)值Jn(z).針對(duì)不同的變量值,采用不同的歸一化方法.當(dāng)|z|≤25時(shí),Bn(z)歸一化處理為

(4)

當(dāng)|z|>25時(shí),使用S=B0(z)/J0(z)或S=B1(z)/J1(z)歸一化效率更高.使用時(shí)應(yīng)避開(kāi)J0(z)和J1(z)的零點(diǎn).

從以上討論可知,后向遞推法從兩個(gè)任意取值Bq(z)和Bq+1(z)開(kāi)始計(jì)算,不僅需要計(jì)算起始階數(shù)q以及從階數(shù)q到階數(shù)N的函數(shù)Bn(z),還需對(duì)Bn(z)做歸一化處理.因此,為了使計(jì)算更加高效,在使用后向遞推方法時(shí),建議采用Taylor級(jí)數(shù)展開(kāi)法預(yù)先計(jì)算Bq(z)和Bq+1(z)的精確值.

3 Bessel函數(shù)的Taylor級(jí)數(shù)展開(kāi)

整數(shù)階Bessel函數(shù)Jn(z)的Taylor級(jí)數(shù)展開(kāi)式為:

(5)

令

(6)

(7)

則有遞推關(guān)系:

(8)

(9)

(10)

對(duì)于復(fù)變量的一般情況z=|z|eiδz,Taylor級(jí)數(shù)展開(kāi)形式寫成如下形式:

(11)

(12)

(13)

(14)

公式(13)中的求和項(xiàng)在計(jì)算過(guò)程中出現(xiàn)正負(fù)交替的情況,如果kpeak比較大,必然導(dǎo)致數(shù)值比較大的數(shù)字相減并導(dǎo)致有效數(shù)字丟失從而計(jì)算結(jié)果變差.反之,公式(14)中的求和項(xiàng)一直為正,不會(huì)出現(xiàn)公式(13)的情況.因此,在這種情況下,即使kpeak比較大,也能得到合理的計(jì)算結(jié)果,但求和項(xiàng)比較多、收斂較慢.

(15)

在采用Taylor級(jí)數(shù)展開(kāi)式時(shí),并不一定要求嚴(yán)格滿足(|z|/2)2

相位角δz=0時(shí)有效數(shù)字丟失最嚴(yán)重,因此在該相位角情況下滿足S-T<18的所有階數(shù)ncal對(duì)所有的相位角計(jì)算均是可信的(見(jiàn)圖5).可以看出,ncal≈|z|2/40-10表示滿足n>ncal時(shí),對(duì)于任意相位角δz的計(jì)算結(jié)果保證至少7位有效數(shù),這與n>|z|2/4-1的要求降低了大約10倍.

圖6給出了|z|=600時(shí)階數(shù)n和相位角δz的關(guān)系,以曲線為界,上部區(qū)域表示Taylor級(jí)數(shù)展開(kāi)式適用的范圍,下部區(qū)域中計(jì)算得到的結(jié)果其有效位數(shù)低于7位.顯然,如果提高對(duì)有效位數(shù)的要求,則圖中的曲線將會(huì)向上移動(dòng),對(duì)應(yīng)的可計(jì)算范圍縮小.

在20<|z|<104范圍內(nèi)對(duì)多個(gè)不同的復(fù)變量z進(jìn)行了數(shù)值計(jì)算,并進(jìn)行曲線擬合,得到如下經(jīng)驗(yàn)公式:

(16)

a=ln(|z|-0.5232|z|0.3+20.49)-2.997

(17)

(18)

4 復(fù)數(shù)的大數(shù)計(jì)算方法

在復(fù)變量情況下,Bessel函數(shù)值及其Taylor級(jí)數(shù)展開(kāi)式的求和項(xiàng)在一個(gè)很大的動(dòng)態(tài)范圍內(nèi)變化.為了保證計(jì)算的順利進(jìn)行,需要對(duì)中間計(jì)算量的數(shù)值范圍進(jìn)行擴(kuò)展.為此,所有變量均采用復(fù)數(shù)形式儲(chǔ)存并參與運(yùn)算.

(19)

5 數(shù)值計(jì)算結(jié)果和分析

為了驗(yàn)證本文提出的算法,我們將所得結(jié)果與du

表1 Jn(z)數(shù)據(jù)比較

6 結(jié)論

本文介紹了第一類整數(shù)階Bessel函數(shù)的Taylor級(jí)數(shù)展開(kāi)算法.為了保證高速收斂并且具有較高的計(jì)算精度,Taylor級(jí)數(shù)展開(kāi)應(yīng)在n≥|z|2/4-1條件下使用;在保證7位有效數(shù)計(jì)算精度的情況下,給出了Taylor級(jí)數(shù)展開(kāi)式適用范圍的經(jīng)驗(yàn)公式;此外,在算法中采用了指數(shù)形式以擴(kuò)展計(jì)算數(shù)值范圍.計(jì)算結(jié)果與現(xiàn)有的軟件進(jìn)行了比較,表明該算法具有較好的效果.

遞推算法與級(jí)數(shù)展開(kāi)算法結(jié)合使用可有效提高計(jì)算效率.此外,本文中給出的算法也可以應(yīng)用于計(jì)算非整數(shù)階Bessel函數(shù)以及改進(jìn)的Bessel函數(shù).

[1]G Gouesbet,G Gréhan.Generalized Lorenz-Mie Theories[M].Springer-Verlag Berlin Heidel-berg,2011.

[2]JA Lock.Improved Gaussian beam-scattering algorithm[J].Applied Optics,1995,34 (3):599-570.

[3]GN Watson.A Treatise on the Theory of Bessel Functions[M].Cambridge University Press,Cambridge,1958.

[4]M Abramowitz,IA Stegun.Handbook of Mathematical Functions[M].10th ed.,National Bureau of Standards,Washington,DC,1972.

[5]S Zhang,J Jin.Computation of Special Functions[M].Wiley,New York,1996.

[6]M Goldstein,RM Thaler.Recurrence techniques for the calculation of Bessel functions[J].Mathematical Tables & Other Aids to Computation,1959,13(66):102-102.

[7]W Gautschi.Computational aspects of three-term recurrence relations[J].Society for Industrial and Applied Mathematics,1967,9(1):24-28.

[8]DE Amos.Algorithm 644:A portable package for Bessel functions of a complex argument and nonnegative order[J].Acm Transactions on Mathematical Software,1986,12:265-273.

[9]CF du Toit.The numerical computation of Bessel functions of the first and second kind for integer orders and complex arguments[J],IEEE Antennas & Propagation Magazine,1990,38(9):1341-1349.

[10]CF du Toit.Bessel functions Jn(z) and Yn(z) for integer order and complex argument[J].Computer Physics Communications,1993,78(1-2):181-189.

[11]CF du Toit.Evaluation of some algorithms and programs for the computation of integer-order Bessel functions of the first and second kind with complex arguments[J].IEEE Antennas & Propagation Magazine,1993,35(3):19-25.

[12]KG Valeev,OY Kostinskii.Calculation of Bessel functions by using continued fractions[J].Ukrainian Mathematical Journal,1997,47(12):1949-1950.

[13]HA Yousif,R Melka.Bessel function of the first kind with complex argument[J].Computer Physics Communications,1997,106(3):199-206.

[14]EJ Rothwell.Computation of the logarithm of Bessel functions of complex argument[J].Communications in Numerical Methods in Engineering,2005,21(10):597-605.

[15]LW Cai.On the computation of spherical Bessel functions of complex arguments[J].Computer Physics Communications,2011,182(3):663-668.

[16]張善杰,唐漢.任意實(shí)數(shù)階復(fù)宗量第一類和第二類Bessel函數(shù)的精確計(jì)算[J].電子學(xué)報(bào),1996,24(3):77-80.

Zhang Shanjie,Tang Han.Accurate computation of Bessel functions of the first and second kinds with arbitrary real order and complex argument[J].Acta Electronica Sinica,1996,24(3):77-80.(in Chinese)

馮佳計(jì) 男，碩士研究生，1990年5月生于江西省撫州市東鄉(xiāng)縣.主要從事光學(xué)測(cè)試方面的研究.

E-mail：jiaji19900505@163.com

沈建琪(通信作者) 男，1965年7月出生，浙江桐鄉(xiāng)人，博士生導(dǎo)師，分別于1987年和1990年在華東師范大學(xué)獲得理學(xué)學(xué)士和碩士學(xué)位，2003年在德國(guó)Cottbus工大獲得工學(xué)博士學(xué)位.主要從事光散射顆粒測(cè)試技術(shù)研究.

E-mail:jqshenk@163.com

Computation of the Integer Order Bessel Functions of First Kind with Small Arguments

FENG Jia-ji,JIA Xiao-wei,SHEN Jian-qi

(CollegeofScience,UniversityofShanghaiforScienceandTechnology,Shanghai200093,China)

Algorithm based on the backward recurrence for computing the integer order Bessel functions of the first kind is stable and efficient.However,the orders of the starting points should be high enough and the normalization is required.In this paper,we introduce an algorithm based on the Taylor series expansion (TSE),in which all the quantities involved are expressed in the exponential format so as to expand the numeric range of calculation.Comparison against du Toit’s algorithm as well as MATLAB and Mathematica shows that our algorithm is stable and reliable.

Bessel function;Taylor series expansion;exponential scaling

2015-11-24;

2016-02-25;責(zé)任編輯：藍(lán)紅杰

國(guó)家自然科學(xué)基金(No.NSFC 51476104)

TP301.6

A

0372-2112 (2016)11-2720-06

??學(xué)報(bào)URL:http://www.ejournal.org.cn

10.3969/j.issn.0372-2112.2016.11.022