- §3.6(i) Introduction
- §3.6(ii) Homogeneous Equations
- §3.6(iii) Miller’s Algorithm
- §3.6(iv) Inhomogeneous Equations
- §3.6(v) Olver’s Algorithm
- §3.6(vi) Examples
- §3.6(vii) Linear Difference Equations of Other Orders

Many special functions satisfy second-order recurrence relations, or difference equations, of the form

3.6.1 | $${a}_{n}{w}_{n+1}-{b}_{n}{w}_{n}+{c}_{n}{w}_{n-1}={d}_{n},$$ | ||

or equivalently,

3.6.2 | $${a}_{n}{\mathrm{\Delta}}^{2}{w}_{n-1}+(2{a}_{n}-{b}_{n})\mathrm{\Delta}{w}_{n-1}+({a}_{n}-{b}_{n}+{c}_{n}){w}_{n-1}={d}_{n},$$ | ||

where $\mathrm{\Delta}{w}_{n-1}={w}_{n}-{w}_{n-1}$,
${\mathrm{\Delta}}^{2}{w}_{n-1}=\mathrm{\Delta}{w}_{n}-\mathrm{\Delta}{w}_{n-1}$, and $n\in \mathbb{Z}$. If
${d}_{n}=0$, $\forall n$, then the difference equation is *homogeneous*;
otherwise it is *inhomogeneous*.

Given numerical values of ${w}_{0}$ and ${w}_{1}$, the solution ${w}_{n}$ of the equation

3.6.3 | $${a}_{n}{w}_{n+1}-{b}_{n}{w}_{n}+{c}_{n}{w}_{n-1}=0,$$ | ||

with ${a}_{n}\ne 0$, $\forall n$, can be computed recursively for
$n=2,3,\mathrm{\dots}$. Unless exact arithmetic is being used, however, each step of
the calculation introduces rounding errors. These errors have the effect of
perturbing the solution by unwanted small multiples of ${w}_{n}$ and of an
independent solution ${g}_{n}$, say. This is of little consequence if the wanted
solution is growing in magnitude at least as fast as any other solution of
(3.6.3), and the recursion process is *stable*.

But suppose that ${w}_{n}$ is a nontrivial solution such that

3.6.4 | $${w}_{n}/{g}_{n}\to 0,$$ | ||

$n\to \mathrm{\infty}$. | |||

Then ${w}_{n}$ is said to be a *recessive* (equivalently, *minimal* or
*distinguished*) *solution* as $n\to \mathrm{\infty}$, and it is unique
except for a constant factor. In this situation the unwanted multiples of ${g}_{n}$
grow more rapidly than the wanted solution, and the computations are
*unstable*. Stability can be restored, however, by *backward
recursion*, provided that ${c}_{n}\ne 0$, $\forall n$: starting from ${w}_{N}$ and
${w}_{N+1}$, with $N$ large, equation (3.6.3) is applied to generate
in succession ${w}_{N-1},{w}_{N-2},\mathrm{\dots},{w}_{0}$. The unwanted multiples of ${g}_{n}$
now decay in comparison with ${w}_{n}$, hence are of little consequence.

The values of ${w}_{N}$ and ${w}_{N+1}$ needed to begin the backward recursion may be available, for example, from asymptotic expansions (§2.9). However, there are alternative procedures that do not require ${w}_{N}$ and ${w}_{N+1}$ to be known in advance. These are described in §§ 3.6(iii) and 3.6(v).

Because the recessive solution of a homogeneous equation is the fastest growing solution in the backward direction, it occurred to J.C.P. Miller (Bickley et al. (1952, pp. xvi–xvii)) that arbitrary “trial values” can be assigned to ${w}_{N}$ and ${w}_{N+1}$, for example, $1$ and $0$. A “trial solution” is then computed by backward recursion, in the course of which the original components of the unwanted solution ${g}_{n}$ die away. It therefore remains to apply a normalizing factor $\mathrm{\Lambda}$. The process is then repeated with a higher value of $N$, and the normalized solutions compared. If agreement is not within a prescribed tolerance the cycle is continued.

The normalizing factor $\mathrm{\Lambda}$ can be the true value of ${w}_{0}$ divided by its trial value, or $\mathrm{\Lambda}$ can be chosen to satisfy a known property of the wanted solution of the form

3.6.5 | $$\sum _{n=0}^{\mathrm{\infty}}{\lambda}_{n}{w}_{n}=1,$$ | ||

where the $\lambda $’s are constants. The latter method is usually superior when the true value of ${w}_{0}$ is zero or pathologically small.

For further information on Miller’s algorithm, including examples, convergence proofs, and error analyses, see Wimp (1984, Chapter 4), Gautschi (1967, 1997b), and Olver (1964a). See also Gautschi (1967) and Gil et al. (2007a, Chapter 4) for the computation of recessive solutions via continued fractions.

Similar principles apply to equation (3.6.1) when ${a}_{n}{c}_{n}\ne 0$, $\forall n$, and ${d}_{n}\ne 0$ for some, or all, values of $n$. If, as $n\to \mathrm{\infty}$, the wanted solution ${w}_{n}$ grows (decays) in magnitude at least as fast as any solution of the corresponding homogeneous equation, then forward (backward) recursion is stable.

A new problem arises, however, if, as $n\to \mathrm{\infty}$, the asymptotic behavior of ${w}_{n}$ is intermediate to those of two independent solutions ${f}_{n}$ and ${g}_{n}$ of the corresponding inhomogeneous equation (the complementary functions). More precisely, assume that ${f}_{0}\ne 0$, ${g}_{n}\ne 0$ for all sufficiently large $n$, and as $n\to \mathrm{\infty}$

3.6.6 | ${f}_{n}/{g}_{n}$ | $\to 0,$ | ||

${w}_{n}/{g}_{n}$ | $\to 0.$ | |||

Then computation of ${w}_{n}$ by forward recursion is unstable. If it also happens
that ${f}_{n}/{w}_{n}\to 0$ as $n\to \mathrm{\infty}$, then computation of ${w}_{n}$ by backward
recursion is unstable as well. However, ${w}_{n}$ can be computed successfully in
these circumstances by *boundary-value methods*, as follows.

Let us assume the normalizing condition is of the form ${w}_{0}=\lambda $, where $\lambda $ is a constant, and then solve the following tridiagonal system of algebraic equations for the unknowns ${w}_{1}^{(N)},{w}_{2}^{(N)},\mathrm{\dots},{w}_{N-1}^{(N)}$; see §3.2(ii). Here $N$ is an arbitrary positive integer.

3.6.7 | $$\left[\begin{array}{ccccc}\hfill -{b}_{1}\hfill & \hfill {a}_{1}\hfill & \hfill \hfill & \hfill \hfill & \hfill 0\hfill \\ \hfill {c}_{2}\hfill & \hfill -{b}_{2}\hfill & \hfill {a}_{2}\hfill \\ \hfill \hfill & \hfill \mathrm{\ddots}\hfill & \hfill \mathrm{\ddots}\hfill & \hfill \mathrm{\ddots}\hfill \\ \hfill \hfill & \hfill \hfill & \hfill {c}_{N-2}\hfill & \hfill -{b}_{N-2}\hfill & \hfill {a}_{N-2}\hfill \\ \hfill 0\hfill & \hfill \hfill & \hfill \hfill & \hfill {c}_{N-1}\hfill & \hfill -{b}_{N-1}\hfill \end{array}\right]\left[\begin{array}{c}\hfill {w}_{1}^{(N)}\hfill \\ \hfill {w}_{2}^{(N)}\hfill \\ \hfill \mathrm{\vdots}\hfill \\ \hfill {w}_{N-2}^{(N)}\hfill \\ \hfill {w}_{N-1}^{(N)}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {d}_{1}-{c}_{1}\lambda \hfill \\ \hfill {d}_{2}\hfill \\ \hfill \mathrm{\vdots}\hfill \\ \hfill {d}_{N-2}\hfill \\ \hfill {d}_{N-1}\hfill \end{array}\right].$$ | ||

Then as $N\to \mathrm{\infty}$ with $n$ fixed, ${w}_{n}^{(N)}\to {w}_{n}$.

To apply the method just described a succession of values can be prescribed for the arbitrary integer $N$ and the results compared. However, a more powerful procedure combines the solution of the algebraic equations with the determination of the optimum value of $N$. It is applicable equally to the computation of the recessive solution of the homogeneous equation (3.6.3) or the computation of any solution ${w}_{n}$ of the inhomogeneous equation (3.6.1) for which the conditions of §3.6(iv) are satisfied.

Suppose again that ${f}_{0}\ne 0$, ${w}_{0}$ is given, and we wish to calculate ${w}_{1},{w}_{2},\mathrm{\dots},{w}_{M}$ to a prescribed relative accuracy $\u03f5$ for a given value of $M$. We first compute, by forward recurrence, the solution ${p}_{n}$ of the homogeneous equation (3.6.3) with initial values ${p}_{0}=0$, ${p}_{1}=1$. At the same time we construct a sequence ${e}_{n}$, $n=0,1,\mathrm{\dots}$, defined by

3.6.8 | $${a}_{n}{e}_{n}={c}_{n}{e}_{n-1}-{d}_{n}{p}_{n},$$ | ||

beginning with ${e}_{0}={w}_{0}$. (This part of the process is equivalent to forward elimination.) The computation is continued until a value $N$ ($\ge M$) is reached for which

3.6.9 | $$\left|\frac{{e}_{N}}{{p}_{N}{p}_{N+1}}\right|\le \u03f5\underset{1\le n\le M}{\mathrm{min}}\left|\frac{{e}_{n}}{{p}_{n}{p}_{n+1}}\right|.$$ | ||

Then ${w}_{n}$ is generated by backward recursion from

3.6.10 | $${p}_{n+1}{w}_{n}={p}_{n}{w}_{n+1}+{e}_{n},$$ | ||

starting with ${w}_{N}=0$. (This part of the process is back substitution.)

The difference equation

3.6.11 | $${w}_{n+1}-2n{w}_{n}+{w}_{n-1}=0,$$ | ||

$n=1,2,\mathrm{\dots}$, | |||

is satisfied by ${J}_{n}\left(1\right)$ and ${Y}_{n}\left(1\right)$, where ${J}_{n}\left(x\right)$ and ${Y}_{n}\left(x\right)$ are the Bessel functions of the first kind. For large $n$,

3.6.12 | $${J}_{n}\left(1\right)\sim \frac{1}{{(2\pi n)}^{1/2}}{\left(\frac{\mathrm{e}}{2n}\right)}^{n},$$ | ||

3.6.13 | $${Y}_{n}\left(1\right)\sim {\left(\frac{2}{\pi n}\right)}^{1/2}{\left(\frac{2n}{\mathrm{e}}\right)}^{n},$$ | ||

(§10.19(i)). Thus ${Y}_{n}\left(1\right)$ is dominant and can be computed by forward recursion, whereas ${J}_{n}\left(1\right)$ is recessive and has to be computed by backward recursion. The backward recursion can be carried out using independently computed values of ${J}_{N}\left(1\right)$ and ${J}_{N+1}\left(1\right)$ or by use of Miller’s algorithm (§3.6(iii)) or Olver’s algorithm (§3.6(v)).

The Weber function ${\mathbf{E}}_{n}\left(1\right)$ satisfies

3.6.14 | $${w}_{n+1}-2n{w}_{n}+{w}_{n-1}=-(2/\pi )(1-{(-1)}^{n}),$$ | ||

for $n=1,2,\mathrm{\dots}$, and as $n\to \mathrm{\infty}$

3.6.15 | ${\mathbf{E}}_{2n}\left(1\right)$ | $\sim {\displaystyle \frac{2}{(4{n}^{2}-1)\pi}},$ | ||

3.6.16 | ${\mathbf{E}}_{2n+1}\left(1\right)$ | $\sim {\displaystyle \frac{2}{(2n+1)\pi}};$ | ||

see §11.11(ii). Thus the asymptotic behavior of the particular solution ${\mathbf{E}}_{n}\left(1\right)$ is intermediate to those of the complementary functions ${J}_{n}\left(1\right)$ and ${Y}_{n}\left(1\right)$; moreover, the conditions for Olver’s algorithm are satisfied. We apply the algorithm to compute ${\mathbf{E}}_{n}\left(1\right)$ to 8S for the range $n=1,2,\mathrm{\dots},10$, beginning with the value ${\mathbf{E}}_{0}\left(1\right)=-\mathrm{0.56865\hspace{0.33em}\hspace{0.17em}663}$ obtained from the Maclaurin series expansion (§11.10(iii)).

In the notation of §3.6(v) we have $M=10$ and $\u03f5=\frac{1}{2}\times {10}^{-8}$. The least value of $N$ that satisfies (3.6.9) is found to be 16. The results of the computations are displayed in Table 3.6.1. The values of ${w}_{n}$ for $n=1,2,\mathrm{\dots},10$ are the wanted values of ${\mathbf{E}}_{n}\left(1\right)$. (It should be observed that for $n>10$, however, the ${w}_{n}$ are progressively poorer approximations to ${\mathbf{E}}_{n}\left(1\right)$: the underlined digits are in error.)

$n$ | ${p}_{n}$ | ${e}_{n}$ | ${e}_{n}/({p}_{n}{p}_{n+1})$ | ${w}_{n}$ | ||||
---|---|---|---|---|---|---|---|---|

0 | $\mathrm{0.00000\hspace{0.33em}000}$ | $-\mathrm{0.56865\hspace{0.33em}663}$ | $-\mathrm{0.56865\hspace{0.33em}663}$ | |||||

1 | $\mathrm{0.10000\hspace{0.33em}000}$ | $\mathrm{\times 10\xb9}$ | $\mathrm{0.70458\hspace{0.33em}291}$ | $\mathrm{0.35229\hspace{0.33em}146}$ | $\mathrm{0.43816\hspace{0.33em}243}$ | |||

2 | $\mathrm{0.20000\hspace{0.33em}000}$ | $\mathrm{\times 10\xb9}$ | $\mathrm{0.70458\hspace{0.33em}291}$ | $\mathrm{0.50327\hspace{0.33em}351}$ | $\mathrm{\times 10\u207b\xb9}$ | $\mathrm{0.17174\hspace{0.33em}195}$ | ||

3 | $\mathrm{0.70000\hspace{0.33em}000}$ | $\mathrm{\times 10\xb9}$ | $\mathrm{0.96172\hspace{0.33em}597}$ | $\mathrm{\times 10\xb9}$ | $\mathrm{0.34347\hspace{0.33em}356}$ | $\mathrm{\times 10\u207b\xb9}$ | $\mathrm{0.24880\hspace{0.33em}538}$ | |

4 | $\mathrm{0.40000\hspace{0.33em}000}$ | $\mathrm{\times 10\xb2}$ | $\mathrm{0.96172\hspace{0.33em}597}$ | $\mathrm{\times 10\xb9}$ | $\mathrm{0.76815\hspace{0.33em}174}$ | $\mathrm{\times 10\u207b\xb3}$ | $\mathrm{0.47850\hspace{0.33em}795}$ | $\mathrm{\times 10\u207b\xb9}$ |

5 | $\mathrm{0.31300\hspace{0.33em}000}$ | $\mathrm{\times 10\xb3}$ | $\mathrm{0.40814\hspace{0.33em}124}$ | $\mathrm{\times 10\xb3}$ | $\mathrm{0.42199\hspace{0.33em}534}$ | $\mathrm{\times 10\u207b\xb3}$ | $\mathrm{0.13400\hspace{0.33em}098}$ | |

6 | $\mathrm{0.30900\hspace{0.33em}000}$ | $\mathrm{\times 10\u2074}$ | $\mathrm{0.40814\hspace{0.33em}124}$ | $\mathrm{\times 10\xb3}$ | $\mathrm{0.35924\hspace{0.33em}754}$ | $\mathrm{\times 10\u207b\u2075}$ | $\mathrm{0.18919\hspace{0.33em}443}$ | $\mathrm{\times 10\u207b\xb9}$ |

7 | $\mathrm{0.36767\hspace{0.33em}000}$ | $\mathrm{\times 10\u2075}$ | $\mathrm{0.47221\hspace{0.33em}340}$ | $\mathrm{\times 10\u2075}$ | $\mathrm{0.25102\hspace{0.33em}029}$ | $\mathrm{\times 10\u207b\u2075}$ | $\mathrm{0.93032\hspace{0.33em}343}$ | $\mathrm{\times 10\u207b\xb9}$ |

8 | $\mathrm{0.51164\hspace{0.33em}800}$ | $\mathrm{\times 10\u2076}$ | $\mathrm{0.47221\hspace{0.33em}340}$ | $\mathrm{\times 10\u2075}$ | $\mathrm{0.11324\hspace{0.33em}804}$ | $\mathrm{\times 10\u207b\u2077}$ | $\mathrm{0.10293\hspace{0.33em}811}$ | $\mathrm{\times 10\u207b\xb9}$ |

9 | $\mathrm{0.81496\hspace{0.33em}010}$ | $\mathrm{\times 10\u2077}$ | $\mathrm{0.10423\hspace{0.33em}616}$ | $\mathrm{\times 10\u2078}$ | $\mathrm{0.87496\hspace{0.33em}485}$ | $\mathrm{\times 10\u207b\u2078}$ | $\mathrm{0.71668\hspace{0.33em}638}$ | $\mathrm{\times 10\u207b\xb9}$ |

10 | $\mathrm{0.14618\hspace{0.33em}117}$ | $\mathrm{\times 10\u2079}$ | $\mathrm{0.10423\hspace{0.33em}616}$ | $\mathrm{\times 10\u2078}$ | $\mathrm{0.24457\hspace{0.33em}824}$ | $\mathrm{\times 10\u207b\xb9\u2070}$ | $\mathrm{0.65021\hspace{0.33em}292}$ | $\mathrm{\times 10\u207b\xb2}$ |

11 | $\mathrm{0.29154\hspace{0.33em}738}$ | $\mathrm{\times 10\xb9\u2070}$ | $\mathrm{0.37225\hspace{0.33em}201}$ | $\mathrm{\times 10\xb9\u2070}$ | $\mathrm{0.19952\hspace{0.33em}026}$ | $\mathrm{\times 10\u207b\xb9\u2070}$ | $\mathrm{0.58373\hspace{0.33em}946}$ | $\mathrm{\times 10\u207b\xb9}$ |

12 | $\mathrm{0.63994\hspace{0.33em}242}$ | $\mathrm{\times 10\xb9\xb9}$ | $\mathrm{0.37225\hspace{0.33em}201}$ | $\mathrm{\times 10\xb9\u2070}$ | $\mathrm{0.37946\hspace{0.33em}279}$ | $\mathrm{\times 10\u207b\xb9\xb3}$ | $\mathrm{0.44851\hspace{0.33em}3}\underset{\xaf}{87}$ | $\mathrm{\times 10\u207b\xb2}$ |

13 | $\mathrm{0.15329\hspace{0.33em}463}$ | $\mathrm{\times 10\xb9\xb3}$ | $\mathrm{0.19555\hspace{0.33em}304}$ | $\mathrm{\times 10\xb9\xb3}$ | $\mathrm{0.32057\hspace{0.33em}909}$ | $\mathrm{\times 10\u207b\xb9\xb3}$ | $0.49269\underset{\xaf}{383}$ | $\mathrm{\times 10\u207b\xb9}$ |

14 | $\mathrm{0.39792\hspace{0.33em}611}$ | $\mathrm{\times 10\xb9\u2074}$ | $\mathrm{0.19555\hspace{0.33em}304}$ | $\mathrm{\times 10\xb9\xb3}$ | $\mathrm{0.44167\hspace{0.33em}174}$ | $\mathrm{\times 10\u207b\xb9\u2076}$ | $0.327\underset{\xaf}{\mathrm{92\hspace{0.33em}861}}$ | $\mathrm{\times 10\u207b\xb2}$ |

15 | $\mathrm{0.11126\hspace{0.33em}602}$ | $\mathrm{\times 10\xb9\u2076}$ | $\mathrm{0.14186\hspace{0.33em}384}$ | $\mathrm{\times 10\xb9\u2076}$ | $\mathrm{0.38242\hspace{0.33em}250}$ | $\mathrm{\times 10\u207b\xb9\u2076}$ | $0.425\underset{\xaf}{\mathrm{50\hspace{0.33em}628}}$ | $\mathrm{\times 10\u207b\xb9}$ |

16 | $\mathrm{0.33340\hspace{0.33em}012}$ | $\mathrm{\times 10\xb9\u2077}$ | $\mathrm{0.14186\hspace{0.33em}384}$ | $\mathrm{\times 10\xb9\u2076}$ | $\mathrm{0.39924\hspace{0.33em}861}$ | $\mathrm{\times 10\u207b\xb9\u2079}$ | $0.\underset{\xaf}{\mathrm{00000\hspace{0.33em}000}}$ |

Similar considerations apply to the first-order equation

3.6.17 | $${a}_{n}{w}_{n+1}-{b}_{n}{w}_{n}={d}_{n}.$$ | ||

Thus in the inhomogeneous case it may sometimes be necessary to recur backwards to achieve stability. For analyses and examples see Gautschi (1997b).

For a difference equation of order $k$ ($\ge 3$),

3.6.18 | $${a}_{n,k}{w}_{n+k}+{a}_{n,k-1}{w}_{n+k-1}+\mathrm{\dots}+{a}_{n,0}{w}_{n}={d}_{n},$$ | ||

or for systems of $k$ first-order inhomogeneous equations, boundary-value methods are the rule rather than the exception. Typically $k-\mathrm{\ell}$ conditions are prescribed at the beginning of the range, and $\mathrm{\ell}$ conditions at the end. Here $\mathrm{\ell}\in [0,k]$, and its actual value depends on the asymptotic behavior of the wanted solution in relation to those of the other solutions. Within this framework forward and backward recursion may be regarded as the special cases $\mathrm{\ell}=0$ and $\mathrm{\ell}=k$, respectively.