Digital Library of Mathematical Functions
About the Project
NIST
3 Numerical MethodsAreas

§3.6 Linear Difference Equations

Contents

§3.6(i) Introduction

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

3.6.1 anwn+1-bnwn+cnwn-1=dn,

or equivalently,

3.6.2 anΔ2wn-1+(2an-bn)Δwn-1+(an-bn+cn)wn-1=dn,

where Δwn-1=wn-wn-1, Δ2wn-1=Δwn-Δwn-1, and n. If dn=0, n, then the difference equation is homogeneous; otherwise it is inhomogeneous.

Difference equations are simple and attractive for computation. In practice, however, problems of severe instability often arise and in §§3.6(ii)3.6(vii) we show how these difficulties may be overcome.

§3.6(ii) Homogeneous Equations

Given numerical values of w0 and w1, the solution wn of the equation

3.6.3 anwn+1-bnwn+cnwn-1=0,

with an0, n, can be computed recursively for n=2,3,. 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 wn and of an independent solution gn, 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 wn is a nontrivial solution such that

3.6.4 wn/gn0,
n.

Then wn is said to be a recessive (equivalently, minimal or distinguished) solution as n, and it is unique except for a constant factor. In this situation the unwanted multiples of gn grow more rapidly than the wanted solution, and the computations are unstable. Stability can be restored, however, by backward recursion, provided that cn0, n: starting from wN and wN+1, with N large, equation (3.6.3) is applied to generate in succession wN-1,wN-2,,w0. The unwanted multiples of gn now decay in comparison with wn, hence are of little consequence.

The values of wN and wN+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 wN and wN+1 to be known in advance. These are described in §§ 3.6(iii) and 3.6(v).

§3.6(iii) Miller’s Algorithm

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 wN and wN+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 gn die away. It therefore remains to apply a normalizing factor Λ. 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 Λ can be the true value of w0 divided by its trial value, or Λ can be chosen to satisfy a known property of the wanted solution of the form

3.6.5 n=0λnwn=1,

where the λ’s are constants. The latter method is usually superior when the true value of w0 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.

§3.6(iv) Inhomogeneous Equations

Similar principles apply to equation (3.6.1) when ancn0, n, and dn0 for some, or all, values of n. If, as n, the wanted solution wn 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, the asymptotic behavior of wn is intermediate to those of two independent solutions fn and gn of the corresponding inhomogeneous equation (the complementary functions). More precisely, assume that f00, gn0 for all sufficiently large n, and as n

3.6.6 fn/gn 0,
wn/gn 0.

Then computation of wn by forward recursion is unstable. If it also happens that fn/wn0 as n, then computation of wn by backward recursion is unstable as well. However, wn can be computed successfully in these circumstances by boundary-value methods, as follows.

Let us assume the normalizing condition is of the form w0=λ, where λ is a constant, and then solve the following tridiagonal system of algebraic equations for the unknowns w1(N),w2(N),,wN-1(N); see §3.2(ii). Here N is an arbitrary positive integer.

3.6.7 [-b1a10c2-b2a2cN-2-bN-2aN-20cN-1-bN-1][w1(N)w2(N)wN-2(N)wN-1(N)]=[d1-c1λd2dN-2dN-1].

Then as N with n fixed, wn(N)wn.

§3.6(v) Olver’s Algorithm

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 wn of the inhomogeneous equation (3.6.1) for which the conditions of §3.6(iv) are satisfied.

Suppose again that f00, w0 is given, and we wish to calculate w1,w2,,wM to a prescribed relative accuracy ϵ for a given value of M. We first compute, by forward recurrence, the solution pn of the homogeneous equation (3.6.3) with initial values p0=0, p1=1. At the same time we construct a sequence en, n=0,1,, defined by

3.6.8 anen=cnen-1-dnpn,

beginning with e0=w0. (This part of the process is equivalent to forward elimination.) The computation is continued until a value N (M) is reached for which

3.6.9 |eNpNpN+1|ϵmin1nM|enpnpn+1|.

Then wn is generated by backward recursion from

3.6.10 pn+1wn=pnwn+1+en,

starting with wN=0. (This part of the process is back substitution.)

An example is included in the next subsection. For further information, including a more general form of normalizing condition, other examples, convergence proofs, and error analyses, see Olver (1967a), Olver and Sookne (1972), and Wimp (1984, Chapter 6).

§3.6(vi) Examples

Example 1. Bessel Functions

The difference equation

3.6.11 wn+1-2nwn+wn-1=0,
n=1,2,,

is satisfied by Jn(1) and Yn(1), where Jn(x) and Yn(x) are the Bessel functions of the first kind. For large n,

3.6.12 Jn(1)1(2πn)1/2(2n)n,
3.6.13 Yn(1)(2πn)1/2(2n)n,

10.19(i)). Thus Yn(1) is dominant and can be computed by forward recursion, whereas Jn(1) is recessive and has to be computed by backward recursion. The backward recursion can be carried out using independently computed values of JN(1) and JN+1(1) or by use of Miller’s algorithm (§3.6(iii)) or Olver’s algorithm (§3.6(v)).

Example 2. Weber Function

The Weber function En(1) satisfies

3.6.14 wn+1-2nwn+wn-1=-(2/π)(1-(-1)n),

for n=1,2,, and as n

3.6.15 E2n(1) 2(4n2-1)π,
3.6.16 E2n+1(1) 2(2n+1)π;

see §11.11(ii). Thus the asymptotic behavior of the particular solution En(1) is intermediate to those of the complementary functions Jn(1) and Yn(1); moreover, the conditions for Olver’s algorithm are satisfied. We apply the algorithm to compute En(1) to 8S for the range n=1,2,,10, beginning with the value E0(1)=-0.56865  663 obtained from the Maclaurin series expansion (§11.10(iii)).

In the notation of §3.6(v) we have M=10 and ϵ=12×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 wn for n=1,2,,10 are the wanted values of En(1). (It should be observed that for n>10, however, the wn are progressively poorer approximations to En(1): the underlined digits are in error.)

Table 3.6.1: Weber function wn=En(1) computed by Olver’s algorithm.
n pn en en/(pnpn+1) wn
0 0.00000 000 -0.56865 663 -0.56865 663
1 0.10000 000 ×10¹ 0.70458 291 0.35229 146 0.43816 243
2 0.20000 000 ×10¹ 0.70458 291 0.50327 351 ×10⁻¹ 0.17174 195
3 0.70000 000 ×10¹ 0.96172 597 ×10¹ 0.34347 356 ×10⁻¹ 0.24880 538
4 0.40000 000 ×10² 0.96172 597 ×10¹ 0.76815 174 ×10⁻³ 0.47850 795 ×10⁻¹
5 0.31300 000 ×10³ 0.40814 124 ×10³ 0.42199 534 ×10⁻³ 0.13400 098
6 0.30900 000 ×10⁴ 0.40814 124 ×10³ 0.35924 754 ×10⁻⁵ 0.18919 443 ×10⁻¹
7 0.36767 000 ×10⁵ 0.47221 340 ×10⁵ 0.25102 029 ×10⁻⁵ 0.93032 343 ×10⁻¹
8 0.51164 800 ×10⁶ 0.47221 340 ×10⁵ 0.11324 804 ×10⁻⁷ 0.10293 811 ×10⁻¹
9 0.81496 010 ×10⁷ 0.10423 616 ×10⁸ 0.87496 485 ×10⁻⁸ 0.71668 638 ×10⁻¹
10 0.14618 117 ×10⁹ 0.10423 616 ×10⁸ 0.24457 824 ×10⁻¹⁰ 0.65021 292 ×10⁻²
11 0.29154 738 ×10¹⁰ 0.37225 201 ×10¹⁰ 0.19952 026 ×10⁻¹⁰ 0.58373 946 ×10⁻¹
12 0.63994 242 ×10¹¹ 0.37225 201 ×10¹⁰ 0.37946 279 ×10⁻¹³ 0.44851 387¯ ×10⁻²
13 0.15329 463 ×10¹³ 0.19555 304 ×10¹³ 0.32057 909 ×10⁻¹³ 0.49269383¯ ×10⁻¹
14 0.39792 611 ×10¹⁴ 0.19555 304 ×10¹³ 0.44167 174 ×10⁻¹⁶ 0.32792 861¯ ×10⁻²
15 0.11126 602 ×10¹⁶ 0.14186 384 ×10¹⁶ 0.38242 250 ×10⁻¹⁶ 0.42550 628¯ ×10⁻¹
16 0.33340 012 ×10¹⁷ 0.14186 384 ×10¹⁶ 0.39924 861 ×10⁻¹⁹ 0.00000 000¯

§3.6(vii) Linear Difference Equations of Other Orders

Similar considerations apply to the first-order equation

3.6.17 anwn+1-bnwn=dn.

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 (3),

3.6.18 an,kwn+k+an,k-1wn+k-1++an,0wn=dn,

or for systems of k first-order inhomogeneous equations, boundary-value methods are the rule rather than the exception. Typically k- conditions are prescribed at the beginning of the range, and conditions at the end. Here [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 =0 and =k, respectively.

For further information see Wimp (1984, Chapters 7–8), Cash and Zahar (1994), and Lozier (1980).