Miller's recurrence algorithm

Miller's recurrence algorithm is a procedure for calculating a rapidly decreasing solution of a linear recurrence relation developed by J. C. P. Miller.[1] It was originally developed to compute tables of the modified Bessel function[2] but also applies to Bessel functions of the first kind and has other applications such as computation of the coefficients of Chebyshev expansions of other special functions.[3]

Many families of special functions satisfy a recurrence relation that relates the values of the functions of different orders with common argument ${\displaystyle x}$.

The modified Bessel functions of the first kind ${\displaystyle I_{n}(x)}$ satisfy the recurrence relation

${\displaystyle I_{n-1}(x)={\frac {2n}{x}}I_{n}(x)+I_{n+1}(x)}$.

However, the modified Bessel functions of the second kind ${\displaystyle K_{n}(x)}$ also satisfy the same recurrence relation

${\displaystyle K_{n-1}(x)={\frac {2n}{x}}K_{n}(x)+K_{n+1}(x)}$.

The first solution decreases rapidly with ${\displaystyle n}$. The second solution increases rapidly with ${\displaystyle n}$. Miller's algorithm provides a numerically stable procedure to obtain the decreasing solution.

To compute the terms of a recurrence ${\displaystyle a_{0}}$ through ${\displaystyle a_{N}}$ according to Miller's algorithm, one first chooses a value ${\displaystyle M}$ much larger than ${\displaystyle N}$ and computes a trial solution taking initial condition${\displaystyle a_{M}}$ to an arbitrary non-zero value (such as 1) and taking ${\displaystyle a_{M+1}}$ and later terms to be zero. Then the recurrence relation is used to successively compute trial values for ${\displaystyle a_{M-1}}$, ${\displaystyle a_{M-2}}$ down to ${\displaystyle a_{0}}$. Noting that a second sequence obtained from the trial sequence by multiplication by a constant normalizing factor will still satisfy the same recurrence relation, one can then apply a separate normalizing relationship to determine the normalizing factor that yields the actual solution.

In the example of the modified Bessel functions, a suitable normalizing relation is a summation involving the even terms of the recurrence:

${\displaystyle I_{0}(x)+2\sum _{m=1}^{\infty }(-1)^{m}I_{2m}(x)=1}$

where the infinite summation becomes finite due to the approximation that ${\displaystyle a_{M+1}}$ and later terms are zero.

Finally, it is confirmed that the approximation error of the procedure is acceptable by repeating the procedure with a second choice of ${\displaystyle M}$ larger than the initial choice and confirming that the second set of results for ${\displaystyle a_{0}}$ through ${\displaystyle a_{N}}$ agree within the first set within the desired tolerance. Note that to obtain this agreement, the value of ${\displaystyle M}$ must be large enough such that the term ${\displaystyle a_{M}}$ is small compared to the desired tolerance.

In contrast to Miller's algorithm, attempts to apply the recurrence relation in the forward direction starting from known values of ${\displaystyle I_{0}(x)}$ and ${\displaystyle I_{1}(x)}$ obtained by other methods will fail as rounding errors introduce components of the rapidly increasing solution.[4]

Olver[2] and Gautschi[5] analyses the error propagation of the algorithm in detail.

For Bessel functions of the first kind, the equivalent recurrence relation and normalizing relationship are:[6]

${\displaystyle J_{n-1}(x)={\frac {2n}{x}}J_{n}(x)-J_{n+1}(x)}$

${\displaystyle J_{0}(x)+2\sum _{m=1}^{\infty }J_{2m}(x)=1}$.

The algorithm is particularly efficient in applications that require the values of the Bessel functions for all orders ${\displaystyle 0\cdots N}$ for each value of ${\displaystyle x}$ compared to direct independent computations of ${\displaystyle N+1}$ separate functions.