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 [math]\displaystyle{ x }[/math].
The modified Bessel functions of the first kind [math]\displaystyle{ I_n(x) }[/math] satisfy the recurrence relation
- [math]\displaystyle{ I_{n-1}(x)=\frac{2n}{x}I_n(x)+I_{n+1}(x) }[/math].
However, the modified Bessel functions of the second kind [math]\displaystyle{ K_n(x) }[/math] also satisfy the same recurrence relation
- [math]\displaystyle{ K_{n-1}(x)=\frac{2n}{x}K_n(x)+K_{n+1}(x) }[/math].
The first solution decreases rapidly with [math]\displaystyle{ n }[/math]. The second solution increases rapidly with [math]\displaystyle{ n }[/math]. Miller's algorithm provides a numerically stable procedure to obtain the decreasing solution.
To compute the terms of a recurrence [math]\displaystyle{ a_0 }[/math] through [math]\displaystyle{ a_N }[/math] according to Miller's algorithm, one first chooses a value [math]\displaystyle{ M }[/math] much larger than [math]\displaystyle{ N }[/math] and computes a trial solution taking initial condition[math]\displaystyle{ a_M }[/math] to an arbitrary non-zero value (such as 1) and taking [math]\displaystyle{ a_{M+1} }[/math] and later terms to be zero. Then the recurrence relation is used to successively compute trial values for [math]\displaystyle{ a_{M-1} }[/math], [math]\displaystyle{ a_{M-2} }[/math] down to [math]\displaystyle{ a_0 }[/math]. 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:
- [math]\displaystyle{ I_0(x)+2\sum_{m=1}^\infty (-1)^mI_{2m}(x)=1 }[/math]
where the infinite summation becomes finite due to the approximation that [math]\displaystyle{ a_{M+1} }[/math] 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 [math]\displaystyle{ M }[/math] larger than the initial choice and confirming that the second set of results for [math]\displaystyle{ a_0 }[/math] through [math]\displaystyle{ a_N }[/math] agree within the first set within the desired tolerance. Note that to obtain this agreement, the value of [math]\displaystyle{ M }[/math] must be large enough such that the term [math]\displaystyle{ a_M }[/math] 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 [math]\displaystyle{ I_0(x) }[/math] and [math]\displaystyle{ I_1(x) }[/math] 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]
- [math]\displaystyle{ J_{n-1}(x)=\frac{2n}{x}J_{n}(x)-J_{n+1}(x) }[/math]
- [math]\displaystyle{ J_0(x)+2\sum_{m=1}^\infty J_{2m}(x)=1 }[/math].
The algorithm is particularly efficient in applications that require the values of the Bessel functions for all orders [math]\displaystyle{ 0 \cdots N }[/math] for each value of [math]\displaystyle{ x }[/math] compared to direct independent computations of [math]\displaystyle{ N+1 }[/math] separate functions.
References
- ↑ 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.0 2.1 Olver, F.W.J. (1964). "Error Analysis of Miller's Recurrence Algorithm". Math. Comp. 18 (85): 65–74. doi:10.2307/2003406.
- ↑ Németh, G. (1965). "Chebyshev Expansions for Fresnel Integrals". Numer. Math. 7 (4): 310–312. doi:10.1007/BF01436524.
- ↑ Hart, J.F. (1978). Computer Approximations (reprint ed.). Malabar, Florida: Robert E. Krieger. pp. 25–26. ISBN 978-0-88275-642-4.
- ↑ Gautschi, Walter (1967). "Computational aspects of three-term recurrence relations". SIAM Review 9: 24–82. doi:10.1137/1009002. https://www.cs.purdue.edu/homes/wxg/selected_works/section_09/029.pdf.
- ↑ Arfken, George (1985). Mathematical Methods for Physicists (3rd ed.). Academic Press. p. 576. ISBN 978-0-12-059820-5. https://archive.org/details/mathematicalmeth00arfk/page/576.
 |  |
