Download Three-Term Recurrence Relations for Special Functions and more Study notes Mathematics for Computing in PDF only on Docsity!
SIAM REVIEW Vol. 9, No. 1, January, 1967
COMPUTATIONAL ASPECTS OF THREE-TERM RECURRENCE
RELATIONS*
WALTER GAUTSCHIf
Introduction. Recurrence relations are one of the basic mathematical tools of
computation. There is hardly a computational task which does not rely on
recursive techniques, at one time or another. The widespread use of recurrence relations can be ascribed to their intrinsic constructive quality, and the great
ease with which they are amenable to mechanization. On the other hand, like
most recursive processes, recurrence relations are susceptible to error growth. Each cycle ot a recursive process not only generates its own rounding errors, but also inherits the rounding errors committed in all the previous cycles. If con-
ditions are unfavorable, the resulting propagation of error may well be dis-
astrous. It is this aspect of recurrence relations--the possibility and the preven-
tion of numerical instability--that will be of concern to us.
The problem of numerical instability has been studied extensively for differ-
ence equations arising in the numerical solution of ordinary and partial dif- ferential equations. In the seemingly much simpler context of a single linear
difference equation, the problem has received only sporadic attention, even
though such difference equations, particularly of the second order, occur promi-
nently in many branches of pure and applied mathematics. We mention, e.g.,
the recurrence relations satisfied by large classes of special functions of mathe-
matical physics and statistics, the three-term recurrence relations that lie at the
heart of continued fraction theory and the theory of orthogonal polynomials,
and the miscellaneous recurrence relations one encounters when constructing
series expansions, asymptotic or otherwise, to solutions of linear differential
equations. We believe, therefore, that a systematic review of some of the compu- tational problems attending recurrence relations might be of value. In the follow-
ing we attempt to present such a survey, restricting attention to the special case
of three-term recurrence relations. The kind of instability we are concerned with, may be described as follows. Consider a^ three-term^ recurrence relation of the form (0.1) y,+l q-^ a,y,, q-.^ b,y,_ O, n^ 1, 2, 3, ...,
where am, b are given sequences of real or complex numbers, and b 0. The
general solution of (0.1) can be spanned by any pair f, g of^ linearly independent
solutions. We are interested in the^ special case where there^ exists such a pair
having the property
-00 g
- (^) Received by the editors (^) February 17, 1966, and in (^) revised form July (^) 18, 1966. Computer Sciences Department, Purdue (^) University, Lafuyette, Indiana, and Argonne National (^) Laboratory, Argonne, Illinois. This work was (^) performed in (^) part under the (^) auspices of the United States Atomic Energy Commission. 24
THREE-TERM RECURRENCE RELATIONS 25
Serious problems then arise if one attempts to compute the solution (^) f or any constant multiple of (^) f.
To see^ this, we first note that (0.2) implies
(0.3) lira f^0 y
for any solution y not proportional to (^) f. Such a solution, indeed, is represent- able in the form y, af, (^) + bg, with b 0, and therefore
lira lim
f/g
O.
y b (^) + a(f/g)
If we now generate (^) f by (0.1), using only approximate initial values yo f0, y (^) f (due to rounding, for example), but recurring with infinite precision, we obtain a solution y (^) which, in general, is linearly independent of (^) f. Therefore, by (^) (0.3), we will have
as
.e., the relative^ error^ of^ ,^ the^ intended^ approximation^ to^ ,^ becomes^ arbitrarily large. Therefore, this^ srahtforward method^ of^ computing^ is^ uterly^ in- effective. Observe that the set of all solutions (^) f. of (^) (0.1) having the property indicated in (0.3) forms^ a^ one-dimensional^ subspce^ of^ the^ spce^ of^ all^ solutions.^ (There cn (^) be no two linearly dependent solutions (^) f, (^) f enjoyg this property, (^) since, otherwise, f/f,’ and^ f/f,^ would^ both^ have^ the^ limit^ zero, as^ n^ ,^ which^ is
absurd.) We^ call^ the^ solutions^ of^ this^ subspace^ minimal^ at^ infinity,^ or^ briefly
minimal. A^ nonminimal^ solution^ will^ be^ referred^ to^ as^ dominant.^ Euch^ dominant solution is^ asymptotically proportional^ to^ g.^ Note^ thut, in^ contrast^ to^ dominant
solutions, minimal^ solution^ is^ uniquely^ determined^ by^ one^ initiul^ value.
To illustrate^ the^ difficulty^ of^ calculatg^ minal^ solutions, consider^ the^ prob-
lem of generating Bessel^ functions^ of^ the^ first^ kind, J,(x), for^ fixed^ x, and^ n 0, 1, 2, ....^ As^ is^ well-known, these functions^ (of n)^ obey^ the^ three-term^ re- currence relation
2n (0.4) y+ y^ + y.-^ 0.
From tbles of^ Bessel functions we^ find, e.g.,^ thut^ for^ x^ 1, J0(1) .7651976866,
J(1) .4400505857, accurute^ to^ ten^ figures.^ Generating^ the^ next^99 values^ of
J(1) on^ digital^ computer^ by^ straightforward^ recursion,^ we^ obtain^ the^ results
The notion of a minimal solution ppears to have first been^ introduced^ by^ Pincherle^ in connection with his generalization of continued frctions [44]. Pincherle^ called^ it^ "dis- tinguished" solution (soluione distint). In the^ theory of^ linear^ differential^ equations^ the term "principal" solution is also in^ use [24]. The minimal^ solution^ can^ often^ be^ identified with a solution "of type II" in terminology of^ Schfke^ [51]. This example is well known, and^ has^ received^ considerable attention^ in^ the^ literature.
See the references in 5.
THREE-TERM RECURRENCE RELATIONS 27
backward recurrence algorithm of J. C. P. Miller. While this algorithm is widely
regarded as (^) iust a "trick of the trade," our presentation will show that it derives naturally from rather elegant results of classical analysis. In 4 two alternate
algorithms are described which are more flexible, but more elaborate, than the first
algorithm. The remaining paragraphs discuss a number of applications, mostly to the computation of higher transcendental functions such as Bessel functions
(5), associated Legendre functions (6), regular Coulomb wave functions
(7), and other miscellaneous functions (^) (8-10). 11 contains an application
to a Sturm-Liouville boundary value problem on an infinite interval.
The extent to which our algorithms are affected by rounding errors will not be
discussed in detail. Our experience seems to indicate, however, that none of the
algorithms is sensitive to rounding, unless the series (0.5) is subject to cancella-
tion of (^) terms. The rigorous analysis of error propagation is an interesting, though
diificult, outstanding problem in this area. A significant contribution in this
direction is due to Olver, who recently analyzed the error accumulation in
Miller’s algorithm [38].
In principle, there are other stable procedures that could be used to calculate
minimal solutions of (0.1). For example, we could set up the boundary value
problem of finding the solution yn of (0.1) which satisfies
y0 (^) f0, y (^) =f
for some sufficiently large N. Clearly, this amounts to solving the linear system
al (^1) Yl
b2 a.^1
b3 a3 1 y
of (^) equations
(^0) aN--
whose matrix is^ tridiagonal. Any^ of the standard methods, such^ as triangular
decomposition methods, may be used to solve (0.6). Unfortunately, the pro-
cedure requires two^ values, (^) f0 and f, of^ the desired solution to^ be^ known in advance. Either^ one may be^ difficult, or time-consuming, to obtain. The problem of computing minimal solutions is^ clearly^ not peculiar to three-
term recurrence^ relations. It^ may equally arise in connection with^ other func-
tional equations, such^ as^ linear homogeneous difference and^ differential^ equa-
tions of arbitrary order, and^ systems of such equations. Whenever^ the space S
of all solutions is the direct sum S $1^ @ S. of two subspaces $1^ and^ S, and every solution s S dominates, in some appropriate sense, over^ every solu- tion s. $2, we may consider S as the set of^ minimal^ solutions^ with^ respect^ to
The author is indebted to Dr. M. E. Rose for pointing this out.
28 WALTER GAUTSCHI
the decomposition S S1 (^ $2,^ (There may be^ several^ such^ decompositions.)
The problem of computing minimal^ solutions^ iu^ this^ sense^ has not^ been^ thor-
oughly studied, though the^ work^ of^ Clenshaw^ [7] and^ Schifke^ [51] suggests that
effective computational methods^ may exist^ also in^ this^ more^ general^ context.
1. Three-term recursion and^ continued^ fractions.^ It^ is^ well-known^ that^ the
concepts of^ three-term recursion^ and^ continued^ fraction^ are^ closely^ related.^ To
every continued fraction^ we^ may in^ fact^ associate^ a^ three-term^ recurrence relation, namely the^ fundamental^ recurrence^ formula^ for^ the^ numerators^ and denominators. Vice^ versa, every^ three-term^ recurrence^ relation^ may^ be^ inter-
preted as^ the^ fundamental^ recurrence^ formula^ for^ some^ continued^ fraction.^ The
first (^) point of view is useful for computing continued (^) fractions, the second for
computing the^ minimal^ solution.^ We^ begin^ by^ considering^ several^ methods^ of
calculating a continued fraction.
Suppose we are given the continued fraction
1.1 al^ a^ aa bl b. (^) -t- b-t-
where the^ partial^ numerators^ am and^ partial^ denominators^ b are^ real^ or^ complex
numbers. Denote^ its^ nth^ numerator^ and^ nth^ denominator^ by^ Am and^ B,,^ re- spectively, so^ that
(1.2) al^ a^ an^ A
bl-f- b-+- b Bn"
The value of the continued fraction^ (1.1), if^ it^ exists, is^ defined^ as^ the^ limit
limn..o A,/B,.^ The^ quantities^ An,^ Bn^ satisfy^ the^ fundamental^ recurrence^ formulas
(see, e.g.,^ [59,^ p.^ 15])
A, bnAn_ -4- a,An_
(1.3)
B, bnBn_ anBn_2,
n 1, 2, 3, ...,
where (1.4) A_^ 1, A0 0;^ B_t^ 0,^ B0 1. This shows that^ am An-1 and^ fin B-1 constitute^ a^ pair^ of^ linearly^ independ-
ent solutions of^ the^ three-term^ recurrence^ relation
(1.5) (^) Yn+ bnYn anyn--1^ 0,^ n^ 1,^ 2,^ 3,^ ....
A first^ method^ of^ computation^ flows^ directly^ from^ these^ fundamental^ recur-
rence relations. Thus, one generates the A’s^ and^ B’s^ recursively,^ by^ means^ of
(1.3) and^ (1.4), and^ concurrently^ the^ ratios^ An/Bn,^ until^ the^ latter^ converge within the^ required^ tolerance.^ As^ An and^ B are^ likely^ to^ grow^ rapidly^ with^ n,
some care^ must^ be^ exercised^ if^ this^ method^ is^ used^ on^ a^ digital^ computer.^ Initial
scaling, and^ possibly repeated^ subsequent^ scaling,^ may^ be^ necessary^ to^ avoid
overflow.
A second method, which^ avoids^ the^ necessity^ of^ scaling,^ consists^ in^ evaluating
30 VALTER GAUTSCItI
the initial values being
(1.10) al
ul 1, vl w
bl"
None of the disadvantages noted in the previous two methods are present
here.
We have seen that the continued fraction (1.1) leads naturally to the three-
term recursion (^) (1.5). Suppose (^) now, conversely, that we are given a three-term recurrence relation
(1.11) y+l -t- any^ -t- by_^ O, bn 0, n^ 1, 2, 3, ....
Define (^) an, n to be the special solutions of (^) (1.11) with initial values
(1.12) a0 1, al 0; fl0 0, fl 1.
Then, evidently,^ A (^) On+l and^ B (^) +1 are^ the^ numerators^ and^ denominators, respectively, of the continued fraction
(1.13) -b^ b^ ba
al a-- a--
which is equivalent to the continued fraction
(1.14) bl^ b^ b^ ....
a- a2-- a3-
We may formally arrive at this continued fraction also in the^ following way.
Let us introduce the ratios
Yn+l
rn
Yn
Dividing (1.11) by yn then gives
r+a+
bn
rn-- from which
rn-
an
rn
n 0,1,2, ....
Applying this^ formula^ repeatedly, with^ n^ successively^ increasing,^ we^ get
yn --bn b+
(1.15) (^) r-
Yn-1 an--^ an+l-- an+2--
In particular, when^ n^ 1,
y --b^ b^ b Y0 al--^ a2--
This derivation indicates that the continued fraction (1.14), and^ similarly
THREE-TERM RECURRENCE RELATIONS 31
the continued fractions in (1.15), are related to ratios of consecutive values
for some solution y. The argument, however, neither insures us of the con-
vergence of these continued fracliions, nor does it tell us for what particular solu-
tion the ratios are to be formed. These matters are clarified by the following theorem. THEOREM 1.1 (Pincherle [45]). The continued (^) fraction (1.14) converges (^) "if and only (^) if the recurrence relation (1.11) possesses a minimal solution (^) f with (^) fo O. In case (^) of convergence, moreover, one has
f, --b bn+l bn+
n 1,2,3,...,
f-i^ an--^ an+l-- an+2--
provided (^) f O (^) for n 0,1,2, .... Proof. (a)^ Assume^ the^ continued^ fraction^ in^ (1.14)^ converges.^ Then^ so^ does
the equivalent continued fraction (1.13). Therefore
C
where an, fin are the solutions of (1.11) defined by the^ initial values (1.12), and
c is some constant. Let
(1.17) (^) f a- (^) c. Take any other solution of (1.11), say y aa (^) + b. Then^ ac^ + b^ 0, and
lim f^ =lim a-c^ lim (a/)^ -c =0. y (^) n aa (^) + b, (^) a(a./) (^) + b This shows^ that^ the^ solution^ f defined^ in^ (1.17) is^ a^ minimal^ solution^ of^ (1.11). Moreover, (^) f0 a0 0. (b) Assume^ now^ that^ (1.11) possesses a^ minimal^ solution, (^) f say,^ for^ which f0 0.^ Then
We note^ that^ @ is^ not^ a^ constant^ multiple^ of^ f, since^ f0 0.^ Therefore, (^) f
being minimal,
lira (^) ,f^ f0 lira a+f 0,
and so
lirn a^ _f f0"
This establishes convergence of the continued fraction (1.13), and thus of that
in (^) (1.14), and also proves (^) (1.16) for n 1. To prove (^) (1.16) for general n (^) > 1, we need only observe that (^) zm f+m-1, considered as a function of (^) m, is a minimal solution of
z,+ na^ an+z,^ -P bn+m_lZm^ O,^ n^ 1,^ 2,^ 3,
THREE-TERM RECURRENCE RELATIONS 33
having again used (1.20). The second ease is verified similarly, using (1.20)
in the form
() 1 () () V (a+v+^ v+.;. bn+l
Since we already established (^) (1.21) for n <_- (^) 1, it now follows by induction
that the result holds for n /c,/c 0, 1, 2, 3, that is, for all integers
n, ,.^ Theorem^ 1.2^ is^ proved. We note that relation (^) (1.21), for (^) > n, can also be obtained from the known fact (^) (el. [43, vol. (^) I, p. (^) 3]) that (^) /-z()^ nd ()re the numerators nd
denominators, respectively, of^ the^ continued^ fraction
where
b’-i
am 1
am b’ b
Alternatively, Theorem 1.2^ may^ be^ obtained, as^ a^ special^ case, from^ the^ known
result that "multipliers" of^ a^ linear^ difference^ equation^ satisfy^ the^ adjoint
difference equation (cf. [35, 12.6]).
- Some results from^ the asymptotic theory of^ ear^ second^ order^ difference
equations. In^ applications of^ Theorem^ 1.1, it^ is^ in^ general^ easier^ to^ recognize^ a
given solution^ of^ a^ three-term^ recurrence^ relation^ to^ be^ minimal^ than^ to^ establish
convergence of^ the^ corresponding continued^ fraction.^ One^ is^ aided^ in^ this^ by
classicM results^ from^ the^ asymptotic theory^ of^ derence^ equations,^ notably^ by^ a
theorem of^ Poincar6, and^ by^ refinements^ and^ extensions^ thereof^ due^ to^ Perron
and Kreuser. For^ convenience^ of^ the^ reader, we^ are^ recalling^ here^ these^ theorems
for the special case^ of^ a^ second-order^ difference^ equation
(2.1) y+^ ay^ by_^ 0, n^ 1, 2, 3, ....
We assume, throughout,^ that
(2.2) b0, n=^ 1,2,3, ....
We begin with the case where the coefficients a and b in (2.1) have finite
limits
(2.3) a a, b b,^ n^ ,
not (^) excluding that b 0. One then calls (^) (2.1) a Poincar$^ difference equation, and calls
(2.4) (t) (^) + at^ + b
the characteristic^ polynomial of^ (2.1). As^ may^ be^ expected,^ the^ solutions^ of (2.1) behave silarly, for^ large^ n, to^ the^ solutions^ of^ the^ difference^ equation (2.1) with^ cstant^ coefficients^ a a, b b.^ This^ is^ borne^ out^ by^ the^ following two theorems.
34 WALTER GAUTSCHI
TtIEOREM 2.1 (^) (Poincar [46]). (^) If the characteristic (^) polynomial (2.4) (^) of (2.1) has zeros (^) tl t2 (^) of distinct (^) moduli,
(2.5) (^) It1 > It2 !,
then for every nontrivial (^) solution y, (^) of (2.1) we (^) have
(2.6) lim Y+-----^ tr, r^ 1, or^ r^ 2. Yn THEOREM 2.2 (^) (Perron [41]). Under the assumption (^) of Theorem 2.1 there exist two (^) linearly independent solutions y, and y,2 (^) of (2.1) such that
(2.7) lim^ Y"+’r^ #,
n--> Yn,r Theorem 2.2 (^) implies that
f
is a minimal solution of (2.1). To see this, choose r and r such that
r= 1,2.
which under the assumption (2.5) is certainly possible. By^ (2.7) we then have,
for n sufficiently large,
Yn+l, (^) <
2 n^ >^ no
Hence
and
Y,I This (^) shows that
/ \n--no
n>no.
from which the assertion (^) follows.
We also^ note^ that^ in^ (2.6) one^ has^ r^2 for^ the^ minimal^ solution, and^ r^1 for
any other^ solution.
We shall require a generalization of Theorem 2.2 relating to a difference equa-
tion (2.1) whose coefficients satisfy
(2.8) a,-an", b,bn^ ,^ ab^ 0; a, real; n-.
The (^) asymptotic structure of the solutions now (^) depends on the (^) Newton-Puiseux
diagram formed with^ the^ points P0(0, 0), P(1, a), P(2, ). This^ is^ the broken
line (^) PoPP if (^) P is above the (^) straight line (^) joining P0 with (^) P2" otherwise it
36 WALTER GAUTSCHI
(3.5)
Similarly, we have
proof of Perron’s theorem is given in [14], and reproduced in [34]. Far-reaching
generalizations, and simplified proofs, of all these theorems, including Kreuser’s
theorem, were^ recently^ obtained^ in^ [51].
3. A first algorithm for computing the minimal solution. We assume now that
the recurrence relation
(3.1) (^) Yn+I a,y (^) bnYn-1 O, n 1, 2, 3,’’’, has a nonvanishing minimal (^) solution, (^) fn. We wish to (^) calculate (^) f for n (^) 0, 1, 2, N. In order to (^) specify (^) f uniquely, we (^) can impose one (^) condition, for example prescribe the value of (^) f0 For later (^) applications, we consider the more general normalization
(3.2) k,cm s, s 0,
m
where s nd 0, h, re given quantities, and the series is known to converge.
We do not exclude that (^) X 0 for 11 m (^) > 0, in which cse (^) (3.2) mounts (^) to prescribing (^) f In (^) sense, (3.2) represents the most general linear condition that (^) my be im-
posed. A class of nonlinear conditions will lso be considered briefly.
To introduce the lgorithm, let
_f+ 1 (3.3) (^) rn (^) f .% (^). Xf.
Suppose first that r, s are known for some value n N. The desired
solution f, n 0( 1 )N, can then be obtained as follows.
From Theorem 1.1 we know that
-b b+ b+
(3.4) r- n 1, 2, 3, .-..
Hence, we^ can^ generate the^ ratios^ r for 0 N n (^) < as in (1.6)-(1.8) by
rn-1 u, u^ 1,^ ...,^ 1. an rn
so that
( )
m=n-t--
(3.6) (^) 8n--1 rn--l()n (^) 21- 8n), ,^ ,^ 1,’’’, l.
Hence, also^ the^ quantities s for^0 -<_^ n (^) < ,^ and^ thus^ in particular so, can be ob-
The assumption of (^) f to be nonvanishing is no serious restriction from the (^) practical point of view. This is further discussed at the end of this paragraph.
THREE-TERM RECURRENCE RELATIONS 37
tained recursively. Using (3.2) we now have
.. Xf,^
(s (^) ,ofo)
and so
8 X0 (^) + So
This gives us the initial value of the desired solution. The remMning wlues cn
now be obtained immediately from
fn rn-’n-1,^ n^ 1,^ 2,-..,^ N.
The actual algorithm follows this procedure very closely, except that for the
infinite continued fraction, and the infinite series, representing rn_l and sn, re-
spectively, we now substitute truncated continued fractions, and truncated series.
More precisely, we define
(3.7) (^) r () O, r ()n-^ --bn^ bn+^ b 1 _-<n-<v, an-- (^) an+l-- a
(3.8) s ()
0, sn
() (^) Xrn()r (^) ()n+ r()- (^0) m----n+l =<^ n^ <^.
One then^ verifies readily that^ the formulas (3.5), (3.6) continue to hold if (^) rn is replaced by (^) rn (), and^ s by sn ()^ throughout. Hence the following set^ of recursions arises naturally,
r
() O rn- an (^) + rn ()^ n v, v 1,..., 1, () (^) r() ()
(3.9) s
() O, sn- n-l(M (^) + sn ),
fo(V)
8 fn(v) (v) rn-lan-l 1,^ 2,^ ", N. ho (^) + So () While our initial procedure gave us the^ exact values (^) fn of^ the^ minimal solution, the quantities (^) fn (’)^ now^ derived^ are^ at^ best^ approximations to^ fn. It^ remains^ to successively improve (^) fn (’)^ by repeating (3.9) for^ a^ sequence of^ increasing values^ of
. (^) The complete algorithm for (^) computing the minimal solution (^) may thus be de- fined (^) as follows: Step 1:^ Select an^ integer (^) => N, and^ let^ n ()
0, n^ 0, 1, N.
Step 2:^ Calculate (^) fn (), n 0, 1, N, according to^ the^ formulas^ in (3.9). Step 3: If the N^ + 1 values of^ fn ()^ obtained^ in Step 2 do^ not^ agree with^ the current values of (^) Cn (’)^ to within the^ desired accuracy, then^ redefine^ (’) by n (’) fn (’),^ n^ 0,^ 1,^ N,^ increase^ by^ some^ fixed^ integer,^ say^ 5,^ and^ repeat Step 2; otherwise accept (^) fn ()^ as the final approximations to fn, n 0, 1, N. We note that in the^ special case k0 1, (^) k X 0, all sn (’)^ vanish, so
that the recursion for s
_ in (3.9) may be omitted. Moreover, s f0, and there-
THREE-TERM RECURRENCE RELATIONS 39
While algorithm (3.9) and Miller’s algorithm are mathematically equivalent,
they have^ different^ computational characteristics. In many cses, e.g., the^ qunti-
ties (^) v() grow rapidly as (^) increases, and my cause "overflow" on digital com- purer. In contrast (^) to this, the (^) quantity (^) r ()^ in (3.9) converges to a finite limit as ,
--
, (^) and (^) so does (^) s ()^ if the (^) algorithm converges at all. We now use (3.12) to discuss convergence s ,
of the lgorithm (3.9).
Let g denote any solution of the difference equation (3.1) other thn f, so that
(3.14) lira
--
gt Clearly, () a^ ()f (^) + b ()^ gn
for some constants a (), b ().^ By (3.11) we must have
The first of these relations gives b ()^ (f+z/gz)a (), so that
() a () (f
f+l ) --
g+g
Substituting in (3.12), and^ simplifying, we^ obtain
f(l
f+ g) (3.15) f() g+
1-! (^2) xf f+^ Xgm
In view^ of^ (3.14) nd^ the^ convergence of^ the^ infinite^ series^ in^ (3.2), it^ is^ clear
that (^) lim f(’) (^) f if^ and^ only^ if
(3.16) lira^ f+^ hg 0.
We hve proved the following theorem.
TonE 3.1. Suppose the recurrence relation (3.1) has a nonvanishing minimal
solution, (^) f for which^ (3.2) holds.^ Let^ g^ be^ any^ other^ solution^ of (3.1). Then^ the
algorithm (3.9) converges^ in^ the^ sense
if and^ only^ if (3.16)^ is^ satisfied.
Condition (3.16) holds, e.g., if the k’s are uniformly bounded, and
g+l (^) ___> (^) tl f+l g,
Itl >^ Itl,^ t,,.^ <^.
40 WALTER GAUTSCHI
If all but a finite number of the },’s are zero, then (3.16) is a consequence of
(3.14). Theorem^ 3.1, in^ this^ case, has^ been^ noted^ previously in^ [16].
It is useful to observe that (^) convergence of the (^) algorithm (3.9), in the sense of
Theorem 3.1, implies that
(3.17) r ()^ -+^ rn, s ()
--
s,
,
where rn, s are the quantities defined in (3.3). The first of these relations follows
directly from (^) (3.4) and (3.7). The second follows by induction on n. (^) Indeed, if n (^) 0, we (^) have from (^) the third line in (3.9), () () s^00 s^ 0f So
f0(
-- f
So,
--
.
Assuming now (^8) n--
S-I, we^ get from^ the^ second^ line^ in^ (3.9), and^ from^ (3.6), that 8() ,) rn-1 rn- In case of convergence of the algorithm (3.9), we may obtain^ from^ (3.15) the
following approximate expression for the^ relative error, valid for^ sufficiently
large,
f() (^) f 1 Xmfm -t-^ f+__L (^) Xm gm f+---2 g--
It is interesting to examine what effect the^ modification (3.10) of^ algorithm (3.9) will have upon the relative error (3.18). We^ assume that
where r, s are defined by (3.3), and e, w are^ small^ numbers.^ Then^ a^ simple
computation will show that in place of (3.15) we^ now have
() =A
8g+l m=O
Since pg/g+lis usually substantially smaller^ than^1 (at least^ for^ large^ ,) we
see that the modification (3.10) reduces the^ relative^ error^ of^ f(’)^ effectively^ by^ a factor of ,^ I, or (^) w I, whichever is^ lurger. Hence, our statement^ made^ earlier^ that the convergence of (^) f ()^ to (^) f is^ fster^ the better p approximates r, and^ ,^ approx- imates s, is clearly vindicated.
It is tempting to try substitution of^ the^ type
(3.19) F c,f,,^ c,^ 0, to exert influence upon the convergence criteria (3.14) and^ (3.16). We^ note, how- ever, that^ these^ criteria^ are^ invariant^ with^ respect^ to^ any^ linear^ substitution^ of the form (3.19).
42 WALTER GAUTSCHI
culties. By the following, admittedly superficial, considerations we wish to show
that the presence, or proximity, of such zeros need be^ of no great concern. Suppose, indeed, that (^) f-i is very small in modulus, compared to (^) f. For definiteness, let n (^) > 1. Then, by (3.3), rn_l is very large, and so is r (^) n- I, when , (^) is sufficiently large. From the first line in (3.9) it follows that (^) an (^) + rn()[ must be very small compared to (^) b (^) I. Since neither (^) an nor (^) r ()^ will normally be small, this means that many digits will cancel when the sum (^) a (^) + r ()^ is formed, and so () r_ is not only very large, but also very inaccurate in terms of significant digits. Consequently, (^) r-() will be very small, and also inaccurate. However, r_
--b_/(a_ ()
- r_^ (if^ n^2 will^ again^ be^ accurate,^ since^ a_^ in^ the^ denomi- nator picks up lost accuracy, r
_
being normally much smaller than a._. Later
on, in the formation of the final results,()^ r()^ ()
_ -_ will come out very small and inaccurate, as one must expect. The really questionable point is the computa- tion of ()^ r () (^) #) ()
n--1 jn--1 since^ r_x^ is^ large^ and^ j-x is^ small,^ and^ both^ are^ inac-
curate. We note, however, that
() (^) b n--1 rn--2 (^) In--2 ’n--1 (^) () a.. 1/r^ () () (^) () which shows that the lrgeness of (^) rn- SVeS (^) fn from becoming (^) inaccurate, even though ()^ ()^ ()
r- is.^ A^ similar^ reasoning^ applies^ to^ s-, s-.
More serious is possible loss of ccurcy in the clcultion of ()
(^0) () s^ this would ffect 11 subsequent ().^ It could indeed occur that (^0) + s0 is smll in comparison with 0 , so that mny digits cncel when (^0) + s0 ()^ is formed. The resulting vlue of ()^ would then be quite inaccurate. The sme difficulty might 0 rise if 0 0. Suppose, indeed, that (p (^) > 0) is the first nonvnishing co-
efficient in the series (3.2),
X O,^ X^ O,^0 m^ <^ p, nd that () Ik + s^ hppens^ to^ be^ very^ smll^ compared^ to^ lX.^ Then^ s_^ is necessarily inaccurate, nd this in.ccuracy will be transmitted to all subsequent () () (^) () () 8n--i nd^ finally^ to^0 (),^ in^ view^ of^ the^ relations^ s- r-s n^ p^1 p (^) 2, 1, nd (^) f0 ()^ S/So (). Now for large ,,^ and p 0, we have 1 f (^) =+
s
so that (^) I(X (^) + s())/X is small if (^) ls/(Xf,)] is small. (^) Hence, dangerous cancel- lation occurs when s is small in absolute value compared to the (^) first nonvanishing
term Xf, in (3.2), i.e., when cancellation occurs in the series (3.2) itself. For this
reason, some^ care^ must^ be^ exercised^ in^ the^ selection^ of^ the^ identity^ (3.2).
4. Second and third algorithm for comput.g the mimal solution. The effec-
tiveness of our first algorithm (3.9) is somewhat limited if no reasonable estimate
of the starting value of n is known a priori. The recursions in (3.9) must then
be (^) repeated with increasing values of ,, until sufficient agreement is obtained between successive results (^) fn (’), for all n (^) 0, 1, N. This (^) disadvantage can
THREE-TERM RECURRENCE RELATIONS 43
be removed, at the expense of a more complex algorithm, by making use of the
duality theorem 1.2, or, alternatively, by evaluating a sequence of continued
fractions (3.4). The corresponding algorithms will now be developed. The first
of these, though not in the form given here, is due to Shintani [52].
As was noted in the previous paragraph, we can obtain r, s reeursively, for
(^0) =< n (^) < N, and hence also (^) fn for 0 __< n (^) N N, once (^) rn, s are known. In the fol- lowing we derive a method for calculating (^) r, sN reeursively. If (^) f0 is (^) known, the
s will^ not^ be^ required,^ and^ the^ algorithm^ then^ reduces^ to^ one^ suggested^ by^ G.
Blanch ([4, p. 405 ff]) in connection with Bessel functions.
As for rn, we may simply evaluate the continued fraction
(4.1) -b+l^ bN+^ b+ r aN+l aN+2-- ar+a--
by either the first, or third method described in (^) 1. In the first case we have
(4.2) (^) r lim A
k- I^ k
where
A_t (^) 1, A0 0; B_ (^) 0,
(4.3) (^) A ar+A_ (^) br+A
_
B ar+B_ br+B_, In the second case w6 have
B0 1;
It= 1, 2, 3,....
(4.4) r lim^ w,,
where the w’s are generated as follows"
b+l
ul 1 vl Wl aar+l 1 (4.5) u+--^ 1- (^) (bv++/ar+av++)u’
v+ v(u+^ 1), 1,2,3, ....
For the computation of s, we make use of the fact that
s lims
where (^) s (’)^ is defined by (3.8). The quantities s (), N, may be^ obtained re-
cursively as follows. From the definition (3.8) of
’
and from (3.13), we note
that
(4.7) (^) v (’)s^ ()^ Xv,, (") m----N+
