整數(shù)階復(fù)宗量貝塞爾函數(shù)的計(jì)算程序研究

任宏紅 郭迎春 王兵兵

摘要：鑒于目前的算法程序集中沒有現(xiàn)成的計(jì)算復(fù)宗量貝塞爾函數(shù)的程序，本文基于貝塞爾函數(shù)的逆向遞推關(guān)系編寫了計(jì)算整數(shù)階復(fù)宗量第一類貝塞爾函數(shù)的Fortran程序源代碼.與Matlab軟件的計(jì)算結(jié)果比較，兩者至少有12位有效數(shù)字一致.接著運(yùn)用此程序，分析了徐士良的《FORTRAN常用算法程序集》中的純虛宗量的貝塞爾函數(shù)，即變形貝塞爾函數(shù)程序的準(zhǔn)確度，發(fā)現(xiàn)其準(zhǔn)確度為6位有效數(shù)字.最后，對(duì)基于實(shí)宗量貝塞爾函數(shù)和純虛宗量貝塞爾函數(shù)相乘然后用無限求和來計(jì)算復(fù)宗量貝塞爾函數(shù)值的方法的準(zhǔn)確性進(jìn)行了探討.證明其僅能對(duì)有限的貝塞爾函數(shù)進(jìn)行準(zhǔn)確計(jì)算.這是由于當(dāng)求和項(xiàng)中有遠(yuǎn)大于最終的求和項(xiàng)時(shí)，會(huì)導(dǎo)致求和結(jié)果的有效數(shù)字減少甚至完全錯(cuò)誤.

關(guān)鍵詞：復(fù)宗量貝塞爾函數(shù)：Fortran源程序：逆遞推計(jì)算

中圖分類號(hào)：0411.2 文獻(xiàn)標(biāo)志碼：A DOI：10.3969/j.issn.1000-5641.2019.01.009

0引言

整數(shù)階貝塞爾函數(shù)是許多輻射、散射和波導(dǎo)問題的數(shù)學(xué)解，貝塞爾函數(shù)中復(fù)宗量以及虛宗量和材料損耗、漏波等相關(guān).貝塞爾函數(shù)也應(yīng)用于電磁場(chǎng)與物質(zhì)相互作用的相關(guān)現(xiàn)象的數(shù)學(xué)描述，如閾上電離及高次諧波等.所以貝塞爾函數(shù)的精確計(jì)算尤為重要.

針對(duì)整數(shù)階第一類貝塞爾函數(shù)，在常用算法程序書中，例如徐士良的《FORTRAN常用算法程序集》以及《Numerical Recipes：The Art of Scientific Computing》等已經(jīng)有實(shí)宗量的貝塞爾函數(shù)的程序和純虛宗量的貝塞爾函數(shù)即變形的貝塞爾函數(shù)的程序，但是還沒有計(jì)算復(fù)宗量貝塞爾函數(shù)的程序.Du Toit探討了復(fù)宗量的貝塞爾函數(shù)的計(jì)算方法，提出了采用逆向貝塞爾函數(shù)遞推關(guān)系進(jìn)行計(jì)算，文中給出了一些典型的數(shù)值結(jié)果，指出這種方法對(duì)于大的宗量不能計(jì)算，而對(duì)于大的宗量，采用冪級(jí)數(shù)展開截?cái)嗟姆椒ńo出了結(jié)果.魏彥玉等人采用級(jí)數(shù)展開方法對(duì)虛宗量的貝塞爾函數(shù)進(jìn)行了計(jì)算，對(duì)小的宗量計(jì)算結(jié)果很好.張爽等人針對(duì)這種方法的計(jì)算進(jìn)行了誤差分析，指出了此方法不適合大宗量貝塞爾函數(shù)計(jì)算的原因.張善杰等人對(duì)任意實(shí)數(shù)階復(fù)宗量貝塞爾函數(shù)的遞推算法進(jìn)行了研究，并探討了如何驗(yàn)證程序的正確性和分析了計(jì)算結(jié)果的精度.

本文針對(duì)第一類整數(shù)階復(fù)宗量貝塞爾函數(shù)，重新考察了基于逆遞推公式對(duì)貝塞爾函數(shù)的計(jì)算，編寫了適用于任意大小的復(fù)宗量的第一類整數(shù)階貝塞爾函數(shù)Fortran源代碼.本程序計(jì)算的貝塞爾函數(shù)的結(jié)果與Matlab軟件結(jié)果進(jìn)行比較，二者一致到第12位有效數(shù)字.由于常用算法程序書（如文獻(xiàn)）中，很容易找到實(shí)宗量的貝塞爾函數(shù)程序和純虛宗量貝塞爾函數(shù)的程序，基于二者相乘求和，原則上可以計(jì)算復(fù)宗量的貝塞爾函數(shù).本文探討了這種方法計(jì)算貝塞爾函數(shù)的可行性.在這之前還探討了算法程序集中變形貝塞爾函數(shù)的計(jì)算精度.

本文安排如下，第1節(jié)簡(jiǎn)要闡述復(fù)宗量貝塞爾函數(shù)的理論基礎(chǔ)，并舉例比較此程序的結(jié)果與Matlab的結(jié)果;第2節(jié)運(yùn)用給出的程序分析變形貝塞爾函數(shù)程序的準(zhǔn)確度;第3節(jié)分析基于實(shí)宗量貝塞爾函數(shù)和虛宗量貝塞爾函數(shù)相乘無窮求和的方法，計(jì)算復(fù)宗量貝塞爾函數(shù)的可行性;附錄中給出源代碼.

1計(jì)算復(fù)宗量的貝塞爾函數(shù)的原理

文獻(xiàn)指出，編程計(jì)算的誤差有舍入誤差（round off error）、截止誤差（truncationerror）和計(jì)算的穩(wěn)定度誤差（stability error）.舍入誤差是由計(jì)算機(jī)精度所決定的，如雙精度型的變量含16位有效數(shù)字，誤差是最后一位有效數(shù)字的量級(jí).此數(shù)據(jù)精度所引入的誤差為舍入誤差.舍入誤差一般會(huì)隨計(jì)算次數(shù)的增加而增加，一般會(huì)影響結(jié)果的最后兩位有效數(shù)字.在某些特殊情況下會(huì)引入很大的誤差，從而減少有效數(shù)字.如兩個(gè)幾乎相等的同型數(shù)據(jù)相減.截止誤差是由計(jì)算過程中將無窮求和截取為有限求和而引入的誤差.穩(wěn)定度誤差是由于計(jì)算方法不當(dāng)（稱為不穩(wěn)定誤差）使最初階段的舍入誤差被連續(xù)放大而導(dǎo)致的誤差，這種誤差甚至?xí)蜎]真實(shí)結(jié)果.減小這三種誤差是我們程序設(shè)計(jì)的出發(fā)點(diǎn)，文獻(xiàn)已證實(shí)我們所采用的這種逆遞推方法是穩(wěn)定的.不同情況下選取不同的s，目的是減小舍入誤差.q的選取決定了截止誤差，只要q值取得足夠大，截止誤差就會(huì)變小，直到小于舍入誤差.具體表現(xiàn)為，q增大到某值后，如果再繼續(xù)增大，計(jì)算的結(jié)果將保持不變.據(jù)此，為計(jì)算準(zhǔn)確又節(jié)省時(shí)間，我們將q取值為可見，貝塞爾函數(shù)的階數(shù)越大，宗量的模越大，需要的q值就越大.

計(jì)算結(jié)果表明，我們的結(jié)果和Matlab軟件的結(jié)果一致到至少12位有效數(shù)字，另外，我們計(jì)算了文獻(xiàn)表1中的幾個(gè)典型數(shù)據(jù)，再現(xiàn)了文獻(xiàn)的全部結(jié)果，并給出了文獻(xiàn)中沒有計(jì)算出的結(jié)果，列在表1中.對(duì)于小的宗量，小的階數(shù)，由公式（8）可見，q值會(huì)相對(duì)很小，此時(shí)不僅q所決定的截止誤差小于舍入誤差，由于q相對(duì)小，計(jì)算的次數(shù)少，和大階數(shù)大宗量的情況相比，舍入誤差也不會(huì)大，所以計(jì)算的結(jié)果應(yīng)該更為精確.如表1中的前兩行所示，我們計(jì)算的結(jié)果與Matlab的結(jié)果一致到14位有效數(shù)字.文獻(xiàn)強(qiáng)調(diào)這種基于逆向遞推公式的方法不適合計(jì)算大的宗量，認(rèn)為計(jì)算如10 000 000+i333大宗量的貝塞爾函數(shù)是不可行的，該文獻(xiàn)提出運(yùn)用冪級(jí)數(shù)展開截?cái)嗟姆椒▉碛?jì)算大宗量貝塞爾函數(shù)，而我們就是采用逆遞推公式的方法，將如10 000 000+i333這樣大宗量的貝塞爾函數(shù)計(jì)算了出來.得到的結(jié)果與該文獻(xiàn)中基于展開截?cái)嗟姆椒ǖ慕Y(jié)果基本一致（見表1）.我們認(rèn)為文獻(xiàn)中采用這種逆向遞推的方法沒有計(jì)算出來的原因，一方面是我們的逆遞推起點(diǎn)q值比文獻(xiàn)中的大，還有可能是計(jì)算機(jī)的性能的局限造成的.

2變形貝塞爾函數(shù)的準(zhǔn)確度

文獻(xiàn)中給出了計(jì)算變形貝塞爾函數(shù)（即純虛數(shù)的貝塞爾函數(shù)）的程序，下面將采用第l節(jié)的程序來探討此變形貝塞爾函數(shù)的準(zhǔn)確度.

可見程序集中Io（y）的結(jié)果準(zhǔn)確到第6位有效數(shù)字.進(jìn)而確定程序集中的純虛數(shù)的貝塞爾函數(shù)即變形的貝塞爾函數(shù)In（y）的準(zhǔn)確度為6個(gè)有效數(shù)字，遠(yuǎn)小于第1節(jié)提供的12位有效數(shù)字.

3無窮求和計(jì)算貝塞爾函數(shù)

由于程序集中已經(jīng)提供了實(shí)宗量的貝塞爾函數(shù)的程序以及純虛宗量的貝塞爾函數(shù)的程序，很自然想到由下面的無窮求和（式（13））來計(jì)算復(fù)宗量的貝塞爾函數(shù)：所以這一節(jié)，我們將探討基于這種無窮求和來計(jì)算復(fù)宗量貝塞爾函數(shù)的準(zhǔn)確度.

通過與第1節(jié)中提供的程序計(jì)算結(jié)果的比較發(fā)現(xiàn)，采用無窮求和來計(jì)算的結(jié)果僅僅在有限范圍內(nèi)是正確的.一般來講，對(duì)于確定的自變量，階數(shù)變大時(shí)，計(jì)算結(jié)果將會(huì)不準(zhǔn)確.我們定義準(zhǔn)確的結(jié)果指的是計(jì)算結(jié)果的有效數(shù)字與第1節(jié)提供的程序的結(jié)果達(dá)到前6位一致.表3給出了幾個(gè)特定自變量，無窮求和方法能給出正確結(jié)果的貝塞爾函數(shù)的階數(shù)范圍，如：對(duì)于Jn（50+i100），當(dāng)一21

無窮求和的方法會(huì)給出不好的結(jié)果的原因在于式子（13）左邊的貝塞爾函數(shù)值小于求和中的個(gè)別的項(xiàng)值，從而減少了最終求和中準(zhǔn)確的有效數(shù)字的個(gè)數(shù).例如

4結(jié)論

我們基于逆向的貝塞爾函數(shù)的遞推關(guān)系，給出了計(jì)算第一類復(fù)宗量整數(shù)階貝塞爾函數(shù)的Fortran源代碼，探討了現(xiàn)有的算法程序集中的純虛數(shù)的貝塞爾函數(shù)的準(zhǔn)確度為6個(gè)有效數(shù)字，針對(duì)運(yùn)用算法程序集中的實(shí)宗量的貝塞爾函數(shù)和純虛宗量貝塞爾函數(shù)的程序，采用它們的乘積對(duì)階數(shù)的無窮求和來計(jì)算復(fù)宗量的貝塞爾函數(shù)的方法，討論了其計(jì)算的準(zhǔn)確度，通過與本文所給程序的結(jié)果比較，得出其僅能對(duì)有限的貝塞爾函數(shù)進(jìn)行計(jì)算.這是由于無窮求和項(xiàng)中有大于最終求和的項(xiàng)而導(dǎo)致結(jié)果的準(zhǔn)確的有效數(shù)字減少甚至結(jié)果完全錯(cuò)誤，這也是物理量的求和計(jì)算中需要注意的問題.