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

任宏红 郭迎春 王兵兵

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

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

中图分类号：0411.2 文献标志码：A DOI：10.3969/j.issn.1000-5641.2019.01.009

0引言

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

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

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

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

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

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

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

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

文献中给出了计算变形贝塞尔函数（即纯虚数的贝塞尔函数）的程序，下面将采用第l节的程序来探讨此变形贝塞尔函数的准确度.

可见程序集中Io（y）的结果准确到第6位有效数字.进而确定程序集中的纯虚数的贝塞尔函数即变形的贝塞尔函数In（y）的准确度为6个有效数字，远小于第1节提供的12位有效数字.

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

由于程序集中已经提供了实宗量的贝塞尔函数的程序以及纯虚宗量的贝塞尔函数的程序，很自然想到由下面的无穷求和（式（13））来计算复宗量的贝塞尔函数：所以这一节，我们将探讨基于这种无穷求和来计算复宗量贝塞尔函数的准确度.

通过与第1节中提供的程序计算结果的比较发现，采用无穷求和来计算的结果仅仅在有限范围内是正确的.一般来讲，对于确定的自变量，阶数变大时，计算结果将会不准确.我们定义准确的结果指的是计算结果的有效数字与第1节提供的程序的结果达到前6位一致.表3给出了几个特定自变量，无穷求和方法能给出正确结果的贝塞尔函数的阶数范围，如：对于Jn（50+i100），当一21

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

4结论

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