A recursive algorithm for an efficient and accurate computation of incomplete Bessel functions

Incomplete Bessel functions, which are a computational challenge, were a subject of significant research. These functions appear when describing a plethora of phenomena in hydrology, statistics, and quantum mechanics [1, 2, 3, 4, 5, 6, 7, 8, 9]. Incomplete Bessel functions of zero order are also involved in a numerous applications to electromagnetic waves [10, 11, 12, 13, 14]. By introducing a recursive algorithm for the $G$ transformation for tail integrals and applying it to incomplete Bessel functions, the article [15] has received some attention recently [16, 17, 18, 19, 20, 21]. Most of this attention has been due to the algorithm for incomplete Bessel functions.

Integral representations of the incomplete Bessel functions are given by [22]:

The algorithm for the $G$ transformation [23, 24, 25, 26] applied to incomplete Bessel functions [15] appears as:

where:

where the $A_{r}^{i}$ are the coefficients in the Slevinsky-Safouhi formula I [27] and are given by:

with $(\mu,\nu,m,n)=(-2,-\nu-1,0,0)$.

The numerator $\tilde{\cal N}_{n}(x,y,\nu)$ in (2) is given by:

where the $A_{s}^{i}$ are the coefficients in the Slevinsky-Safouhi formula I with $(\mu,\nu,m,n)=(-2,\nu-1,0,0)$.

Therefore, the original algorithm involves sums over sums over sums, effectively making the complexity of the computation of the numerator an ${\cal O}(n^{3})$ process and the complexity of the computation of the denominator an ${\cal O}(n^{2})$ process. Therefore, calculating a sequence of approximations $\left\{G_{n}\right\}_{n\geq 0}$ results in an ${\cal O}(n^{4})$ algorithm.

In the following, we introduce an algorithm for the $G$ transformation that reduce the calculation to an easily programmable and parallel four-term recurrence relation: with one initialization, it computes the numerator; with another initialization, it computes the denominator. As there are only a finite number of arithmetic operations required in the recurrence relation, the resulting algorithm for calculating a sequence of approximations $\left\{G_{n}\right\}_{n\geq 0}$ is an ${\cal O}(n)$ algorithm.

The inductive proof of this recurrence relation follows and the notation is heavy because of the two variables $x$ and $y$ and the parameter $\nu$ already included in the incomplete Bessel function. As well, with this new algorithm comes the ability to program the $G$ transformation to unprecedentedly high order. We show some numerical results in the form of figures of the relative error for six different values. The reduction in complexity of the original algorithm comes with the benefit of a stable algorithm for high order.

To solve for the unknowns $\beta_{k,i}$, we must set up and solve a system of linear equations.

The $G_{n}^{(m)}$ transformation [23] truncates the asymptotic expansions (10) after $n$ terms each and the system is formed by differentiation. The approximation $G_{n}^{(m)}$ to $\int_{0}^{\infty}f(t)\mathrm{d}t$ is given as the solution of the system of $mn+1$ linear equations [23]:

where it is assumed that $\displaystyle\frac{\mathrm{d}^{l}}{\mathrm{d}x^{l}}G^{(m)}_{n}\equiv 0,\forall l>0$. In the above system (11), $\sigma_{k}=\mathrm{min}(s_{k},k+1)$ where $s_{k}$ is the largest of the integers $s$ such that $\displaystyle\lim_{x\to\infty}x^{s}f^{(k)}(x)=0$ holds, $k=0,1,\ldots,m-1$. Also, $G^{(m)}_{n}$ and $\bar{\beta}_{k,i}$ are the respective set of $mn+1$ unknowns.

The $G^{(1)}_{n}$ transformation can be written as the solution to the linear system (11) with $m=1$. Instead of solving the system of linear equations each time for each order $n$, it would be ideal to resolve each approximation $G^{(1)}_{n}$ in a recursive manner [15].

By considering the equation (11) for $l=0$:

and by isolating the summation on the right hand side, we obtain:

To eliminate the summation, and consequently all of the unknowns $\bar{\beta}_{0,i}$, we apply the $\displaystyle\left(x^{2}\frac{\mathrm{d}}{\mathrm{d}x}\right)$ operator $n$ times, obtaining:

which leads to a recursive algorithm for the $G^{(1)}_{n}$ transformation.

1. 1.

   Set:

   $${\cal N}_{0}(x)=\frac{F(x)}{x^{\sigma_{0}}f(x)}\qquad\textrm{and}\qquad{\cal D}_{0}(x)=\frac{1}{x^{\sigma_{0}}f(x)}.$$

   (15)

2. 2.

   For $n=1,2,\ldots,$ compute ${\cal N}_{n}(x)$ and ${\cal D}_{n}(x)$ recursively from:

   $${\cal N}_{n}(x)=\left(x^{2}\frac{\mathrm{d}}{\mathrm{d}x}\right){\cal N}_{n-1}(x)\qquad\textrm{and}\qquad{\cal D}_{n}(x)=\left(x^{2}\frac{\mathrm{d}}{\mathrm{d}x}\right){\cal D}_{n-1}(x).$$

   (16)

3. 3.

   For all $n$, the approximations $G^{(1)}_{n}(x)$ to $\displaystyle\left(\int_{0}^{x}+\int_{x}^{\infty}\right)f(t)\,\mathrm{d}t$ are given by:

   $$G^{(1)}_{n}(x)=\frac{{\cal N}_{n}(x)}{{\cal D}_{n}(x)}.$$

   (17)

For the incomplete Bessel functions, $\sigma_{0}=0$, and the algorithm for the $G_{n}^{(1)}$ transformation takes the form below. Let:

1. 1.

   Set:

   $${\cal N}_{0}(x,y,\nu,t)=\frac{F(t)}{f(t)}\qquad\textrm{and}\qquad{\cal D}_{0}(x,y,\nu,t)=\frac{1}{f(t)}.$$

   (19)

2. 2.

   For $n=1,2,\ldots,$ compute ${\cal N}_{n}(x,y,\nu,t)$ and ${\cal D}_{n}(x,y,\nu,t)$ recursively from:

   	$\displaystyle{\cal N}_{n}(x,y,\nu,t)$	$\displaystyle=$	$\displaystyle\left(t^{2}\frac{\mathrm{d}}{\mathrm{d}t}\right){\cal N}_{n-1}(x,y,\nu,t)$
   	$\displaystyle{\cal D}_{n}(x,y,\nu,t)$	$\displaystyle=$	$\displaystyle\left(t^{2}\frac{\mathrm{d}}{\mathrm{d}t}\right){\cal D}_{n-1}(x,y,\nu,t).$		(20)

3. 3.

   For all $n\geq 0$, the approximations $\tilde{G}^{(1)}_{n}(x,y,\nu)$ to $K_{\nu}(x,y)$ are given by:

   $$\tilde{G}^{(1)}_{n}(x,y,\nu)=\frac{\tilde{\cal N}_{n}(x,y,\nu,1)}{{\cal D}_{n}(x,y,\nu,1)},$$

   (21)

   where:

   $$\tilde{\cal N}_{n}(x,y,\nu,t)={\cal N}_{n}(x,y,\nu,t)-F(t){\cal D}_{n}(x,y,\nu,t).$$

   (22)

We improve an existing algorithm for incomplete Bessel functions based on the $G$ transformation [15] by developing a recurrence relation the numerator sequence and the denominator sequence whose ratio form the sequence of approximations. By finding this recurrence relation, we reduce the complexity from ${\cal O}(n^{4})$ to ${\cal O}(n)$. We plot relative error showing that the algorithm is capable of extremely high accuracy for incomplete Bessel functions. The stability and convergence appear to be remarkable.