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.

References

1. Bickley, W.G.; Comrie, L.J.; Sadler, D.H.; Miller, J.C.P.; Thompson, A.J. (1952). British Association for the advancement of science, Mathematical Tables, vol. X, Bessel functions, part II, Functions of positive integer order. Cambridge University Press. ISBN 978-0521043212., cited in Olver (1964)
2. Olver, F.W.J. (1964). "Error Analysis of Miller's Recurrence Algorithm". Math. Comp. 18 (85): 65–74. doi:10.2307/2003406. JSTOR 2003406.
3. Németh, G. (1965). "Chebyshev Expansions for Fresnel Integrals". Numer. Math. 7 (4): 310–312. doi:10.1007/BF01436524.
4. Hart, J.F. (1978). Computer Approximations (reprint ed.). Malabar, Florida: Robert E. Krieger. pp. 25–26. ISBN 978-0-88275-642-4.
5. Gautschi, Walter (1967). "Computational aspects of three-term recurrence relations" (PDF). SIAM Review. 9: 24–82. doi:10.1137/1009002.
6. Arfken, George (1985). Mathematical Methods for Physicists (3rd ed.). Academic Press. p. 576. ISBN 978-0-12-059820-5.