Miller's recurrence algorithm

Last updated November 08, 2024

Miller's recurrence algorithm is a procedure for the backward calculation of a rapidly decreasing solution of a three-term 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 $x$ .

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

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

.

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

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

.

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

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

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 $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 $M$ larger than the initial choice and confirming that the second set of results for $a_{0}$ through $a_{N}$ agree within the first set within the desired tolerance. Note that to obtain this agreement, the value of $M$ must be large enough such that the term $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 $I_{0}(x)$ and $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]

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

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 $0\cdots N$ for each value of $x$ compared to direct independent computations of $N+1$ separate functions.

Related Research Articles

Bessel functions, first defined by the mathematician Daniel Bernoulli and then generalized by Friedrich Bessel, are canonical solutions $y (x)$ of Bessel's differential equation $for an arbitrary complex number, which represents the order of the Bessel function. Although and produce the same differential equation, it is conventional to define different Bessel functions for these two values in such a way that the Bessel functions are mostly smooth functions of .$

In mathematics, the gamma function is the most common extension of the factorial function to complex numbers. Derived by Daniel Bernoulli, the gamma function $is defined for all complex numbers except non-positive integers, and for every positive integer, The gamma function can be defined via a convergent improper integral for complex numbers with positive real part:$

In numerical analysis, an $n$ -point Gaussian quadrature rule, named after Carl Friedrich Gauss, is a quadrature rule constructed to yield an exact result for polynomials of degree $2 n - 1$ or less by a suitable choice of the nodes $x i$ and weights $w i$ for $i = 1, ..., n$ .

In mathematics, a recurrence relation is an equation according to which the $th term of a sequence of numbers is equal to some combination of the previous terms. Often, only previous terms of the sequence appear in the equation, for a parameter that is independent of; this number is called the order of the relation. If the values of the first numbers in the sequence have been given, the rest of the sequence can be calculated by repeatedly applying the equation.$

In mathematics, tables of trigonometric functions are useful in a number of areas. Before the existence of pocket calculators, trigonometric tables were essential for navigation, science and engineering. The calculation of mathematical tables was an important area of study, which led to the development of the first mechanical computing devices.

The Chebyshev polynomials are two sequences of polynomials related to the cosine and sine functions, notated as $and . They can be defined in several equivalent ways, one of which starts with trigonometric functions:$

In mathematics, a linear differential equation is a differential equation that is defined by a linear polynomial in the unknown function and its derivatives, that is an equation of the form $where a 0 (x), ..., a n (x) and b (x) are arbitrary differentiable functions that do not need to be linear, and y', ..., y (n) are the successive derivatives of an unknown function y of the variable x .$

In the physical sciences, the Airy function (or Airy function of the first kind) $Ai(x)$ is a special function named after the British astronomer George Biddell Airy (1801–1892). The function Ai(x) and the related function Bi(x), are linearly independent solutions to the differential equation $known as the Airy equation or the Stokes equation .$

In mathematics, Mathieu functions, sometimes called angular Mathieu functions, are solutions of Mathieu's differential equation

A fully polynomial-time approximation scheme (FPTAS) is an algorithm for finding approximate solutions to function problems, especially optimization problems. An FPTAS takes as input an instance of the problem and a parameter ε > 0. It returns as output a value which is at least $times the correct value, and at most times the correct value.$

In queueing theory, a discipline within the mathematical theory of probability, Buzen's algorithm (or convolution algorithm) is an algorithm for calculating the normalization constant G(N) in the Gordon–Newell theorem. This method was first proposed by Jeffrey P. Buzen in his 1971 PhD dissertation and subsequently published in a refereed journal in 1973. Computing G(N) is required to compute the stationary probability distribution of a closed queueing network.

In mathematics, Neville's algorithm is an algorithm used for polynomial interpolation that was derived by the mathematician Eric Harold Neville in 1934. Given n + 1 points, there is a unique polynomial of degree ≤ n which goes through the given points. Neville's algorithm evaluates this polynomial.

In mathematics, and more specifically in analysis, a holonomic function is a smooth function of several variables that is a solution of a system of linear homogeneous differential equations with polynomial coefficients and satisfies a suitable dimension condition in terms of D-modules theory. More precisely, a holonomic function is an element of a holonomic module of smooth functions. Holonomic functions can also be described as differentiably finite functions, also known as D-finite functions. When a power series in the variables is the Taylor expansion of a holonomic function, the sequence of its coefficients, in one or several indices, is also called holonomic. Holonomic sequences are also called P-recursive sequences: they are defined recursively by multivariate recurrences satisfied by the whole sequence and by suitable specializations of it. The situation simplifies in the univariate case: any univariate sequence that satisfies a linear homogeneous recurrence relation with polynomial coefficients, or equivalently a linear homogeneous difference equation with polynomial coefficients, is holonomic.

In mathematical analysis, the Bessel–Clifford function, named after Friedrich Bessel and William Kingdon Clifford, is an entire function of two complex variables that can be used to provide an alternative development of the theory of Bessel functions. If

In numerical analysis, a quasi-Newton method is an iterative numerical method used either to find zeroes or to find local maxima and minima of functions via an iterative recurrence formula much like the one for Newton's method, except using approximations of the derivatives of the functions in place of exact derivatives. Newton's method requires the Jacobian matrix of all partial derivatives of a multivariate function when used to search for zeros or the Hessian matrix when used for finding extrema. Quasi-Newton methods, on the other hand, can be used when the Jacobian matrices or Hessian matrices are unavailable or are impractical to compute at every iteration.

In numerical analysis, Gauss–Legendre quadrature is a form of Gaussian quadrature for approximating the definite integral of a function. For integrating over the interval $[-1, 1]$ , the rule takes the form:

In mathematics, the FEE method, or fast E-function evaluation method, is the method of fast summation of series of a special form. It was constructed in 1990 by Ekaterina Karatsuba and is so-named because it makes fast computations of the Siegel $E$ -functions possible, in particular of $.$

In computational geometry, a maximum disjoint set (MDS) is a largest set of non-overlapping geometric shapes selected from a given set of candidate shapes.

Quantum counting algorithm is a quantum algorithm for efficiently counting the number of solutions for a given search problem. The algorithm is based on the quantum phase estimation algorithm and on Grover's search algorithm.

In mathematics, and especially in numerical analysis, a homogeneous linear three-term recurrence relation is a recurrence relation of the form

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)
1 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.
↑ 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" (PDF). SIAM Review. 9: 24–82. doi:10.1137/1009002.
↑ Arfken, George (1985). Mathematical Methods for Physicists (3rd ed.). Academic Press. p. 576. ISBN 978-0-12-059820-5.

This page is based on this Wikipedia article
Text is available under the CC BY-SA 4.0 license; additional terms may apply.
Images, videos and audio are available under their respective licenses.

[Miller-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)

[Olver-2] 1 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.

[Nemeth-3] Németh, G. (1965). "Chebyshev Expansions for Fresnel Integrals". Numer. Math. 7 (4): 310–312. doi:10.1007/BF01436524.

[Hart-4] Hart, J.F. (1978). Computer Approximations (reprint ed.). Malabar, Florida: Robert E. Krieger. pp. 25–26. ISBN 978-0-88275-642-4.

[Gautschi-5] Gautschi, Walter (1967). "Computational aspects of three-term recurrence relations" (PDF). SIAM Review. 9: 24–82. doi:10.1137/1009002.

[Arkfen-6] Arfken, George (1985). Mathematical Methods for Physicists (3rd ed.). Academic Press. p. 576. ISBN 978-0-12-059820-5.

[1]

[2]

[3]

[4]

[5]

[6]