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 .

The modified Bessel functions of the first kind satisfy the recurrence relation

.

However, the modified Bessel functions of the second kind also satisfy the same recurrence relation

.

The first solution decreases rapidly with . The second solution increases rapidly with . Miller's algorithm provides a numerically stable procedure to obtain the decreasing solution.

To compute the terms of a recurrence through according to Miller's algorithm, one first chooses a value much larger than and computes a trial solution taking initial condition to an arbitrary non-zero value (such as 1) and taking and later terms to be zero. Then the recurrence relation is used to successively compute trial values for , down to . 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:

where the infinite summation becomes finite due to the approximation that 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 larger than the initial choice and confirming that the second set of results for through agree within the first set within the desired tolerance. Note that to obtain this agreement, the value of must be large enough such that the term 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 and 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]

.

The algorithm is particularly efficient in applications that require the values of the Bessel functions for all orders for each value of compared to direct independent computations of 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.
This article is issued from Wikipedia. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.