Fast and accurate computation of Bessel functions with (large) complex order and argument

https://doi.org/10.1016/j.amc.2008.10.057Get rights and content

Abstract

A fast and highly accurate computational scheme to compute a product of Bessel functions with (large) complex order and argument is of utmost importance in some problems in physics involving cylindrical configurations. Since these products can occur in the kernel of multiple integrals, and are subject to mutual cancellations, speed and accuracy of their numerical computation are key. Moreover, the order and argument can vary over a vast range of several orders of magnitude. To deliver a high relative accuracy (10-11) over the entire range, we use the integral representations of the Bessel functions, and introduce a deformation into steepest descent paths. The numerical computation is accelerated by three orders of magnitude, if these paths are replaced by well-chosen piecewise linear segments, while maintaining the set relative accuracy.

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 (10-11) 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 1<Re(ν)<106.

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]:Jν(x)=12πi--πi-+πie-xsinht+νtdt,Hν(1)(x)=1πi--πie-xsinht+νtdt;Hν(2)(x)=1πi-+πie-xsinht+νtdt,where {ν,x}C, and |arg(x)|<π/2. 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]Iν(ξ)=e-12νπiJνξe12πifor-π<arg(ξ)π/2,Kν(ξ)=

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 1<Re(ν)<106 and arguments that occur in practice. For validation purposes, we have computed the absolute value of the product Iν(ξ)Kν(ξ) 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 10-11. For relatively small complex order and argument, the code has been

References (19)

There are more references available in the full text version of this article.

Cited by (8)

  • Quasi-three-dimensional frequency-domain modeling to study size limitations of quantum-dot microdisk lasers

    2010, Optics Communications
    Citation 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)
  • Exact and asymptotic ser of distributed TAS/MRC in MIMO relay networks

    2011, IEEE Transactions on Wireless Communications
View all citing articles on Scopus
View full text