Fast and accurate computation of Bessel functions with (large) complex order and argument
Introduction
In many applications involving cylindrical configurations, products of Bessel functions arise [1], [2], [3]. Moreover, if the pertaining structure has a large radius (of curvature) and allows for radiation, the order of the Bessel functions becomes large and complex. The radiation loss is directly related to the imaginary part of the order. This implies that for small losses, the real part exceeds the imaginary part by several orders of magnitude [4], [5]. The argument of the Bessel functions, on the other hand, can be large and complex as well, e.g. in the case of a bent optical fibre [6]. In such a problem, cancellations of the various Bessel function products occur in the kernel of (multiple) integrals. Hence, computation has to be both fast and highly accurate.
For common ranges of order and argument, it is well-known how to evaluate the Bessel functions numerically, and many routines are available to perform the actual computation. However, the available (double precision) routines fail when it comes to computing a product of Bessel functions with large complex order and argument, especially when a high relative precision () is required [7], [8]. With radiation problems in mind, we have developed a method and a computer code in which high accuracy can be attained, even up to machine precision, for orders, , that vary over the entire range of .
Since asymptotic expansion techniques cannot be guaranteed to provide sufficient accuracy over the total range of [9], we take standard integral representations of the Bessel functions as our starting point [10]. Similar to the approach of Remenets [11] for Schläfli contour integrals, we deform the integration paths into steepest descent paths (SDPs). This deformation leads to an accurate and rapid computation of the Bessel functions, if appropriate scaling factors are introduced to avoid overflow and underflow in the numerical scheme [12]. However, for complex order and argument, analytical expressions of the SDPs do not exist, and hence, the location of the paths have to be determined numerically. This additional computational effort can be avoided through the use of well-chosen piecewise linear path segments (LPSs) for the integration at hand, which reduces computation times by as much as three orders of magnitude. As the integral representations also hold for these modified paths, the accuracy of the results is not influenced by this change.
We will first discuss our numerical scheme with its computational bottlenecks. Numerical results for both small and large complex orders and arguments are computed in Section 3. We have verified the results for the smaller arguments with Mathematica™ [8]. As reference packages for the verification of the results for larger parameters do not exist, we make use of Wronskian relationships [13], which involve products of Bessel functions. Further, we show a comparison of computation times between integration along the SDPs and along the LPSs, which demonstrates that the latter is faster by as much as three orders of magnitude.
Section snippets
Integral representation of the Bessel function
As the starting point of our numerical scheme, we consider the following integral representations for the Bessel functions [10], [13]:where , and . This region of validity can be extended via analytic continuation. Note that modified Bessel functions are covered as well, as the argument of the Bessel function may be complex, i.e. [13]
Numerical results
We have developed a numerical code to compute results for a product of (modified) Bessel functions. The modified Bessel functions are directly related to the Bessel functions by Eq. (2). We are not aware of a single mathematical package which can be used as a reference to validate our results over the entire range of complex orders and arguments that occur in practice. For validation purposes, we have computed the absolute value of the product and settled for a relatively
Conclusions
To compute Bessel functions with (large) complex order and argument over a vast range in a fast and highly accurate way, we have used integral representations. Instead of integration along the non-oscillatory steepest descent paths, we have shown that integration along well-chosen linear path segments reduces computation time considerably (three orders of magnitude), while retaining the same high relative accuracy of . For relatively small complex order and argument, the code has been
References (19)
- et al.
On an integral related to biaxially anisotropic media
J. Comput. Appl. Math.
(2002) Anti-plane wave motion in a cylindrically orthotropic elastic solid containing a stress-free crack
Wave Motion
(2004)- et al.
Coulomb and Bessel functions of complex arguments and order
J. Comput. Phys.
(1986) Computation of Hankel (Bessel) functions of complex index and argument by numerical integration of a Schläfli contour integral
USSR Comput. Math. Math. Phys.
(1973)- et al.
Evaluation of the modified Bessel function of the third kind of imaginary orders
J. Comput. Phys.
(2002) - et al.
Mathematical study of the two-dimensional lasing problem for whispering-gallery modes in a circular dielectric microcavity
Opt. Quant. Electron.
(2004) - et al.
A model impedance-angle formalism: Schemes for accurate graded-index bent-slab calculations and optical fiber mode counting
Radio Sci.
(2003) - et al.
Analytic approach to dielectric optical bent slab waveguides
Opt. Quant. Electron.
(2005) - et al.
Bending loss in optical fibers – a full-wave approach
J. Opt. Soc. Am. B
(2007)
Cited by (8)
Numerical calculation of Bessel, Hankel and Airy functions
2012, Computer Physics CommunicationsAutomatic particular solutions of arbitrary high-order splines associated with polyharmonic and poly-Helmholtz equations
2011, Engineering Analysis with Boundary ElementsQuasi-three-dimensional frequency-domain modeling to study size limitations of quantum-dot microdisk lasers
2010, Optics CommunicationsCitation Excerpt :We perform the computation of the Bessel functions using their integral representations [33] with well chosen integrations paths. The pertaining numerical scheme is described extensively in [34]. The solution of Eq. (11) is defined by the order of the (modified) Bessel function kφ, and consequently this leads to the solution of Eqs. (4) and (3) respectively.
Calculation of rail internal impedance based on tubular conductor model
2020, Zhongnan Daxue Xuebao (Ziran Kexue Ban)/Journal of Central South University (Science and Technology)High-accurate numerical computation of internal impedance of cylindrical conductors for complex arguments of arbitrary magnitude
2014, IEEE Transactions on Electromagnetic CompatibilityExact and asymptotic ser of distributed TAS/MRC in MIMO relay networks
2011, IEEE Transactions on Wireless Communications