整数阶复宗量贝塞尔函数的计算程序研究

引用本文

任宏红, 郭迎春, 王兵兵. 整数阶复宗量贝塞尔函数的计算程序研究[J]. 华东师范大学学报(自然科学版), 2019, (1): 76-82, 92. DOI: 10.3969/j.issn.1000-5641.2019.01.009.

REN Hong-hong, GUO Ying-chun, WANG Bing-bing. Program for calculating the integer order of Bessel functions with complex arguments[J]. Journal of East China Normal University (Natural Science), 2019, (1): 76-82, 92. DOI: 10.3969/j.issn.1000-5641.2019.01.009.

基金项目

国家自然科学基金（61275128，11474348）

第一作者

任宏红, 女, 硕士研究生, 研究方向为理论物理

通信作者

郭迎春, 女, 副教授, 研究方向为强场与物质相互作用.E-mail:ycguo@phy.ecnu.edu.cn

文章历史

收稿日期：2017-12-20

Contents Abstract Full text Figures/Tables PDF

整数阶复宗量贝塞尔函数的计算程序研究

任宏红 ¹, 郭迎春 ¹, 王兵兵 ²

1. 华东师范大学物理与材料科学学院, 上海 200241;
2. 中国科学院物理研究所凝聚态物理国家重点实验室光物理实验室, 北京 100190

收稿日期：2017-12-20

基金项目：国家自然科学基金（61275128，11474348）

第一作者：任宏红, 女, 硕士研究生, 研究方向为理论物理

通信作者：郭迎春, 女, 副教授, 研究方向为强场与物质相互作用.E-mail:ycguo@phy.ecnu.edu.cn

摘要：鉴于目前的算法程序集中没有现成的计算复宗量贝塞尔函数的程序，本文基于贝塞尔函数的逆向递推关系编写了计算整数阶复宗量第一类贝塞尔函数的Fortran程序源代码.与Matlab软件的计算结果比较，两者至少有12位有效数字一致.接着运用此程序，分析了徐士良的《FORTRAN常用算法程序集》中的纯虚宗量的贝塞尔函数，即变形贝塞尔函数程序的准确度，发现其准确度为6位有效数字.最后，对基于实宗量贝塞尔函数和纯虚宗量贝塞尔函数相乘然后用无限求和来计算复宗量贝塞尔函数值的方法的准确性进行了探讨.证明其仅能对有限的贝塞尔函数进行准确计算.这是由于当求和项中有远大于最终的求和项时，会导致求和结果的有效数字减少甚至完全错误.

关键词：复宗量贝塞尔函数 Fortran源程序逆递推计算

Program for calculating the integer order of Bessel functions with complex arguments

REN Hong-hong ¹, GUO Ying-chun ¹, WANG Bing-bing ²

1. School of Physics and Materials Science, East China Normal University, Shanghai 200241, China;
2. Laboratory of Optical Physics, Beijing National Laboratory of Condensed Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China

Abstract: Fortran source code for calculating the integer order of Bessel functions of the first kind with complex arguments is presented. The method is based on the backward recurrence relation of Bessel functions. Values of the Bessel function generated by our program are in-agreement with the values generated by Matlab to at least 12 significant digits. We use the program to calculate the integer order of Bessel functions of the first kind with pure imaginary arguments, provided by Xu Shiliang's Fortran algorithm assembly. The results show that the first 6 significant digits are accurate. We also analyze the algorithm for calculating Bessel functions with complex arguments, which use the infinite sum of the product of the real arguments of the Bessel function and the pure imaginary arguments of the Bessel function, provided by Xu Shiliang's algorithm assembly. The results show that this algorithm does not always get accurate values for Bessel functions with complex arguments. The reason lies with the fact that the term in the sum larger than the function value causes the loss of significant digits.

Keywords: Bessel function with complex arguments Fortran source code backward recurrence relation of Bessel function

0 引言

整数阶贝塞尔函数是许多辐射、散射和波导问题的数学解, 贝塞尔函数中复宗量以及虚宗量和材料损耗、漏波等相关.贝塞尔函数也应用于电磁场与物质相互作用的相关现象的数学描述, 如阈上电离及高次谐波等^[1-3].所以贝塞尔函数的精确计算尤为重要.

针对整数阶第一类贝塞尔函数, 在常用算法程序书中, 例如徐士良的《FORTRAN常用算法程序集》^[4], 以及《Numerical Recipes: The Art of Scientific Computing》^[5]等已经有实宗量的贝塞尔函数的程序和纯虚宗量的贝塞尔函数即变形的贝塞尔函数的程序, 但是还没有计算复宗量贝塞尔函数的程序. Du Toit^[6]探讨了复宗量的贝塞尔函数的计算方法, 提出了采用逆向贝塞尔函数递推关系进行计算, 文中给出了一些典型的数值结果, 指出这种方法对于大的宗量不能计算, 而对于大的宗量, 采用幂级数展开截断的方法给出了结果.魏彦玉等人^[7]采用级数展开方法对虚宗量的贝塞尔函数进行了计算, 对小的宗量计算结果很好.张爽等人^[8]针对这种方法的计算进行了误差分析, 指出了此方法不适合大宗量贝塞尔函数计算的原因.张善杰等人^[9]对任意实数阶复宗量贝塞尔函数的递推算法进行了研究, 并探讨了如何验证程序的正确性和分析了计算结果的精度.

本文针对第一类整数阶复宗量贝塞尔函数, 重新考察了基于逆递推公式对贝塞尔函数的计算, 编写了适用于任意大小的复宗量的第一类整数阶贝塞尔函数Fortran源代码.本程序计算的贝塞尔函数的结果与Matlab软件结果进行比较, 二者一致到第12位有效数字.由于常用算法程序书(如文献[4])中, 很容易找到实宗量的贝塞尔函数程序和纯虚宗量贝塞尔函数的程序, 基于二者相乘求和, 原则上可以计算复宗量的贝塞尔函数.本文探讨了这种方法计算贝塞尔函数的可行性.在这之前还探讨了算法程序集中变形贝塞尔函数的计算精度.

本文安排如下, 第1节简要阐述复宗量贝塞尔函数的理论基础, 并举例比较此程序的结果与Matlab的结果; 第2节运用给出的程序分析变形贝塞尔函数程序的准确度; 第3节分析基于实宗量贝塞尔函数和虚宗量贝塞尔函数相乘无穷求和的方法, 计算复宗量贝塞尔函数的可行性; 附录中给出源代码.

1 计算复宗量的贝塞尔函数的原理

贝塞尔函数满足递推关系:

$ \begin{align} {\rm J}_{n-1}(z)= \frac{2n}{z}{\rm J}_n(z)-{\rm J}_{n+1}(z), \end{align} $

(1)

令

$ \begin{align} {\rm J}_n(z)=\frac{B_n(z)}{S}, \end{align} $

(2)

从而有

$ \begin{align} B_{n-1}(z)= \frac{2n}{z}B_n(z)-B_{n+1}(z). \end{align} $

(3)

因为J$_n(z)$随$|z|$的增加而减小, 所以采用逆向的递推关系, 即令$B_q (z)=1$, $B_{q+1}(z)=0$, 忽略大于$q$阶的所有项, 反复运用上式计算可得到$B_n (z)$.

下面计算式(2)中的$S$.因为贝塞尔函数有下列性质:

$ \begin{align} 1={\rm J}_0(z)+2\sum\limits_{k=1}^{ \infty } {\rm J}_{2k}(z), \end{align} $

(4)

$ \cos z={\rm J}_0(z)+2\sum\limits_{k=1}^{ \infty } {\rm J}_{2k}(z)(-1)^{k}, $

(5)

所以$S$满足：

$ \begin{align} S=B_0(z)+2\sum\limits_{k=1}^{ \infty } B_{2k}(z) \approx B_0(z)+2\sum\limits_{k=1}^{q/2} B_{2k}(z), \end{align} $

(6)

$ S\cos z=B_0(z)+2\sum\limits_{k=1}^{ \infty } B_{2k}(z)(-1)^{k} \approx B_0(z)+2\sum\limits_{k=1}^{q/2} B_{2k}(z)(-1)^{k}. $

(7)

当${\rm Im} (z)\leq 1$时, 采用式(6)计算$S$, 当${\rm Im} (z) > 1$时, 采用式(7)计算$S$.这样做可以保障$S$中无穷求和的每一项都小于$S$值, 从而不会带来有效数字的减少.

文献[5]指出, 编程计算的误差有舍入误差(round off error)、截止误差(truncation error)和计算的稳定度误差(stability error).舍入误差是由计算机精度所决定的, 如双精度型的变量含16位有效数字, 误差是最后一位有效数字的量级.此数据精度所引入的误差为舍入误差.舍入误差一般会随计算次数的增加而增加, 一般会影响结果的最后两位有效数字.在某些特殊情况下会引入很大的误差, 从而减少有效数字.如两个几乎相等的同型数据相减.截止误差是由计算过程中将无穷求和截取为有限求和而引入的误差.稳定度误差是由于计算方法不当(称为不稳定误差)使最初阶段的舍入误差被连续放大而导致的误差, 这种误差甚至会淹没真实结果.减小这三种误差是我们程序设计的出发点, 文献[6]已证实我们所采用的这种逆递推方法是稳定的.不同情况下选取不同的$S$, 目的是减小舍入误差. $q$的选取决定了截止误差, 只要$q$值取得足够大, 截止误差就会变小, 直到小于舍入误差.具体表现为, $q$增大到某值后, 如果再继续增大, 计算的结果将保持不变.据此, 为计算准确又节省时间, 我们将$q$取值为

$ \begin{align} q=n+\sqrt{40n}+100+2|z|. \end{align} $

(8)

可见, 贝塞尔函数的阶数越大, 宗量的模越大, 需要的$q$值就越大.

计算结果表明, 我们的结果和Matlab软件的结果一致到至少12位有效数字, 另外, 我们计算了文献[4]表 1中的几个典型数据, 再现了文献的全部结果, 并给出了文献中没有计算出的结果, 列在表 1中.对于小的宗量, 小的阶数, 由公式(8)可见, $q$值会相对很小, 此时不仅$q$所决定的截止误差小于舍入误差, 由于$q$相对小, 计算的次数少, 和大阶数大宗量的情况相比, 舍入误差也不会大, 所以计算的结果应该更为精确.如表 1中的前两行所示, 我们计算的结果与Matlab的结果一致到14位有效数字.文献[4]强调这种基于逆向递推公式的方法不适合计算大的宗量, 认为计算如10000000+i333大宗量的贝塞尔函数是不可行的, 该文献提出运用幂级数展开截断的方法来计算大宗量贝塞尔函数, 而我们就是采用逆递推公式的方法, 将如10000000+i333这样大宗量的贝塞尔函数计算了出来.得到的结果与该文献中基于展开截断的方法的结果基本一致(见表 1).我们认为文献中采用这种逆向递推的方法没有计算出来的原因, 一方面是我们的逆递推起点$q$值比文献中的大, 还有可能是计算机的性能的局限造成的.

表 1 J$_n(z)$的计算举例 Tab.1 Numerical examples of J$_n(z)$

$n$	$z$	本文结果	[1]中基于逆向递推的结果
5	18+0i	$-$0.155370098779049	$-$0.155370098779049
5	0+18i	i3.05782771756610$×$10$^{6}$	i3.05782771756610$×$10$^{6}$
50	$-$180+100i	3.235775172706374$×$10$^{40}$$-$i2.054938087958093$×$10$^{40}$	3.23577517270637$×$10$^{40}$$-$i2.05493808795809$×$10$^{40}$
600	$-$4000$-$300i	$-$2.767691224061821$×$10$^{126}$+i3.172834828322395$×$10$^{126}$	$-$2.76769122406189$×$10$^{126}$+i3.172834828322395$×$10$^{126}$
1000	600$-$200i	5.373339700916123$×$10$^{-104}$+i1.658289459267489$×$10$^{-103}$	5.37333970091597$×$10$^{-104}$+i1.65828945926745$×$10$^{-103}$
0	10000000+333i	$-$1.810334213201858$×$10$^{140}$$-$i4.938482320281676$×$10$^{140}$	Not computable with this algorithm
4090	10000000+333i	4.878566679363567$×$10$^{140}$$-$i1.965697327877443$×$10$^{140}$	Not computable with this algorithm
3300	3000+10i	$-$3.031921315712407$×$10$^{-42}$$-$i2.477259163766169$×$10$^{-41}$	$-$3.03192131571261$×$10$^{-42}$$-$i2.47725916376618$×$10$^{-41}$

$n$	$z$	[1]中基于幂级数展开截断的结果	Matlab的结果
5	18+0i	$-$0.155370098779049	$-$0.155370098779049
5	0+18i	$-$3.38362801581971$×$10$^{-8}$+i3.05782771756610$×$10$^{6}$	1.872379463330692$×$10$^{-9}$+i3.057827717566103$×$10$^{6}$
50	$-$180+100i	3.23577517271025$×$10$^{40}$$-$i2.05493808795888$×$10$^{40}$	3.235775172706450$×$10$^{40}$$-$i2.054938087958026$×$10$^{40}$
600	$-$4000$-$300i	$-$2.76769122367143$×$10$^{126}$+i3.172834828322395$×$10$^{126}$	$-$2.767691224059726$×$10$^{126}$+i3.172834828323653$×$10$^{126}$
1000	600$-$200i	Not computable with this algorithm	5.373339700915873$×$10$^{-104}$ +i1.658289459267151$×$10$^{-103}$
0	10000000+333i	$-$1.81033421565001$×$10$^{140}$$-$i4.93848231938389$×$10$^{140}$	$-$1.810334213201893$×$10$^{140}$$-$i4.938482320281608$×$10$^{140}$
4090	10000000+333i	4.878566815077401$×$10$^{140}$$-$i1.96569732864925$×$10$^{140}$	4.878566679363453$×$10$^{140}$$-$i1.965697327877465$×$10$^{140}$
3300	3000+10i	Not computable with this algorithm	$-$3.031921315712416$×$10$^{-42}$$-$i2.477259163766451$×$10$^{-41}$

表 1 J$_n(z)$的计算举例 Tab.1 Numerical examples of J$_n(z)$

2 变形贝塞尔函数的准确度

文献[4]中给出了计算变形贝塞尔函数(即纯虚数的贝塞尔函数)的程序, 下面将采用第1节的程序来探讨此变形贝塞尔函数的准确度.

变形的贝塞尔函数I$_n(y)$与纯虚数的贝塞尔函数J$_n({\rm i}y)$是等价的, 它们的关系是

$ \begin{align} {\rm J}_n(iy)=i^n{\rm I}_n(y), \end{align} $

(9)

程序集中I$_n (y)$的计算和第1节一样, 也是采用逆向的贝塞尔函数递推公式, 即

$ \begin{align} {\rm I}_{n-1}(y)= \frac{2n}{y}{\rm I}_n(y)-{\rm I}_{n+1}(y). \end{align} $

(10)

令

$ \begin{align} {\rm I}_n(y)=\frac{B_n}{S}, \end{align} $

(11)

并令$B_q(y)=1$, $B_q(y)=0$, 忽略大于$q$阶的所有项, 反复运用上式计算得到$B_q(y)$, 与第2部分中我们的做法不同的是, $S$不是采用公式(6)和(7)求得, 而是运用拟合的方法首先计算出${\rm I}_0(y)$, 从而有

$ \begin{align} S=\frac{B_0 (y)}{{\rm I}_0(y)}. \end{align} $

(12)

所以${\rm I}_n (y)$的准确度不仅取决于$q$是否足够大, 还取决于${\rm I}_0 (y)$的准确度.所以, 我们采用第1节提供的程序对${\rm I}_0 (y)$进行了计算, 结果列于表 2.

表 2 ${\rm I}_0(y)$的计算举例 Tab.2 Numerical examples of ${\rm I}_0(y)$

可见程序集中${\rm I}_0(y)$的结果准确到第6位有效数字.进而确定程序集中的纯虚数的贝塞尔函数即变形的贝塞尔函数${\rm I}_n(y)$的准确度为6个有效数字, 远小于第1节提供的12位有效数字.

3 无穷求和计算贝塞尔函数

由于程序集[4]中已经提供了实宗量的贝塞尔函数的程序以及纯虚宗量的贝塞尔函数的程序, 很自然想到由下面的无穷求和(式(13))来计算复宗量的贝塞尔函数:

$ \begin{align} {\rm J}_n (z=x+{\rm{i}}y)=\sum\limits_{k= - \infty }^{ + \infty } {\rm J}_k(x){\rm I}_{n-k}(y)({\rm{i}})^{n-k}. \end{align} $

(13)

所以这一节, 我们将探讨基于这种无穷求和来计算复宗量贝塞尔函数的准确度.

通过与第1节中提供的程序计算结果的比较发现, 采用无穷求和来计算的结果仅仅在有限范围内是正确的.一般来讲, 对于确定的自变量, 阶数变大时, 计算结果将会不准确.我们定义准确的结果指的是计算结果的有效数字与第1节提供的程序的结果达到前6位一致. 表 3给出了几个特定自变量, 无穷求和方法能给出正确结果的贝塞尔函数的阶数范围, 如:对于${\rm J}_n (50+{\rm{i}}100)$, 当$-21 < n< 21$时, 无穷求和能够给出正确的结果准确到6位有效数字.

表 3 运用程序集[4]中的现有程序结合公式(13)能准确计算${\rm J}_{n}(x+{\rm i}y)$的$n$的范围 Tab.3 The range of the orders of ${\rm J}_n(x+{\rm i}y)$ which can be accurately

无穷求和的方法会给出不好的结果的原因在于式子(13)左边的贝塞尔函数值小于求和中的个别的项值, 从而减少了最终求和中准确的有效数字的个数.例如

$ \begin{align} {\rm J}_{35} (50+{\rm{i}}40)=\sum\limits_{k= - \infty }^{ + \infty } {\rm J}_k(50){\rm I}_{35-k}(40)({\rm{i}})^{35-k}. \end{align} $

(14)

${\rm J}_{35} (50+i40)$无穷求和方法给出的结果是: $-$5.14206$-$i2.85989, 相较于第1节给出的结果$-$5.14164$-$i3.22674, 在第一位有效数字就有误差, 原因在于求和项中, 很多项远远超过求和后的结果, 如: $k$为45时, ${\rm J}_k(50){\rm I}_{35-k}(40)(i)^{35-k}=-i8.554382770939035 \times 10^8$, 由第2节的讨论可知, 此数据前6位有效数字是准确的, 十位上的数据是不准确的, 所以求和后的结果, 十位以后都是不准确的.而$J_{35}(50+i40)$的结果的虚部最大在个位上, 所以造成了结果的完全错误.

4 结论

我们基于逆向的贝塞尔函数的递推关系, 给出了计算第一类复宗量整数阶贝塞尔函数的Fortran源代码, 探讨了现有的算法程序集^[4]中的纯虚数的贝塞尔函数的准确度为6个有效数字, 针对运用算法程序集中的实宗量的贝塞尔函数和纯虚宗量贝塞尔函数的程序, 采用它们的乘积对阶数的无穷求和来计算复宗量的贝塞尔函数的方法, 讨论了其计算的准确度, 通过与本文所给程序的结果比较, 得出其仅能对有限的贝塞尔函数进行计算.这是由于无穷求和项中有大于最终求和的项而导致结果的准确的有效数字减少甚至结果完全错误, 这也是物理量的求和计算中需要注意的问题.

附录

计算整数阶复宗量第一类贝塞尔函数的Fortran源代码

function xbessj(nw, xw)!nw:贝塞尔函数的阶数, xw:复宗量

implicit double precision (a-h, o-w)

implicit double complex(x-z)

parameter (iacc=40, bigno=1.d10, bigni=1.d-10)

n=abs(nw)

aw=cdabs(xw)

if(aw.lt.1.e-11)then

xbessj=0.

else

xtox=2./xw

m=2^*((n+int(sqrt(float(iacc^*n))))/2)+100+2^*int(aw)

xbessj=0.

jsum=0

xsum=0.

xbjp=0.

xbj=1.

do 12 j=m, 1, -1

xbjm=j^*xtox^*xbj-xbjp

xbjp=xbj

xbj=xbjm

if(cdabs(xbj).gt.bigno) then

xbj=xbj^*bigni

xbjp=xbjp^*bigni

xbessj=xbessj^*bigni

xsum=xsum^*bigni

endif

if((jsum.ne.0).and. (cdabs(dimag(xw)).le.1.0d0)) xsum=xsum+xbj

if((jsum.ne.0).and. (cdabs(dimag(xw)).gt.1.0d0))then

k=mod((j-1)/2, 2)

xsum=xsum+xbj^*(-1.0d0)^**k

endif

jsum=1-jsum

if(j.eq.n) xbessj=xbjp

12 continue

if(n.eq.0) xbessj=xbj

xsum=2.0d0^*xsum-xbj

if (cdabs(dimag(xw)).gt.1.0d0) then

xcos=(cdexp(dcmplx(0.0d0, 1.0d0)^*xw)+

^*cdexp(dcmplx(0.0d0, -1.0d0)^*xw))/2.0d0

xsum=xsum/xcos

endif

参考文献

[1]	LEWENSTEIN M, BALCOU P, IVANOV M Y, et al. Theory of high harmonic generation by low frequency laser fields[J]. Physical Review A, 1994, 49(3): 2117-2132. DOI:10.1103/PhysRevA.49.2117
[2]	GUO Y, FU P, YAN Z C, et al. Imaging the geometrical structure of the H₂⁺ molecular ion by high order above-threshold ionization in an intense laser field[J]. Physical Review A, 2009, 80(6): 3694-3697.
[3]	JIA X Y, LI W D, FAN J, et al. Suppression effect in the nonsequential double ionization of molecules by an intense laser field[J]. Physical review A, 2008, 77(6): 3195-3199.
[4]	徐士良. FORTRAN常用算法程序集[M]. 2版. 北京: 清华大学出版社, 2012.
[5]	PRESS W H, TEUKOLSKY S A, VETTERLING W T, et al. Numerical Recipes:The Art of Scientific Computing[M]. 3rd ed. Cambridge: Cambridge University Press, 2007.
[6]	DU TOIT C F. The numerical computation of Bessel functions of the first and second kind for integer orders and complex arguments[J]. IEEE Transactions on Antennas and Propagation, 1990, 38(9): 1341-1349. DOI:10.1109/8.56985
[7]	魏彦玉, 宫玉彬, 王文祥. 任意阶复宗量贝塞尔函数的数值计算[J]. 电子科技大学学报, 1998, 27(2): 171-176.
[8]	张爽, 郭欣, 宋立军. 利用贝塞尔函数的级数形式进行数值计算的误差分析[J]. 长春大学学报, 2004, 14(2): 57-59. DOI:10.3969/j.issn.1009-3907-B.2004.02.019
[9]	张善杰, 唐汉. 任意实数阶复宗量第一类和第二类bessel函数的精确计算[J]. 电子学报, 1996, 24(3): 77-81. DOI:10.3321/j.issn:0372-2112.1996.03.020