--page0001--
Editorial Board
A. Jeffrey, Engineering Mathematics,
University of Newcastle upon Tyne (Main Editor)
R. Aris, Chemical Engineering and Materials Science,
University of Minnesota
W. Burger, Institut fur Theoretische Mechanik,
Universitat Karlsruhe
J. Douglas Jr, Mathematics,
University of Chicago
K. P. Hadeler, Institut fur Biologie,
Universitat Tubingen
W. F. Lucas, Operations Research and Industrial Engineering,
Cornell University
J. H. Seinfeld, Chemical Engineering,
California Institute of Technology
Computation with
Recurrence Relations
Jet Wimp
Drexel University, Philadelphia
It
Pitman Advanced Publishing Program
BOSTON • LONDON • MELBOURNE
.
--page0002--
PITMAN PUBLISHING LIMITED
128 Long Acre, London WC2E 9AN
PITMAN PUBLISHING INC
1020 Plain Street, Marshfield, Massachusetts 02050
Associated Companies
Pitman Publishing Pty Ltd, Melbourne
Pitman Publishing New Zealand Ltd, Wellington
Copp Clark Pitman, Toronto
First published 1984
© Jet Wimp 1984
AMS Subject Classifications: (main) 65D15
(subsidiary) 65D99
library of Congress Cataloging in Publication Data
Wimp, Jet.
Computation with recurrence relations.
(Applicable mathematics series)
Bibliography: p.
Includes index.
1. Functional differential equations. 2. Point
mappings (Mathematics). 3. Approximation theory.
I. Title. II. Series.
QA431.W64 1984 515.3'5 83-13226
ISBN 0-273-08508-5
British Library Cataloguing in Publication Data
Wimp, Jet
Computation with recurrence relations.—
(Applicable mathematics)
1. Recursive functions
I. Title II. Series
511.3 QAG.615
ISBN 0-273-08508-5
All rights reserved, No part of this publication may be reprodueed, stored in a retrieval
system, or transmitted, in any form or by any means, electronic, mechanical, photocopying,
recording and/or otherwise without the prior written permission of the publishers. This book
may not be lent, resold, hired out or otherwise disposed of by wiiy of trade in any form of
binding or cover other than that in whieh it is published, without the prior consent of the
publishers.
Filmset and printed in Northern Ireland at The Universities Press (Belfast) Ltd, and bound
at the Pitman Pretts. Bath, Avon.
Dedication
Dedicated with profound gratitude to that loving and anonymous fellowship
whose gifts are spirit, sanity, life itself.
.
--page0003--
Acknowledgement
I wish to thank Drexel University, which has been generous in its support
of the writing of this book, and Alison Chandler, whose typing skills
resulted in such a beautiful manuscript. Janet Heilman proofread both
manuscript and page proofs, and her expertise saved me from much of
the travail that usually attends the production of a book. Mark D. Cain
has developed the graphics illustrating strange attractors.
Contents
Preface xi
List of symbols xiii
1 Introduction 1
2 General results on the forward stability of recursion relations 4
2.1 Background 4
2.2 Homogeneous systems 6
2.3 Nonhomogeneous systems 11
2.4 The first-order scalar case; forward vs. backward recursion 12
2.5 The computation of successive derivatives 14
2.6 Scalar equations of higher order: minimal and dominant
solutions 17
3 First-order equations used in the backward direction: the Miller
algorithm 23
3.1 Introduction: the algorithm 23
3.2 Convergence and error analysis 25
4 Second-order homogeneous equations: the Miller algorithm 29
4.1 The algorithm 29
4.2 Reduction of the error by the use of asymptotic
information 37
4.3 Error analysis, the simplified algorithm: case of negative
coefficients 38
4.4 Error analysis, the general algorithm 45
4.5 Algorithms based on continued fractions 46
4.6 Minimal solutions and orthogonal polynomials 53
4.7 The Clenshaw averaging process 57
5 Applications of the Miller algorithm to the computation of special
functions 61
5.1 The confluent hypergeometric function (a, c;x) 61
5.2 The confluent hypergeometric function ty(a,c\\x) 63
.
--page0004--
viii Contents
5.3 The Gaussian hypergeometric function 70
5.4 Associated Legendre functions 73
5.5 The Legendre function Q?(z) 75
5.6 The Jacobi polynomials P^T1\"iTl)(-itD) 75
5.7 Bessel functions 78
5.8 Zeros of Bessel functions 82
5.9 Eigenvalues of Mathieu's equation 82
6 Second-order nonhomogeneous equations: the Olver algorithm 86
6.1 Introduction: the algorithm 86
6.2 Solution by forward elimination (Method A) 90
6.3 The method of averaging (Method B) 94
6.4 The LU-decomposition (Method C) 96
6.5 Adaptation to general normalizing conditions 98
6.6 Conclusions 103
7 Higher-order systems: homogeneous equations 104
7.1 Observations on higher-order equations 104
7.2 The Miller algorithm 105
7.3 The matrix formulation: stability and weak stability 114
7.4 The Clenshaw averaging process 122
7.5 Topics for future research: infinite systems 128
7.5.1 Basic series for functions satisfying functional
equations 128
7.5.2 Stieltjes moment integrals 131
8 Higher-order systems (continued): the nonhomogeneous case 135
8.1 Introduction 135
8.2 The Wimp-Luke method 137
8.3 The Lozier algorithm 139
9 The computation of 3F2A) 152
9.1 The recursion 152
9.2 The algorithm; truncation error 154
9.3 Computing the Beta function 157
9.4 Another 3F2A) 158
10 Computations based on orthogonal polynomials 161
10.1 Preliminaries: properties of some orthogonal
polynomials 161
10.1.1 Chebyshev polynomials 161
10.1.2 Jacobi polynomials 163
10.2 Evaluation of finite sums of functions which satisfy a linear
homogeneous recurrence 16.S
Contents ix
11
12
13
10.2.1 The algorithm 165
10.2.2 Error analysis; three-term recurrence 167
10.2.3 Converting one expansion into another 172
Series solutions to ordinary differential equations 177
11.1 Taylor series solutions 178
11.2 The construction of general recurrence relations for the
coefficients of Gegenbauer series 186
11.2.1 Introduction and basic formulas 186
11.2.2 The algorithms of Clenshaw and Elliott 188
11.2.3 The Lewanowicz construction 196
11.3 Chebyshev series solutions for nonlinear differential
equations 202
Multidimensional recursion algorithms; general theory 209
12.1 Background: convergence properties of complex
sequences 209
12.2 Introduction; invariants 211
12.3 Divergence; strange attractors 219
12.4 Mean values 223
Two-dimensional algorithms 229
13.1 General remarks: invariant curves 229
13.2 Evaluation of certain infinite products 231
13.3 Carlson's results 237
13.4 An algorithm of Gatteschi 245
13.5 Tricomi's algorithm 248
13.6 Solutions of linear functional equations 249
14 Higher-dimensional algorithms 254
14.1 The computation of a class of trigonometric integrals 254
14.2 Incomplete elliptic integrals 260
Appendix A The general theory of linear difference equations 265
Appendix B The asymptotic theory of linear difference
equations 270
B.I General theory 270
B.2 The construction of formal series solutions 272
B.3 The Olver growth theorems 282
Appendix C Recursion formulas for hypergeometric functions 289
Index of higher mathematical functions discussed 293
Bibliography 295
Index 307
.
--page0005--
Preface
The purpose of this book is to present applied mathematicians, numerical
analysts, engineers, physicists and computer scientists with an in-depth
study of that vast body of computational techniques based on the use of
recurrence relations.
These methods can be traced back to the dawn of mathematics. The
Babylonians used such a technique to compute the square root of a
positive number, and the Greeks to approximate it. Much later, Lagrange
A789) used a computational scheme based on a two-dimensional
nonlinear recurrence to compute an elliptic integral.
Current interest in the subject, though, can be attributed to some
almost offhand remarks made by J. C. P. Miller in an introduction to a
table of Bessel functions he had compiled for the British Association for
the Advancement of Science. This book, which was published in 1952,
coincided with the vaulting growth of large-scale digital computers;
although Miller's work was prefigured by some observations of Lord
Rayleigh A910) on the calculation of spherical Bessel functions, it was
only with the advent of large-scale computers that the great promise of
these algorithms could be brought to fruition.
What Miller perceived was that in a second-order linear recurrence
which has solutions sufficiently differentiated asymptotically, there is a
solution that may be uniquely characterized by one initial value and a
knowledge of its growth. This led to an algorithm for computing certain
solutions of the equation which required only a scant knowledge of their
pointwise values.
Very rapidly Miller's technique was generalized to a whole host of
computational problems. By now, virtually every important special
function of mathematical physics may be computed this way. In addition,
similar methods have been developed for computing the zeros of
functions, eigenvalues of certain differential operators, and the coefficients for
the expansion of functions in Taylor's series and in series of orthogonal
polynomials. Further, the linearity of the underlying recurrences leads to
very elegant statements about the growth of roundoff error and
truncation error in the algorithm.
The theory of nonlinear recurrences, of which Lagrange's scheme
furnishes an example, pursued a different path. Such algorithms are also
.
--page0006--
xii Preface
quite useful, for instance in the evaluation of infinite products, solutions
of functional equations, and a wide variety of definite integrals. In this
situation error analyses are more difficult, but the convergence is usually
extraordinarily rapid, much more so than for algorithms based on linear
recurrences.
Roughly the first three-quarters of the book is devoted to linear
algorithms, and the remainder to the nonlinear case. I have given
throughout a very large number of examples for, as I see it, one of the
primary purposes of this book is to put the reader interested in computer
software in touch with the many specific successful applications of the
algorithms. To this end, I have included an index of special functions
whose computation is discussed.
I have also tried to make the development self-contained. Thus the
book has substantial appendices incorporating the available material on
the asymptotic theory of linear difference equations and the construction
of recursion formulas for hypergeometric functions.
August 1983
Jet Wimp
List of symbols
z
n, k
Z°
Nr(a)
dS
x Ay
x vy
dp
dfcP
/:=g
X
lin [xl5 x2,
Q(p; n)
s(p; n)
r
integers: 0, ±1, ±2,...
always integers
integers: 0,1,2,...
Euclidean p-space
m1
positive (>0) reals
space of complex p-tuples
space of real sequences
space of complex sequences
{z e % | \\z - a\\ < r} (open disk)
boundary of S, S <= %
1, fc=0; 2, fc>0
Kronecker 5
minimum of x and y, x, y e S%
maximum of x and y, x, y e S%
degree of p, p a polynomial
coefficient of tk in p, p a polynomial
/ as denned by g
objects in a linear space
set spanned by xt,x2,... ,x(T
asymptotic form in Birkhoff series
asymptotic form in Birkhoff series
asymptotic form in Birkhoff series
sum in which first term is to be halved
finite sum in which first and last terms are to be halved
.
--page0007--
1 Introduction
Many of the most important techniques of numerical analysis developed
since the beginning of this century can be subsumed under the theory of
the vector system (or equation),
x(n + l)=f[x(n), n], A.1)
where x(n)e(t): for instance, that 4>(t) satisfies a differential equation. The
components of x(n) may then, for example, represent the coefficients for
the expansion of --page0008--
2 Introduction
known to be X. For linear equations this is never the case since the limit
B) always depends in a trivial algebraic manner on the initial or boundary
conditions. When f is nonlinear, however, X may depend in a
transcendental fashion on the initial condition even for the simplest equations.
The equation then serves as an algorithm to compute X. These kinds of
algorithms have been known for many years, the most famous one
bearing Gauss' name, though it is probably originally due to Lagrange.
However, these algorithms have not, apparently, hitherto received a
unified treatment in book form.
When f is nonlinear the equation A) is usually constructed so that X is
a root of an equation
is an invariant for the algorithm, i.e., a function with
= F[x(n)].
where F:%v
the property
F[x(n + l
These invariants sometimes take the form of integrals or solutions of
functional equations. The algorithms constructed this way are among the
most rapidly convergent in all numerical analysis. In fact in Chapter 14 I
construct an algorithm to compute a very general trigonometric integral
and show that the convergence is quadratic in every case.
Cases where f is nonlinear but where attention focuses on the behavior
of the trajectory rather than on the computation of X (which may not
even exist) have come to the forefront recently in the study of discrete
time dynamical systems, systems arising in multifarious disciplines: theory
of turbulence, population growth, biosystems. One may ask whether x(n)
has limit points or whether it is dense on some surface in %v or how it
changes as the initial point x is varied. The research going on in this area
nowadays can only be described as frenzied and I have the space to do
little more than offer the reader a tantalizing glimpse (in Chapter 12) of
some of the ideas currently being pursued. However the mathematicians
involved in this effort continue to be very productive and the reader after
completing the chapter will find himself with a useful list of names for
monitoring progress in the field.
In those cases where f is linear computing the trajectory itself is of
primary interest. The reader may at first glance fail to understand how
this could be a problem. Doesn't the equation offer an explicit format for
computing x(n)? The point is, of course, that x@) is often not known,
only asymptotic information about x(n), or perhaps the asymptotic
behavior of x(n) compared to other solutions of the equation. It is an
amazing fact, first articulated by J. C. P. Miller in 1952 (although the
origins of the idea can be traced back at least to Lord Rayleigh), that x(n)
can often be computed efficiently with no knowledge of x@) whatsoever
provided that ||x(n)||/|ly(n)||-»() as n—•* for any solution y(n) not a
Introduction 3
constant multiple of x(n). This statement, strictly speaking, applies only to
the homogeneous case, f[0, n] = 0. However, powerful techniques (due
primarily to F. W. J. Olver) based on the same principle have been
devised to handle the nonhomogeneous case also.
The most exciting work on the linear case has taken place in the last 20
years. Among the mathematicians productively involved are F. W. J.
Olver, C. W. Clenshaw, W. Gautschi, D. Lozier, J. R. Cash, J. Oliver and,
I immodestly add, myself. By now a lexicon of techniques has been
developed that can be used to do almost everything but peel apples, as
Norbert Weiner once said of the Fourier transform. With these techniques
we can now compute an array of important integrals and large classes of
higher transcendental functions, solve eigenvalue problems, find zeros of
functions, solve differential equations and determine the coefficients for
the expansion of commonly occurring functions in series of orthogonal
functions. It has been my purpose to give as many examples as I can of
the use of these algorithms throughout this book.
My approach is to start with the very simplest linear cases, first p = 1,
then p = 2. The p = 1 case, simple as it may be, nevertheless presents
many interesting features. The p = 2 case is of course more significant and
three full chapters are devoted to it. Although I generally disapprove of
mathematical exposition which moves from the specific to the general
there is, I think, in this situation a very good reason for it. The cases p = 1
and 2 are tractable and lend themselves to an elegant convergence and
error analysis because then the solutions of the system adjoint to A) can
be written in a simple way in terms of the solutions to (i) itself. When
p>2 this breaks down, and much of the rich theory characterizing the
simpler cases cannot be extended. Consequently the theory for p>2,
though equally important, is aesthetically less satisfying.
The nonlinear case is treated in the second half of the book.
The book contains three substantial appendices, all of which are quite
intrinsic to developments in the book and they should probably be
scanned before even Chapter 2 is attempted. Throughout the book
referents preceded by an 'A', 'B' or 'C refer to material in the
appendices, for example, A-X, Theorem A.3, (A.12), etc. In Appendix A I
have included a self-contained and virtually complete account of the
theory of linear difference equations. Appendix B contains an exposition
of the asymptotic theory of linear difference equations, first the classical
results such as the theorems of Poincare and Perron, next the very
important research of Birkhoff and Trjitzinsky on the expansion of
solutions in generalized asymptotic series and finally the recent work of
Olver, which provides growth theorems allowing one to connect given
solutions with their asymptotic forms. Appendix C contains a summary of
the contiguous function relationships for generalized hypergeometric
functions referred to earlier and which form the basis of many examples
given in the text.
.
--page0009--
2 General results on the
forward stability of recursion
relations
2.1 Background
In this chapter I shall discuss, computation based on the use of the
recurrence
B.1)
i/=0
in the forward direction. I will assume that the functions Av satisfy the
conditions given in Section A.I so that the equation possesses a unique
solution satisfying arbitrary initial conditions and such that the related
homogeneous equation has a fundamental set. Generally speaking the
problem is this: if w(n) is a solution of A) satisfying prescribed initial
data, what error can one expect to encounter when w is computed from
A) for subsequent n? This problem is different in nature from the
computational problems encountered in the more sophisticated Miller
algorithm to which much of this book is devoted but it is important
nonetheless. Many experienced numerical analysts are not fully aware of
the hazards posed by such computations. Consider the following very
simple first-order scalar equation provided by Gautschi A972a tr. 1973).
Example 2.1 Compute the numbers
w(n) = n![ex-en(x)],
where
fc=O
and x>0.
Note that w(n) satisfies the recurrence
l)y(n)-y(n
_ vn+l
B.2)
B.3)
We start with w@) = e\" -1. Taking x = 1 and rounding all computations
to 4 SF produces Table 2.1. With the seventh entry one realizes
something has gone awry, since, obviously, w(n) must be positive. Things get
1
Background 5
Table
2.1
n w(n) (approx.)
0
1
2
3
4
5
6
7
1.7180
0.7180
0.4360
0.3080
0.2320
0.1600
-0.0400
-1.2800
worse; in fact wA0) is computed to be —1023. This disastrous
accumulation of round-off error originated with the single rounding error in
w@). ¦
The situation encountered in this example can be also found when the
general difference equation i?[y(n)] = /(n) is used in the forward
direction. Whether or not small arithmetical errors introduced at one stage in
the computation will tend to grow with succeeding computation until they
overwhelm (in the sense of relative error) the true solution depends, as
we shall see, on the relative growth properties of the desired solution
versus the growth properties of other solutions of the original as well as
the related homogeneous equation i?[y(n)] = 0.
These are among the matters to be investigated in this chapter.
As is well known, the nonhomogeneous equation i?[y(n)] = f(n) can be
written as a system
y(n) + A(n)y(n + l) = f(n),
where A(n) is the crXcr matrix
A(n) =
A0(n)
0
0
Ao(n)
0
0
0 0 ¦¦¦ -1
y(n) = [y(n), y (n +1),..., y(n +
K\") = [/(\")/A0(n), 0, 0,..., 0]T
B.4)
B.5)
This matrix formulation can be quite useful for analyzing the stability
properties of the scalar equation. Also we will have occasion to consider
.
--page0010--
6 Forward stability of recursion relations
the matrix equation in itself, not necessarily arising from some given
scalar equation. A(n) then will be some given o-xo- matrix, nonsingular
for n5=0.
When A(n) is singular for certain values of n problems arise in the
existence and uniqueness of solutions which must often be resolved on an
ad hoc basis. Existence and uniqueness questions for difference equations
are discussed in great detail in the standard works, Milne-Thomson
A960), Meschkowski A959), and K. S. Miller A968), the latter being
particularly good on matrix systems.
2.2 Homogeneous systems
The analysis in this and subsequent sections draws heavily on the work of
Gautschi A967, 1972a tr. 1973).
We start with the system
y(n)+A(n)y(n + l) = 0,
B.6)
where A(n) is a 0, i.e.,
w(k): = w(k) + ||w(fc)|| e, ||e|| = 8. B.7)
Here || • || denotes any convenient vector norm. Also for any matrix B, ||B||
will be the matrix norm induced by || • ||.
Let ^e(fc) denote the set of all such vectors. We have, trivially,
SUp
B.8)
Now consider G) to be an initial condition. Then each vector w(fc) e
^e(fe) generates a 'perturbed' solution w(n) of F). By analogy with (8), I
wish to investigate the quantity
sup
t fixed, t+k, r>0.
Let
be a fundamental matrix of solutions of F); see Miller A968).
For every solution of F), in particular for *(n), one can find a constant
Homogeneous systems 7
vector c so that
w(n) = Y(n)c.
From the condition G) we find
(n) = Y(n) Y
e)
B.9)
so
sup
f.(fc) HwV«;il ll\"V«;il ||e||=e
The last term on the right can be computed in terms of the matrix norm
of Yir^Yiky1, so we get
Mk)\\\\
B.10)
This provides an indication of how the relative error at n compares with
the initial relative error at k.
Let
B.11)
Let Z: = {zA), zB),..., z(o°} be another fundamental matrix for F). Then it
is readily verified that for some constant matrix K
so
Z=YXK,
Z(n)Z-\\k)=Y(n)Y-\\k).
Thus the matrix on the right does not depend on the choice of the
solutions {y(h)(n)}, only on the system.
It is thus natural to adopt the following as a measure of stability for the
computation of w(n) by forward recursion from F).
Definition 2.1 The index of stability for the forward computation of w(n)
at k is
a(fc): = sup a(fc, n).
// a(fc)<°° F) is said to be stable for w(n) at the point k.
.
--page0011--
8 Forward stability of recursion relations
This property can be shown to be equivalent to one sometimes called
weak stability (see also Section 7.3) denned as follows:
Definition 2.2 // for each K > 0 there exists a constant C = C(K) > 0 such
that
sup a(k, n)= C<°°,
n>k
then the recurrence F) is said to be weakly stable for w(n). ¦
To establish the equivalence, 'stability at some kOweak stability' we
proceed as follows:
«(k + !,„)== l
\\\\Y(n)Y-\\k +
\\\\A-\\k)\\\\ \\\\A(k)\\\\ a{k, n) = Cond A(k)a(k, n).
1!! is called the condition number of B; see Isaacson
and Keller A966, p. 37).)
Similarly,
a(k,n) =
|w(n)
so
a{k + 1, n) =
1 a(k, n).
But
||w(k)H||A(k)||||W(k
so
\\\\Mk + DH.
-\\\\A(k)\\\\-\\
Combining the above results gives
[Cond Aik^aik, n)*Sa(k +1, n)=?[Cond A(k)]a(k, n),
and taking a sup over n gives
«(k+ 1)< [Cond A(k)]a(k),
which shows the equivalence (and much more).
Homogeneous systems 9
An important global property of stability is formulated in
Definition 2.3 If
sup a(k, n)= C<°°,
n,k»O
n>k
then the recurrence is said to be stable for w(n). ¦
Obviously stable ^weakly stable.
Definition 2.4 When neither Definition 1 nor 3 is satisfied, the term stable
is to be replaced by the term unstable in either situation. ¦
By the equivalence of vector norms in finite-dimensional spaces the
properties of stability and instability are norm-independent.
The definition of stability in Definition 2 is rather stringent. However,
even if a system is theoretically stable, if C is very large the desired
solution may be difficult to compute accurately. Conversely though the
system may be unstable for w(n) the function can usually be computed
accurately on small subsets of Z°.
Example 2.2 Second-order systems will be discussed in more detail later
on but it seems a good idea to offer an example to clarify these
definitions. The system
has a fundamental matrix
n + 1 1
We have ||yA) = 2n +1, ||yB)(w)||i = 2. A straightforward computation
shows that for w(n) = yA)(n)
Thus
a(k) = 2
Bk-
:k + i
hl)Bn-2k + l)
Bn + l)
ks=0.
The recursion is stable for the computation of yA)(n) at any point k and
also weakly stable, but not stable. ¦
.
--page0012--
10 Forward stability of recursion relations
It is easily shown that the system with fundamental matrix
is stable for ym(n) (and hence weakly stable and stable at any point)
provided lA.^lA.J.
The following provides a criterion for instability.
Theorem 2.1 (Gautschi) Let 0-3=2 and let the recursion relation F)
have a solution y*(n) such that
lim||w(n)||/||y*(n)|| = 0. B.12)
n—»°°
Then the recursion is unstable for w(n) at every point.
Conversely if F) is unstable for w(n) in the following manner:
lim a(k, n) = °° for some fixed k^O, B.13)
n—»oo
then a solution y*(n) of F) exists having the property A2).
Proof Note that for any nonsingular constant matrix B and any matrix
W(n)
lim ||W(n)B\\\\ =
because
B.14)
lim \\\\W(n)\\\\ = °°,
and
Now assume A3) holds. The solutions w(n) and y*(n) are linearly
independent. Therefore they may be completed to form a fundamental
matrix
Y(n): = [w(n), y*(n), yC)(n),. .., y(CT
Let's assume (without loss of generality) that we are working with the 11
norm. We have (see Isaacson and Keller A966, p. 9))
B.15)
a(k, n) = *> for every kss(). Thus the recursion is unstable
Now applying the result A4) to the matrices
W(n):= Y(n )/||w(n)||, B: = ||w(fc)|| Y(fc)-\\
shows Iim
for w(n).
Conversely assume A3) holds. Let
Nonhomogeneous systems 11
be a fundamental matrix for F). The relation A4), again applied to the
matrices A5), yields
Let Y(n):=[y(n
and the theorem follows. ¦
2.3 Nonhomogeneous systems
We now consider
y(n)+A(n)y(n + l) = f(n), ns=0.
1 will assume once again that Det A(n) f 0, and that f(n) #0.
B.16)
.
--page0013--
12 Forward stability of recursion relations
Let
where w is the solution to be computed. I will assume K is infinite.
The general solution of A6) is
y(n) = Y(n)c+w(n),
where Y(n) is a fundamental matrix of solutions of the homogeneous
equation. As before, if w(n) is a solution of A6) whose initial vector w(k)
approximates w(k) to within e, in accordance with G), one finds that the
equation (9),
w(n)=w(n) + ||w(k)|| Y(n)Y(k)~1e,
holds also for the nonhomogeneous system. Thus equation A0) holds and
Definition 1 of the previous section still makes sense with one minor
change, that is, n, keK.
Definition 2.5 The recurrence A6) is said to be stable for w(n) if there
exists a constant C>0 such that
sup a(k, n) = C<°°,
n>Tc
n,keK
a(k, n) as in A1). Otherwise A6) is said to be unstable for w(n). ¦
The following result follows immediately from Theorem 1.
Theorem 2.3 Leto-^2 and let the associated homogeneous recursion A6)
have a solution y*(n) such that
Then the nonhomogeneous system A6) is unstable for w(n) at any
point. ¦
2.4 The first-order scalar case: forward vs. backward recursion
For o-=l the difference equation F) is
y(n) + a(n)y(n + i)=f(n), n&0.
We find
|w(k)y*(n)| p(n)
a(k, n) =
w(n)y*(k)
p(k)'
B.17)
p(n):=|w@)y*(n)/w(n)|.
First-order scalar case 13
where
is a solution of the homogeneous equation. (Empty products are to be
interpreted as 1.) Going from k to n we find the error is increased if
p(n)^p(k) and decreased if p(n)0.
B.19)
-• B.20)
Since a solution of the homogeneous equation is y(n) = l, p(n)= 1 and
a(k, n) = l. Thus B0) is stable for the computation of 0, n>0.
B.21)
.
--page0014--
14 Forward stability of recursion relations
w(n) = En+1 satisfies the equation
^y(n) + y(n + l) = ^.
p(n) is found to be
B.22)
It can be shown that if x =?0.61006 . .., then p(n) decreases
monotonically from 1 to 0 as n increases so that B2) is stable for En. If
x>0.61006 ..., then p(n) first increases to a maximum at n = n0 and
then tends monotonically to 0 so again B2) is stable for En. If the error
amplification p(n0) is acceptable B2) may be safely applied to the
computation of En{x) for all n. ¦
2.5 The computation of successive derivatives
This topic provides a very interesting application of the material in the
preceding sections. Sources for this material are Gautschi A966, 1972a)
and Gautschi and Klein A970).
Let f(z) be analytic on some set A; it is desired to compute the
quantities
zgA-{°>- B.23)
From the identity
dznJ \"dz\" Vz \"I dzr
we see that (n) depends
only on the value of / at 0 provided / is analytic in a large enough disk
containing 0 and z.
Theorem 2.4 Let /(?) be analytic in the disk
NT(z): = {t\\\\l-z\\*r), r>\\z\\,
and assume /(?)/? is not a polynomial. Then B4) is stable for \\z\\ so
z
w: =
?edNr.
Thus
lim
z4>{n)
y*(n)
Also
B.25)
B.26)
4=: Since f(C)/? is not a polynomial the set K = {n | <^>(n) j= 0} must be
infinite. Suppose /@) ^ 0. Because of B5) there are constants C, C and
n()>0 such that
0(n)
Thus
y*(n)
P(n)
P(k)
n>n0.
zcf>(k) /zcf>(n)
y*(k)/ y*(n)
c
for k, n sufficiently large, and this shows stability.
=>: By contradiction. Assume /@) = 0. Then lim,,^^ p(n) = oo; So by
A8) the equation is unstable for w(n). ¦
Example 2.6 (Gautschi and Klein) For fixed x ^ 0 compute
.
--page0015--
16 Forward stability of recursion relations
The recurrence relation is
, ,xv(n + l) Imete)
y(n)+i^nr=~(^+ir-
Using the method described in A-IX this equation may be solved and we
find
en as in equation B). Thus
sin x
p(n) =
Since e^e^-ix) -»lasn-»°° this confirms the instability asserted by the
theorem. It might be thought that (n) could be computed by backward
recursion starting, say, with 0 such that
n>n0.
Although these estimates are not uniform in x they provide a good
indication of what actually happens. For a given x the curves of p(n) start
out decreasing, attaining a minimum value near n = x, and then begin to
grow essentially monotonically. Although the recursion is unstable for
(n), it is stable in the region O^n ^|x|. Clearly, however, the recursion
cannot be used in an unrestricted fashion in either the forward or
backward direction. ¦
The derivatives (n) of the previous example arise in the theory of
cardinal interpolation. Let f:R—*C and a real number h >0 be given.
The series
zeC.
Scalar equations of higher order 17
is called the cardinal series off (with respect to h). Obviously, Fh(mh) =
f(mh), meZ. Further, it is known that if / is an analytic function
bounded in the strip R x [-ia, ia] and
f°°
|/(T±ia)|2dT<°°,
then the series converges uniformly on compact subsets of the strip and,
further, limh_^0 Fh(z) = f(z); see Kress A972). For some functions the
representation is exact, e.g., when / is entire with
|/(z)|=?ceplIm21,
for some cs=0, 0=?p<7r/h, then Fh(z)=f(z) with uniform convergence
on compact subsets of C. Often the derivatives of Fh(z) are required,
which means one needs the quantities
_d^_ . h^
dz\"Sm (z-mh)
d\" sinx
dxn
:=?(z-mh).
Gautschi A972a) has observed that this series arises in communication
theory and, for h = 1, has to do with the reconstruction of a function with
a limited frequency band when its sample values at equidistant points are
known.
2.6 Scalar equations of higher order; minimal and dominant
solutions
The stability of the nonhomogeneous scalar equation i?[y(n)] = /(n) is
rather easy to appraise. Basically this is due to the fact that in the scalar
case the fundamental matrix for the related homogeneous system assumes
a simple form, i.e.,
V(n) =
yB)(n)
yB)(n
y(--page0016--
18 Forward stability of recursion relations
where w(n) is the desired solution of the nonhomogeneous equation A).
Then i?[y(n)] = /(n) is unstable for w(n) at every point.
Proof Let y*(n): = [y*(n), y*(n + 1),..., y*(n + (or is dominated by «/>). ¦
Definition 2.7 A nontrivial solution (x(n) of the homogeneous difference
equation i?[y(n)] = O is called a minimal solution if /x(n) is dominated by
each of the other solutions in some fundamental set of solutions containing
Scalar equations of higher order 19
(This means, of course, that linv^oo fx(n)/y(n) = O for any other
solution y(n) linearly independent of (x(n).)
Definition 2.8 A nontrivial solution l).
A corollary to the previous theorem can now be stated in terms of
these definitions.
Corollary Let (x(n) be a minimal solution of i?[y(n)] = O. Then this
recurrence is unstable for (x(n). I
Example 2.7 Consider the difference equation of Example B.3:
xg(n) + (n + 2)g(n + l)-2g(n + 3) = 0, x>0.
As shown there, the equation has a fundamental set {y(h)(n)} with the
property
n ' n I
n n
Thus the equation possesses a minimal solution y°\\n).
Example 2.8 The equation
l)(n+2)y(n+2) = 0, -7r--page0017--
20 Forward stability of recursion relations
possesses no minimal solution, for a fundamental set is yA)(n)=l,
y(
yB)(n) =
There seem to be no simple necessary and sufficient conditions which
will guarantee i?[y(n)] = 0 has a minimal solution. However, in many
cases the asymptotic theory of linear difference equations developed in
Appendix B can be used to determine the existence of such solutions. Of
particular importance are the Theorems B.I and B.2.
Demonstrating stability is generally more difficult than demonstrating
instability and usually requires more asymptotic information than
furnished by either of the preceding theorems. The construction of the
matrix Y(k) and the estimation of it can be a formidable problem,
even in the scalar case. I give an example to show that in some cases the
computations may actually be performed.
Example 2.10 (Repeated integrals of the error function) The functions
m(n): = i\" erf c x : = -?- f ^\"f* e~'2 df, x > 0, B.28)
V77 Jx n!
'2df, x>0, B.29)
V V\"\"~V,rJL.
satisfy the difference equation
y(n)-2xy(n
+ 2) = 0. B.30)
Neither of the preceding theorems suffices to indicate whether this
equation has a minimal solution. However, it is known that
The result E.8) shows
Tf
-V©-
Gautschi also gives the lead term of this result A961b, p. 230). By the
discussion in Section B.2 it follows that there is a second solution of the
equation y*(n) with the property
L n'n J
(n) may be identified as a constant (independent of n, that is)
Scalar equations of higher order 21
multiple of y*(n) since (-l)>B)(n; -x) = Owe have
A direct computation shows
B) I\" r(n) 1
LBn) r(n) ~Bn)
where
r(\") • = i B1/ \\ :
1)ne-2xVBn))_ q^
Let us first examine stability of the recurrence for the computation of
a\\n). A little algebra yields
in the 1-norm. Thus by C1) the recurrence C0) is stable for n)(n) is a minimal solution, the equation is unstable at any point
for --page0018--
22 Forward stability of recursion relations
|| • \\l norm
I
for some C>0.
(ii) If
lim
< 1, lim
yB)(n
yB)(n)
B.32)
then the equation is stable for the computation of yA)(n).
Proof (i) We have
a(k, n) *? (sup jSpl) \\\\Y(k)-l\\\\ ||yA)(k)||,
and the result follows easily. (Note D(k) $¦ 0.)
(ii) Choose N so that
., yB)(n +
Then
yB)(n)
: = |yA)(k)yB)(k +1) - y B)(k)y A)(k + DI
yA)(k)yB)(k + l
n>N.
1-
k>N.
Thus we have
¦ A),
|y U)(n +
i ,l
l
0-yB)(n + :
max [|y°
<4C, n>k>N.
Since D(kOt0, a(k, n) is finite on the grid 0=?fc--page0019--
24 First-order equations used in the backward direction
Table 3.1
I
N
12
16
0
1
2
3
4
2.327 386
-0.758 774
0.832 854
-0.587 601
0
2.350 395
-0.735 767
0.878 869
-0.449 555
0.552 182
2.350 402
-0.735 759
0.878 885
-0.449 507
0.552 373
2.350 402
-0.735 759
0.878 885
-0.449 507
0.552 373
To fix ideas, let x = 1 and take N = 4,18,12,16. We obtain Table 3.1 for
the values of yN(n), 0=?n=?4. ¦
Obviously some kind of convergence is occurring in the example. In
fact we can show easily that limN_^ooyN(n) = w(n). By the theory in
Appendix A I can write
Using the initial value yN(N) = 0 and solving for c yields
w(N)y*(n)
yN(n)=w(n) —
Since
I have
2ex
|yN(n)-w(n)
y*(N)
x>0,
,2e\"n!xN~\"
(N+l)! '
C.5)
and the assertion follows. Of course the accuracy deteriorates with
increasing n but, nevertheless, for n e any fixed subset of Z° and for any
fixed x, w(n) may be computed as accurately as one wishes by this
method, at least theoretically. There is as yet no assurance that some
potentially destructive build-up of round-off error caused by finiteness of
the computer representation of numbers will not occur. I will return to
this topic later.
This kind of an algorithm—by now actually a class of algorithms
applicable to many different kinds of difference equations under many
different conditions—is called a Miller algorithm after J. C. P. Miller who
used a variant of the algorithm in 1952 to construct tables of the Bessel
Convergence and error analysis 25
functions /n(x). Apparently, though, the idea is much older. It can be
traced back at least as far as a 1910 paper of Lord Rayleigh.
The algorithm of the previous example displays two features which
characterize the general class of Miller algorithms:
(i) the difference equation is used in the backward direction;
(ii) no initial (i.e., for n =0) values were required.
The second feature, in the general Miller algorithm, takes the form that
only global information, if any, is required to single out the desired
solution of the equation. Such information may be very general; for
instance, it may only involve knowing the sum of some series of the form
Z c(k)w(k). This is"one of the most attractive features of the algorithm,
since in almost any practical situation the necessary information is
available.
The following notation will be useful.
EN(n):=wN(n)-w(n),
denotes the error of the algorithm and
C.6)
C.7)
the relative error.
3.2 Convergence and error analysis
Consider the general first-order nonhomogeneous equation
(n)y(n + l) = b(n), ns=0. C.8)
We assume as before that a(n)^0. Suppose one wishes to compute a
solution w(n) of (8) by the use of Miller's algorithm, i.e., one sets
yN(N) = 0 and computes yN(n) for n = N-1, N-2,..., 0 from
yN(n) = b(n) - a(n)yN(n +1).
C.9)
There are, of course, two apparently distinct problems associated with
this computational scheme. First, does the algorithm, done in infinite
precision at each step, converge? And, second, how does error in the
algorithm accumulate assuming that at each step random errors are
introduced into the computations?
I shall show that for the simple first-order equation both of these
problems can be treated at the same time. The analysis is essentially that
of Gautschi A961a).
Let {eN(n)} be an infinite matrix of errors. What one actually computes
.
--page0020--
26 First-order equations used in the backward direction
when one applies the Miller algorithm to (8) is a quantity which satisfies
zN(n) = [b(n) — a(n)zN(n +1)~\\[1 + eN(n)], n = N—l, N—2,
C.10)
with zN(N)=eN(N). We rewrite this as
zN(n) + a(n)(l + eN(n))zN(n + 1) = b(n)(l + eN(n)),
and (8) as
w(n) + a(n)(l + eN(n))w(n +1) = b(n) + a(n)eN(n)w(n +1).
Subtracting these two equations and letting
EN(n):=zN(n)-w(n)
gives
EN(n) + a(n)(l + eN(n))EN(n + 1) = eN(n)w(n). C.11)
Solving this equation by the procedure in A-IX gives
C.12)
where y?(n) is any nontrivial solution of the homogeneous equation
associated with A1). Further, (A.24) shows y?(n) may be written
where y*(n) is any nontrivial solution of the homogeneous equation
associated with (8) and e%(n) is an appropriately chosen nontrivial
solution of
y(n)-(l + eN(n))y(n + l) = 0.
We thus obtain
C.13)
C.14)
Equation A4) with eN(n) = 0 produces the following convergence
criterion:
Theorem 3.1 Let yN(n) be computed, as described, from (9) in exact
arithmetic. Then
w(N)
Urn yN(n) = y(n) O^^=o(l),
for every nontrivial solution y*(n) of the homogeneous equation
y(n) + a(n)y(n + 1) = 0.
Convergence and error analysis 27
We then have
|yN(n)-w(n)| = |y*(n)|
w(N)
y*(N)
Suppose a substitution Y(n) = \\(n)y(n) is made in the original
difference equation. Let Y(n) be computed from this new equation by an
application of the Miller algorithm and then y(n) computed from Y(n).
Can \\(n) be chosen to substantially reduce the error [YN(n)/A.(n)]-y(n)?
Gautschi (ibid) has observed that the computations cannot be improved in
such a facile way. Indeed, this is obvious from the representation
A description of the error accumulation follows from A4).
Theorem 3.2 Let w(n)^0. Then
zN(n)-w(n)
w(n)
y*(\")
w(n)
|y*(N)|
411-
y*(j)
C.15)
where
e = sup |eN(n)|.
Proof The theorem results by taking as a solution of equation A3)
Example 3.2 Consider Example 1. Using the value given by D) for
y*(n) and the inequality E) yields
zN(n)-w(n)
w(n)
n\\
xnw(n)
i N!
+ 2eexxn X -^
Using the fact that (j + n +1)! ss/! (n + 1)! produces
zN(n)-w(n)
w(n)
N\\w(n)
2ex \\ 2exne
nea+E)x
.
--page0021--
28 First-order equations used in the backward direction
The latter term represents the contribution to the error produced purely
by accumulated round-off error. Since it is independent of N it cannot be
diminished by taking N large. ¦
Example 3.3 (General moment integrals) Let
f
n):= tnilf(t)dt,
o, neZ°,
where i/zeC'fa, b] and satisfies the linear differential equation
tp'(t) + d\\J/(t) = g(t), d = const.
A direct computation shows that (n) satisfies
d
y(n)-
(n + 1)
y(n + l) = c(n),
C.16)
with
c(n): =
We have
(n + 1)
n!
\\ = max(|a|, |b|), for some M. Thus Miller's algorithm for the
computation of @)y*(n)
l--page0022--
30 Second-order homogeneous equations
Table 4.1 Computation of In(l).
wN(n)
0
1
2
3
4
1.265 927 978
0.565 096 950
0.135 734 072
0.022 160 665
0.002 770 083
1.266 066 460
0.565 159 364
0.135 747 732
0.022 168 435
0.002 737 123
1.266 065 876
0.565 159 103
0.135 747 669
0.022 168 425
0.002737120
by
wN(n) =
w@)yN(n)
yN@)
The algorithm so defined is called the simplified Miller algorithm.
Example 4.1 Calculate In(x) from the recurrence
2(n-
D.7)
y(n)-
y(n
and the normalization
1= I ek(~l)kI2k(x).
D.8)
D.9)
(Here ek = 2 — Sk0; see List of Symbols.)
Some sample computations are given in Table 4.1 for x = 1. The entries
in the last column are accurate to the figures given. This i6 essentially the
algorithm given in Miller A952). Another possible normalization relation
for this example is
e* =
D.10)
lc=0
As one might suspect, the success of the method depends on the rate of
growth of the desired solution w(n) compared with other solutions of the
equation. Recall that for any homogeneous linear difference equation the
set M of minimal solutions is a one-dimensional subspace of the solution
set if of the equation. (Any minimal solution ym(n) may be completed to
a basis {yA)(n), yB)(n),..., y(cr)(n)} of Sf. For the case a = 2 the actual
construction of yB)(n) may be carried out by the formula A-X.) Further,
no matter what the order of the equation is a minimal solution is uniquely
determined by one initial value. These observations suggest that minimal
solutions are ideal candidates for the application of the Miller algorithm.
I
while the behavior of a linearly independent solution is
The Miller algorithm 31
If the equation has a minimal solution but w(n) is not a constant
multiple of that solution then forward recurrence (which requires two
initial values) should be used to compute w(n). When the equation
possesses no minimal solution either backward or forward
recurrence may be unstable. However, there are special techniques which
depend on averaging solutions which can then be used (see Section 4.7).
The reason the algorithm works so well on In is that
D.11)
D.12)
Thus In is strongly minimal. A good question is, how may Kn be
computed (without the use of initial values)? This question is addressed in
Sections 5.2 (Example 5.2), 5.7, 7.3.
If only one value of {/„} is required, the Miller algorithm loses some of
its attractiveness. Henrici A979) has shown that in this case it may be
preferable to use the fast Fourier transform on an appropriate generating
function. At any rate, often many of the values {In} are required, for
instance, when it is necessary to evaluate an infinite series involving the
functions.
Before I discuss the problems of convergence and error propagation I
will give a formalization of the Miller algorithm due to Shintani A965)
which is very useful computationally since it allows one to increase N in
yN(r), r fixed, without recomputing yN(N+l), yN(N),..., yN@).
Theorem 4.1 (Duality theorem) yN(n) satisfies
b(N)yN(n) + a(N + l)yN+1(n) + yN+2(n) = 0,
with the initial conditions
Vn(n)=l, yn+1(n) = -a(n),
and
SN(n):= ? c(fc)yN(fc),
lc = n
satisfies
b(N)SN(n) + a(N + 1)SN+
,..., D.13)
D.14)
D.15)
^in) = c(N+2),
0=?N=n,n + l,n+2,..., D.16)
.
--page0023--
32 Second-order homogeneous equations
with initial conditions
Sn+1(n) = -a(n)c(n) + c(n + l).
D.17)
Proof See Theorem 7.3.
Letting n = 0 gives
Corollary SN satisfies
So = c@), S, = c(l)-c@)a@).
D.18)
The two previous theorems' results are particularly welcome when only
a few values of the sequence {w(n)} are required. For these values of n
yN(n) can be computed for increasing values of N, N = 0,1, 2,..., from
A3) and A4) while from (IB) SN can be computed without the necessity
of computing the sum E) at each step.
The convergence properties of this algorithm are described in the two
following theorems.
Theorem 4.2 Let A) have a minimal solution w(n) for which
S:= ? c(fc)w.
lc = 0
Then the Miller algorithm C)-E) converges to w(n) O
UmTNS? = 0, D.19)
where
D.20)
y(n) any dominant solution of the equation.
Proof Let
UN:=
c(fc)w(fc).
A direct computation using D) and E) gives
, s -STNy(n)+w(n)(UN + TNS%) ,.,.,
wN(n)-w(n)= — r^ . D.21)
S — UN — TjsjOn
If the limit A9) exists, then wN(n) is defined for n sufficiently large
and wN(n) —» w(n). The converse is obvious. ¦
The Miller algorithm 33
The following corollary is useful in many applications.
Corollary Let A) have a minimal solution w(n) with w@)^=0.
Then the simplified Miller algorithm defined by C), G) converges to
w(n).
Conversely, when the algorithm converges in the sense that
hm^yN(n)/yN@) = w(n)
exists, then w(n) is a minimal solution of the equation.
Proof The proof of the first part of the corollary is direct.
To show the second part, choose a fundamental set {y(h)(n)} with the
properties
(i) yB)(n) does not vanish for n>n();
(ii) r(n):= yA)(n)/yB)(n) does not approach a finite nonzero constant C.
By A-IV such a yB)(n) exists. Also, if r(n)-»C, choose instead the set
{yA)(n)-CyA)(n), yB)(n)} which will be linearly independent by
Thomson A960, p. 360).
I have
yN@) yA)@) -
Either r(n) -* ±°° or r(n) —* 0 if convergence is to occur. In the first
case convergence is to yB)(n)/yB)@) which, indeed, is a minimal solution.
In the second case we must have yA)@)^=0, else yA)(n) = 0 and
{y(h)(n)} is not a fundamental set. So convergence is to yA)(n), which is
the minimal solution. ¦
The convergence of the Miller algorithm also depends intimately on the
speed of convergence of the normalization series ?c(fc)w(fc), or, more
accurately, on the decay rate of the coefficients c(fc). In fact if w(n) is a
minimal and y(n) a dominant solution in most instances the relative error
= wN(n)-w(n)
w(n)
is O(w(N+l)c(N)y(N)/y(JV+l)) for the Miller algorithm and
O(w(N+ \\)ly(N+ 1)) for the simplified algorithm, as can be inferred,
under appropriate conditions, from B0), B1). The latter is usually much
smaller. Of course the price to be paid in the simplified algorithm is
that w@) must be known.
When the equation is of Birkhoff type (i.e., a(n), b(n) have appropriate
asymptotic series, see Section B.2) a- much stronger statement about
.
--page0024--
34 Second-order homogeneous equations
convergence can be made. The theory of linear difference equations
discussed in Section B.2 shows that in such a case there is a fundamental
set of solutions having the behavior
y(h)(n)~eQ*(p;n)Sh(p,n), n-^°°.
(See Appendix B for the notation in this equation.)
Wimp A969) demonstrates the following.
Theorem 4.3 Let equation A) be of the Birkhoff type, let w(n) be a
minimal solution of the equation, let
c(fc) = O(eQ(k)fca),
Q(k) of the form (B.ll) and let
eQ(k)fca+1lnfcw(fc) = o(l).
Then the computation of w(n») by backward recursion based on C)-E)
converges.
Proof See Wimp A969). ¦
Sources for much of the material to follow, as well as good general
surveys of many computational aspects of three-term recurrences, are in
Gautschi A967, 1972a).
When the simplified algorithm converges an interesting property of a
related continued fraction results.
Theorem 4.4 Let the limit
1
lim
exist.
Then
yN@)
:=w(n)
-1
b(n) b(n
w(n + l]
w(n) a(n)-a(n + l)~ a(n
D.22)
Proof Note that for each n, yN, N sufficiently large. In fact if
Yn@) = 0 for an infinite number of values of N, then w(n) is not defined;
clearly yN(n) ^ 0, n >0, for N sufficiently large.
Now let rN(n):= yN(n + l)/yN(n). From C) I find for N sufficiently
large, n fixed,
-1
D.23)
The Miller algorithm 35
or
-1
b(n) b(n +
b(N-2)
a(n)- a(n + l)- a
D34)
by repeated application of B3), since tn(N) = 0. (Note a (N-1)^0 or
else yN(n) = 0, 0=sn=?iV-l.) From the above expression, the theorem
follows on letting N —* °°. ¦
A criterion for the existence of a minimal solution can be based on the
continued fraction B2).
Theorem 4.5 (Pincherle A894)) The continued fraction B2) with n = 0
converges A) has a minimal solution jx(n) with ^@)^0.
When the fraction does converge the relationship B2) holds (with
w = fx) for each n for which ju,(n) ^ 0.
Proof ^: The numerator approximants P(n) and the denominator
approximants Q(n) of the continued fraction both satisfy
(n) = a(n - l)y(n - l)-b(n-2)y(n -2),
with
P@) = 0,
Q@)=l,
= a@).
D35)
D.26)
Let {y(h)(«)} be a fundamental set for A). The equation B5) is essentially
the adjoint equation and has solutions
(-ir+h+1y(h)(n + l)/D(n);
see A-VI. Thus the fact that the ratio of two solutions of B5) approaches
a limit implies the ratio of two solutions of A) approaches a limit; again,
call these solutions {y(H)(n)}. Thus
lim
: = c.
Define
H(n):=ym(n)-cya\\n).
The initial conditions B6) guarantee that {y(h)(n)} are linearly
independent. Let y(n) = c1ya)(n) + c2yB)(«) be any other solution of A). Then
= 0,
y(n) n™c1{[yu\\n)/yi2\\n)l-c}+c2+c1c
unless c2 + ctc = 0. But if this is true, y(n) is a multiple of fx(n).
4=: Obvious. ¦
.
--page0025--
36 Second-order homogeneous equations
When the coefficients a(n), b(n) of the equation A) are negative the
continued fraction B2) provides a very clean criterion for convergence of
yN(n)/yN@).
Theorem 4.6 Let a(n), b(n)<0, n 3=0. Then the limit
1
lim^,
N— yN@)
exists for « = 1 o it exists for all n 3= 1 and this happens the series
v a(k)a(k +
Hk)
D.27)
diverges.
Proof By Pringsheim's theorem (Khovanskii 1963, p. 45) the continued
fraction B2) converges •?> theories
I
k=2
k-2)a(n +
b(n + k-2)
1/2
is divergent.
But this series diverges with B7). Thus lim^,,^ yN(n + l)/yN(n) exists
for one neZ°Oit exists for all neZ° and this happensB7) diverges.
(Note that yN(n)^0, 0=?n=?N.) I have
7T7n\\ = 11
If any yN0 + l)/yN0) converges they all do and the quantity on the left
converges. ¦
This theorem does not reveal whether the simplified Miller algorithm
converges to the desired solution. The following is more or less a negative
criterion.
Theorem 4.7 Let
a(n),b(n)<0,
n3=0.
If the simplified Miller algorithm converges it converges to a solution
w(n) of the recurrence having the property
\\a(n)\\
b(n)
a(n + l
w(n) ~~\\a(n)\\'
D.28)
Proof This follows from the representation B2) and the inequality
(Khovanskii A963, p. 7 A.10))). ¦
Reduction of error by use of asymptotic information 37
Example 4.2 For the recursion (8) satisfied by In(x), a(n) = -2(n + l)/x,
b(n) = —1, and the series B7) diverges. From the asymptotic estimates
A1)-A2) we conclude immediately that the algorithm converges to
() m
Example 4.3 Consider the recurrence B.30) satisfied by w(n) =
inerfcx; here a(n) = -2x, b(n) = -2(n + 2). The series B7) diverges so
the algorithm converges. However, unless we have some a priori
knowledge of the asymptotic behavior of w(n) there is no easy way to decide to
which solution the Miller algorithm converges. B
As the previous example shows, it is difficult, lacking any additional
information, to characterize the limit function of the Miller algorithm by a
cursory examination of the difference equation. Even when a great deal
of information is available about the coefficients a(n) and b(n)—for
example, the fact that they are rational functions of n—it may be very
difficult to identify w(n).
The Olver growth theorems given in Section B.3 sometimes allow one
to decide whether a solution having a prescribed initial value is a minimal
solution and thus to conclude that the Miller algorithm converges to that
solution. In effect the theorems allow us to match the given solution with
its asymptotic expansion. Example B.4 demonstrates how the theorems
can be used.
4.2 Reduction of the error by the use of asymptotic information
Gautschi A967) and Scraton A972) have observed that when asymptotic
information is available about the desired solution of the difference
equation this information may often be put to advantage.
Suppose
where --page0026--
38 Second-order homogeneous equations
Table 4.2 Error in the tabulation of w(n) by Miller's algorithm (upper part)
and the revised algorithm (lower part).
10
15
Miller N 20
25
30
10
15
Revised N 20
25
30
n = l
-38 060
-5 427
-1035
-239
-63
-206
-18
-2
0
0
n = 2
-76 120
-10 854
-2 070
-477
-126
-412
-36
-5
-1
0
Error xlO8
n = 3
-123 695
-17 638
-3 364
-776
-205
-669
-58
-8
-1
0
n = 4
-183 957
-26230
-5 002
-1153
-305
-995
-87
-12
-2
0
n = 5
-259 680
-37 028
-7 062
-1628
-430
-1404
-123
-17
-3
-1
satisfies
(n+?)y(n)-2(n +
and has the behavior
w(n) = Ke\"VBn)[l
see Section 5.2; w(n) is minimal. Since
w(n)/w(n +1) = l+Bn)~m + Dnyl + O)
D.30)
(as is easily determined by direct substitution of a series of the form B.13
with p = 2 into the equation) we take
yN(N+l):=l,
Table 4.2 is self-explanatory.
4.3 Error analysis, the simplified algorithm: case of negative
coefficients
The truncation error of the algorithm is the error in exact arithmetic, i.e.,
|w(n)-wN(n)|. When the difference equation has real negative coefficients
the error study of the simplified Miller algorithm is very elegant. The
analysis and numerical data to follow are due to Tait A967). Note that if
a(n)>0 while b(n)<0 then the substitution y(n) = (-l)nz(n) brings the
equation into the desired form. Most of the second-order homogeneous
Error analysis, the simplified algorithm 39
difference equations that arise in practical applications are amenable to
this analysis.
I shall actually consider a slight generalization of the algorithm defined
by G), i.e., let yN(n) be defined by
n = N,
andi?[yN(n)] = 0, n = JV-1, N-2,.. .,0. Allowing yN(N+l) to be other
than 0 can give rise to certain computational advantages.
Again, let w(n) be the solution to be computed. Define, as usual,
EN(n) := wN(n) ~ w{n),
wN(n) as in G).
D.31)
Theorem 4.8 Let {y&V)}, h = l,2,bea fundamental set for if[y
with the property
n=N+l;
l, n = JV+l;
Assume yjJ(n)^0, 0=?n=?N, and wN@)^0. Then
„ , , Gw(N)-w(N+l))
yN@)
._ b(n-k-l)b(n-k)-- -b(N-l)
lc=0
where
Proof We have
w(n) = w(N)yW(n) + w(N+ l)yg>(n),
D.32)
D.33)
so by direct substitution
However, by (A. 12),
y(s(n) y%\\n -1) b{n - l)ft(n) • • • b(N-1)
so the result follows by induction.
.
--page0027--
40 Second-order homogeneous equations
Corollary
@
(ii)
T
1 ^
1 t
<1O
>\\b)\\,
Gw(N)-w(N+l))y»>@)
0=?fc=?n-2.
I now make the following assumptions:
(A) w(n) is positive nonincreasing;
(B) a(n), b(n)<0, O^
Note that (B) guarantees that y^(n), yg)(n)>0. The use of equation^
A) with n replaced by (n-fc-2) shows that Corollary (ii) holds.
Obviously Tk alternates in sign and is decreasing in magnitude.
For our error studies it is sufficient to take 7 = 0 for the formula C2)
displays explicitly how the error depends on 7. Thus y^M
Now EN@) = 0,
and from the difference equation
EnB)=
M
By induction the following may be established.
Theorem 4.9 Under the stated conditions
J=n-i
(ii)
D.34)
D.35)
w(N)yN(*)yN(n-l) '
Remark Note that (i) implies that EN(n) alternates in sign (with n).
Proof Both yJJXn) are nondecreasing as n decreases from IVtom-1.
Error analysis, the simplified algorithm 41
Assume
Then
D.36)
J=n-1
yN(n-DyN(n-2) = -a(n-2)yN(n-lJ-b(n-2)yN(n)yN(n-l)
s* -a(n-2)yN(n - l)[-a(n - l)yN(n)]
-b(n-2)yN(n)yN(n-l)
= [a(n -2)a(n - 1) - b(n-2)]yN(n)yN(n -1),
and this provides the statement C6) with n —*¦ n — 1. Also the statement is
true for n=N— 1. Thus the induction is completed to demonstrate C6)
for m-l«n«N-l. Now w(N)s*-a(N)w(N+l). Using this and C6)
in Theorem 9 (ii) gives the theorem. ¦
Corollary Let the conditions of Theorem 9 hold and a(n)=?-l for
n 3= m - 1 for some m 3= 1 and let b(N) < 0. Then
a(N)a(N-
-1 [a(j)a(j + l) - b(j)Y
m-l=?n=?N-2. ¦ D.37)
The usefulness of equation C7) results from the fact that the quantity
on the right depends only on the known coefficients of the equation.
Example 4.5 (Bessel functions) Consider the computation of In(x)
from equation (8). We get
x2/4N2
D.38)
Denote the denominator of the right-hand side by V(N).
A simple computation shows V(N) satisfies
1
Thus convergence is
)
°((tOn(niJ)-
The data in Table 4.3 illustrate a typical computation.
There is one more useful result concerning EN(n),
.
--page0028--
42 Second-order homogeneous equations
Table 4.3 Computation of In(x); x = 2; N= 10.
win)
bound |E3(n)|
|E5(n)l
11
10
9
8
7
6
0
1
10
91
738
5 257
3.043 890744 (-6) 3.044185 903 (-6)
2.769 940 577 (-5) 2.769 936951 (-5)
2.246 391369 (-4) 2.246 391420 (-4)
1.600173 364 (-3) 1.600173 364 (-3)
1.099 (-4)
1.506 (-6)
2.641 (-8)
6.142 (-10)
0.9695 (-4)
1.309 (-6)
2.207 (-8)
0 7 489 051 2.279 585 302
2.279 585 302
Theorem 4.10 ^
(i) Let (A), (B) hold and a(n)^b(n). Then
(ii) Let (A), (B) hold. Then
|Ej(n-l)|>?|E5(n + l)|, l=
Proof Left to the reader. ¦
When a(n), b(n) are not necessarily negative then, as might be
expected, the analysis is not so simple nor the convergence so rapid. The
following result is, however, useful. Again we take 7 = 0 and assume
0. Let
a: = sup
Theorem 4.11 Let
-a(fc)>l + |b(fc)|,
and
-a(fc) >max [1 + \\b(k)\\, 2 \\a(k)\\],
Then
D.39)
Proof See Tait, ibid.
The round-oflF error can be handled similarly. It is not so easy to discuss
the relative round-oflF error as was done for the first-order equation, cf.
Error analysis, the simplified algorithm 43
Theorem 3.2. I shall instead adopt the tactic of examining the
propagation of a single error, e, introduced at some stage in the computations.
Let m be an integer m,
zN@)
Tait shows that
w(n)
. \"v1 b(n-k- l)b(n -fc) • • • b(m - 1)
\\=o yN(n-fc)yN(n-fc-l)
where
KN(n): =
Vn@)
2n@)
.
--page0029--
44 Second-order homogeneous equations
Since
' yN(oV '
KN(n) is known.
Example 4.6 Consider the previous example with x = 2, N=10. Write
Using the previous results it is readily established that
, Je| yN(m + l) \\b(t^- 1) • • ¦ ft(m-l)| (l + |E*(n)|)
Now let the value 10 for yN(9) be replaced by 9, i.e., take e\"=—1. Table
4.4, also taken from Tait, displays the results obtained by using the
bounds
\\W(n)\\-
10(l + |E*(n)|)
= 9
(/\" D)
and
The above analysis can be generalized considerably. For instance by the
additivity of particular solutions we find that if errors eN(n) are
introduced at each stage of the analysis, then
= - X
and the total contribution of the error to yN(n) can be assessed by a
Table 4.4
zN(n)
zn(n)w@)/zN@)
|EN(n)|
|EN(n)| bound
11
10
9
8
7
6
0
1
9
82
665
4 737
3.040 228 070 (-6) 1.3001 (-3) 1.332 (-3)
2.769 985 575 (-5) 1.755 4 (-5) 1.824 (-5)
2.246 390 741 (-4) 3.022 6 (-7) 3.200 (-7)
1.600 173 374 (-3) 6.249 3 (-9) 7.440 (-9)
6 748 266
2.279 585 302
Error analysis, the general algorithm 45
straightforward summation. The contribution to the error in wN(n) is not
so clear because the normalization of yN(n) makes the process nonlinear.
Note that the alternative algorithm with normalization S =?c(fc)w(fc)
can easily be incorporated into the previous analysis provided only
c(k)EN(k) =
lc=0
see Tait A967) for details.
4.4 Error analysis, the general algorithm
The general Miller algorithm with unrestricted coefficients is more of a
challenge. Olver A964), however, has shown that by making a few simple
estimates of the solutions of the equation useful error bounds may still be
obtained.
Let the equation have a minimal solution u(n)= yA)(n), the solution to
be computed, and let yB)(n) be any solution dominating w(n). Define
p:=|w(N+l)/w(N)|,
c:= max |c(fc)|,
c(fc)w(fc)
lcvv(JV) k=N+1
a:=c max |w(fc)|.
Let o- be the least quantity such that the inequalities
yB)(n)
=?0-\"
hold and let
SN:= ? c(fc)yN(fc).
k-0
Then if cr, p are small (such that, e.g., pcr--page0030--
46 Second-order homogeneous equations
Olver A964) provides simple estimates for round-off error and gives
numerical examples.
The previous estimates may, depending on the problem, be very
imprecise but to do better one must work much harder. Oliver A967),
Gautschi A967, 1972a), and Mattheij and van der Sluis A976) have all
studied the error of this algorithm. The last study is the most ambitious to
date and the authors attack the problems of rounding and truncation
error for both the scalar equation and the 2x2 system with indications of
how their studies may be generalized to higher-order systems. Since the
estimates are very technical the reader is referred to the paper itself for
details.
4.5 Algorithms based on continued fractions
When N becomes large in the Miller algorithm so, in general, do the
quantities yN(n) and machine overflow becomes a very real possibility.
However, the quantities yN(n + l)/yN(n) are usually moderate in size, so
it is preferable to compute these ratios.
Gautschi A967) has analyzed an algorithm for doing this which is
mathematically equivalent to the Miller algorithm but for which overflow
can be avoided.
Let
: = ^ ?
w(n)
where w(n) is the desired solution of the equation.
From the difference equation we get
-1
(n)r(n + l)'
under appropriate conditions. Also
D.41)
so
Also
w(n)
(n) = r(n)(c(n +
? c(fc)w(fc)l,
D.42)
w@) = S/(c@) + s@)). D.43)
The algorithm results from truncating the infinite continued fraction and
Algorithms based on continued fractions 47
infinite series representing r(n) and s(n+ 1). In fact let
,,„ . , , -1 Hn) b(N-2)
rN(N): = 0; rN(n): = -— ,
a(n)-a(n + l)- a(N-l)
sN(N): = 0; sN(n): = —^-r ? c(fc)yN(fc)
Yn^J k=n + l
N
= X c(fc)rN(n)rN(n + l)---rN(fc-l).
k=n+l
The algorithm is then defined by the following computational format:
-1
rN(N) = 0; rN(n) = -
sN(N) = 0;
y -=N-l,N-2,.. . ,0;
c@) + sN@)
= rN(n-l)wN(n-l),
n = 1, 2, ...,
D.44)
As with the Miller algorithm an initial value of N is chosen for the
computations. The algorithm D4) is particularly useful when large values
of N are required to define w(n) accurately.
In the simplified Miller algorithm (c@) = l, c(fc) = 0, fc >0) all the
sN(n) vanish so the computation of sN(n) in D4) may be omitted. In this
case, the value of w@) must be known. It is clear that this algorithm is
precisely equivalent to the usual Miller algorithm and thus converges
when the conditions of Theorem 2 are satisfied.
If no a priori knowledge is available about a starting value of N, the
previous procedure has to be repeated for increasing values of N until the
required accuracy is achieved. One might think there would be some way
of employing the duality Theorem 1 to avoid repeating similar
computations. Indeed there is and the corresponding algorithm is due jointly to
Gautschi A967) and Shintani A965).
Note that one can compute r(n), s(n) by recursion for 0=?n=sK, and
thus w(n), for O^n^K, once r(K), s(K) are known for some value of K.
This is accomplished by means of the relationships D1), D2), D3),
respectively.
Let's first discuss the computation of r(K) for some fixed value of K.
To do this we exploit the connection between continued fractions and
series. The necessary formulas can be found in Wall A948, p. 17ff.).
Suppose we are given the continued fraction
K: =
aB) aC)
CB)+
.
--page0031--
48 Second-order homogeneous equations
Let
1
Then
p(k),
where p(k) satisfies
l+p(k + l) = —
with
The above formulas can be used to generate approximants to the
continued fraction, P(n)IQ(n).
Let
u(l):=l, u(k):=l + p(k), k^2,
V(k):= I»@,
i = l
Then V(k) = P(k)IQ(k) and
i?(k + l) = »(k)[u(k+
V(k + 1)= V(k) + o(k
with
We now apply this to the continued fraction B2).
Take K fixed ^ 0, and let
-l)u(k)
V(k
»(k)[u(k + l)-l
V(fc) + »(fc + l).
fc»l.
D.45)
Algorithms based on continued fractions 49
Then r(K) may be computed to the desired accuracy, assuming the
conditions of Theorem 1 are satisfied, by taking k sufficiently large, i.e.,
lim V(k) = r(K).
Now consider the computation of s(X). The quantities sN(K), N^K,
may be computed by recursion also. To do this I use the formulas B2) and
B4). Define
pN(K):=yN(K)lyN+1(K).
A straightforward application of these equations gives the format:
pK_1(X) = 0; yK(K) = l, yK_1(X) = 0;
= -[a(N -1) + b(N- 2)pN_2,
N=X
= -a(N-l)yN_1(X)-b(N-2)yN_2(X),
sN(K) = -pN_1(X)[a(N- l)sN_!
D.46)
For N sufficiently large sN(K) will be as close as one desires to s(K),
i.e.,
lim sN(K) = s(K).
To review the procedure I tabulate the steps involved in the algorithm:
A) An integer K > 0 is chosen; this means the desired solution w(n) will
be computed for the K+2 successive values 0=?n=?K + l;
B) In equation D5) take k = 1, 2,. .., until V(k) defines r(K) to the
desired number of places;
C) Take N=K + 1, K + 2,..., in equation D6) until sN(K) defines
s(K) to the desired number of places;
D) Compute r(K -1), r(K- 2),..., r@) from D1);
E) Compute s(K~ 1), s(K-2),..., s@) from D2);
F) Compute w@) from D3);
G) Compute w(l), wB),..., w(X + l) from w(n + l) = r(n)w(n).
Van der Cruyssen A981) has analyzed error propagation in this
algorithm.
.
--page0032--
50 Second-order homogeneous equations
Gautschi and Slavik A978) and Amos A974) have discussed the
computation of the modified Bessel function ratios rv(x): = Iv(x)/Iv_1(x)
by means of continued fraction algorithms, the former reference
comparing the convergence properties of two different continued fraction
representations for rv,
rv(x): =
1
2v 2(v + l)
1 1
XX
x, v>0,
h
and
rv(x): =
2x(v+D
x, v>0,
known as the Gauss and Perron continued fractions, respectively.
Asymptotically the Gauss continued fraction converges twice as fast as
Perron's. Unless very high accuracies are desired, it is the initial
convergence behavior that matters and the authors conclude that if x » v the
Perron fraction has distinct computational advantages.
In the second reference it is shown, for example, that
rv+i(x)
:rv(x),
and this estimate can be used to determine an efficient value of rv
which to begin backward recursion in the algorithm
with
2v
— + rv+1(x)
which is the formalization of Gauss' continued fraction.
Gautschi A970) discusses the application of continued fractions to the
complex error function.
The literature on computation with continued fractions is vast. The
survey by Blanch A964) is probably the best reference to consult for an
overview of the subject. Significant advances have also been made in
analyzing the truncation error for continued fractions, for instance,
Gautschi A982), Field A978), Field and Jones A972), Jones and Thron
A968, 1971), Jones and Snell A969), Sweezy and Thron A967), Henrici
and Pfluger A966), and Merkes A966).
Van der Cruyssen A979a) has generalized this algorithm by using the
concept of a generalized continued fraction. The resulting algorithm is
useful for solving recurrences of order o\\ cr3e2, for which there exist
cr— 1 linearly independent solutions that are dominated by any solution
not lying in the subspace spanned by these solutions. (This means,
Algorithms based on continued fractions 51
according to Definition 2.8, that the equation possesses a dominant
solution.) His definitions and main results are easy to state.
Definition 4.1 A generalized continued fraction (GCF) of dimension , AB>,..., A(, aB>,. . ., a(o°; b) is said to converge if all the limits
exist. The value K(am, aB>,... ,a(;b) is then defined to be the a--tuple
(d(\\) dB)d())Sr ¦
(This concept is due to de Bruijn. For a further discussion of GCF, see
ibid, the previous references.)
Van der Cruyssen establishes
Theorem 4.12 A GCF converges iff there exist a linearly independent
lti (h)() lh^v, of D7) such that:
solutions y(h)(n), l
(i) yih)(n),
-, is dominated by any solution
= o.
Now suppose we wish to compute (approximate) values of a function
.
--page0033--
52 Second-order homogeneous equations
given the initial values w@), w(l),..., w(o--l). Pick N»0 and define
W r, by
r$(N) = r%\\N) =¦¦¦= r%\\N) = 0;
ft\\n) = aA)(n)l(b(n) + r^(n + 1)), n = N- 1, N-2,..., 0;
1)), 1
Now put
wN@) = w@),
and
),..., wN(cr-1) = w(o— 1),
Van der Cruyssen establishes that if the conditions of Theorem 12 are
satisfied and, in addition,
ya\\k)
then
lim wN(n) = w(n), n
Example 4.7 Consider the recurrence
, Dn + 22)(n + 3) , ^ En+29)(n + 2)
(n+4)(n
2(n + l)(n + 6)
+ (n+4)(n + 5)
which has a fundamental set
y(w
y(n),
yA>(n)=l, yB)(n) = -^T, yC)(n) = 2\".
Here a = 2. The conditions are satisfied, and the previous algorithm will
compute successfully any solution in lin [1, nl(n + 1)]. ¦
An algorithm analogous to the one based on a three-term recurrence that
utilizes normalization relationships rather than initial values can be
formulated in a straightforward manner.
Minimal solutions and orthogonal polynomials 53
4.6 Minimal solutions and orthogonal polynomials
Gautschi A981) has observed that there is an intimate connection
between the existence of minimal solutions for the three-term recurrence
A) and the determinacy of a related moment problem. Since his work has
important implications for Gaussian quadrature, I shall present his results
in detail.
To maintain consistency with the standard treatises on orthogonal
polynomials and continued fractions I shall write the recurrence
= (z-a(n))y(n)-b(n)y(n-l),
D.48)
Obviously A) may be written in the form D8) in many ways, the role of
the variable z being discretionary.
Let /x(t) be a nondecreasing function with infinitely many points of
increase (see Henrici A977, vol. 2, p. 579)) defined on an interval [a, b],
finite or infinite, and let the moments
n):=f tndm(t), n3*0,
m(n): =
be defined.
Let pn(z) be the monic polynomials orthogonal on [a, b] with respect to
fb
Pn@Pm@dM,(t) = 0, n^m.
As is well known pn satisfies the recurrence
pn+1(z) = (z-a(n))pn(z)-b(n)pn_i(z), n5=0, j
p_!, po(z):=l, a(n) real, b(n)>0.1
The polynomials
D.49)
D.50)
are called the polynomials associated with the pn(z). They also satisfy
D9) for ns=l, but with
qo(z): = 0, q1(z):=m@).
Clearly pn(z), qn(z) are linearly independent solutions of D8).
Now let
Z
D.51)
.
--page0034--
54 Second-order homogeneous equations
F(z) then possesses the 'formal' series in powers of 1/z,
and corresponding to this series is the continued fraction
6@) 6A) 6B)
D.52)
F(z)~
z-a@)-z-a(l)-z-aB)-
6@) = m@), D.53)
where a(n), b(n) are the coefficients in D9); see Perron A957, Theorem
4.1). By corresponding I mean that when the nth convergent to the
continued fraction E3) is developed in a power series in 1/z it agrees with
the series E2) to n terms; see Henrici A977, vol. 2, p. 518).
Furthermore, the nth convergent to E3) is
6@)
6A)
z-a@)-z-a(l)-\"\"\"z-a(n-l) pn(z)'
n:
= 1.
I wish to investigate the case where the continued fraction E3) actually
converges to F(z),
F(z):= lim
Pn(z)
If E4) holds, then
/n(z):=F(z)pn(z)-qn(z),
is a minimal solution of D8) since E4) implies
D.54)
D.55)
Because of E0) and E1) fn can be written
D.56)
If we define 6@) = m@), q_i(z) = -l, then /n(z) satisfies D8) for
and the starting value is
Any condition which guarantees E4) will be sufficient for the existence
of a minimal solution to D8), fn(z), in fact. There are a number of such
conditions. For instance if [a, b] is finite E4) will hold by Markov's
theorem. If [a, b] = [0, <»] E4) will hold if the corresponding Stieltjes
moment problem is determined. (These statements can be found in
Perron A957).) A sufficient condition for the latter is due to Carleman
Minimal solutions and orthogonal polynomials 55
(Henrici A977, vol. 2, p. 608)):
For the interval [—00,00] E4) is assured if the corresponding Hamburger
moment problem is determined. Sufficient conditions for this, also due to
Carleman, are
mBk)
-1/2fc =
or
Since the moment problem for a finite interval is always determined
(Shohat and Tamarkin A943), p. 9) all cases are exhausted and I
summarize the above observations in
Theorem 4.13 The recurrence D8) possesses a minimal solution, given by
E6), whenever the moment problem for m(n) is determined. ¦
Example 4.8 (Special Gaussian quadrature rules) Let d/x(t) be a
distribution with finite or infinite support [a, 6], let x be real, x^[a, 6], and
consider the new distribution
do-(t): =
x-t
te[a,b],
[a, 6].
Suppose we wish to construct the monic orthogonal polynomials {<7rn} with
respect to do-(t) and the corresponding Gaussian quadrature rules, all of
which exist and are unique.
To do this one first determines the coefficients a(n), |3(n) in the
recurrence for the desired polynomials
irn+i(z) = (z-a rn(z) - P(n)-n-n_1(z), n 2= 0,
7T_1(z): = 0, <77-0(z):=l,
in terms of the a(n), b(n) of the given orthogonal polynomials and in
terms of the modified moments of do-(t),
Ja IA 11
7rr(z) may then be interpreted as the characteristic polynomial of the rth
leading principal minor matrix of the Jacobi matrix of order n
Hn): =
@)
L @)
V/3(n-l)
.
--page0035--
56 Second-order homogeneous equations
and, as is well-known, the required Gaussian nodes ?nv are the zeros of
•nn(z), hence the eigenvalues of J(n). The Gaussian weights (Christoffel
numbers) Anv can be written either in terms of the {irn(z)}
Anv = h(n -
where
^ v =? n,
fb
h(n-l):= irS_i(x) d--page0036--
58 Second-order homogeneous equations
When the situation is such that every fundamental set for the equation
consists of solutions having roughly the same rate of growth, so that, if
there is a minimal solution it is only weakly minimal, backward
recurrence, which is so desirable because it requires no initial values, will
fail. In this chapter we outline an approach first, apparently, used by
Clenshaw A962) in constructing the Chebyshev coefficients for the
exponential integral. In such situations the method can be used to generate
any solution of the equation. (Other techniques can be used when the
solutions of the equation are not strongly differentiated; I shall deal with
these in later chapters.)
Before I formulate the algorithm I will provide an example to show
where such a technique would be desirable.
Consider the recurrence
(n + a)(n + a + 1 - c)y (n) - (n + l)[2(n + a +1) - x - c]y (n +1)
+ (n + l)(n + 2)y(n + 2) = 0, D.58)
which has as solutions
(a)n(a + l-c)
w(n): =
n!
n\\
The recurrence has a basis of solutions {yih)(n)} having the behavior (see
Section 5.2)
ym(n)« n<4a-2c-3>/4 cos BV(nx)), yB)(n)« nDa\022c-3>/4 sin BV(nx)).
Thus the equation has no minimal solution and the attempt to generate a
solution by means of the Miller algorithm will fail. Note that E8) is a
rather significant equation. Important special cases of w@) are error,
exponential and logarithmic integrals.
To describe the Clenshaw procedure, I let {y(h\\n)} be a fundamental
set for E7) and assume, without loss of generality, that I wish to compute
ym(n).
Let constants c(h\\n) be given, h = 1, 2, ca\\n): = c(n) in equation B)
and define
S(M>:= ? ce)(k)y(h)(k), l^h, j«2.
lc=O
R(h)(n):=T('l)(n)/TA)(n) (see A-VI), h=l, 2.
The method supposes that two normalization relations are known for
yA>(n) and that the solutions {y(h\\n)} behave similarly.
The algorithm is described by the following theorem.
The Clenshaw averaging process 59
Theorem 4.14 Let
lim
N—«o
with
Let yN(n) be defined as in C), nonzero for N sufficiently large, 0=s n =sN,
SA'1)yN(n)/f cA)(k)yN(k),
lc=O
and let ya\\n) + 0 for n sufficiently large.
Then for Nt sufficiently large we can determine N2 > Nt such that the
system of equations
k=0
has a unique solution -nx, 7r2 (depending, of course, on N1; N2).
Let \\Rih\\Nj)\\l be bounded away from zero as Ni-^oo. Then
Proof The proof is deferred to Chapter 7 where the algorithm is given
for the general o-th-order equation. ¦
A minor modification of the Clenshaw method is achieved by taking
Ni = N2 and computing two sequences y%\\n), y%\\n) corresponding to
different initial values, say
and then computing an appropriate average of the sequences to satisfy
the two given normalizations. Although this modification makes the
computations more efficient it may cause numerical instabilities because
the equation for -nx generally becomes ill-conditioned as N gets large.
Example 4.10 Compute the solution w(n) of the recurrence
n +1
having the behavior
w@)=l; ?
.
--page0037--
60 Second-order homogeneous equations
Table 4.5
n
1
1
1
0 1
1 0.241721434 5 0.241657 804 2 0.241654 527 6 0.241654 506 4
2 -0.258 278 565 5 -0.258 342195 8 -0.258 347 236 7 -0.258 345 493 6
3 -0.505 519 043 6 -0.505 561463 9 -0.505 563 648 2 -0.505 563 662 4
4 -0.564 569 6414 -0.564 585 548 9 -0.564 586 368 1 -0.564586 373 4
To do this by the Clenshaw method we take
cm(k)=8k0, cB)(k) = ^.
Note this is the recurrence E8) with the specializations a = c = x = 1.
In the computations (Table 4.5) I have taken N2 = IN^ The entries in
the last column are accurate to all figures given. ¦
For the second-order recursion the Clenshaw method does not have
anything special to recommend it since a similar algorithm, based on
forward recursion, is just as stable. One simply generates two linearly
independent solutions by using the equation in the forward direction and
then determines an appropriate linear combination of these solutions to
satisfy the constraints on the desired solution.
The Clenshaw method is most useful in its application to higher-order
equations; see Section 7.4.
Wimp A967) has shown that when the equation is such that a(n), b(n)
possess asymptotic series in powers of n~v<°, o> an integer, and when
members of a certain fundamental set for the equation have similar
behavior, then the Clenshaw averaging process always converges.
5 Applications of the Miller
algorithm to the computation
of special functions
Since its inception the Miller algorithm has found a use in a wide variety
of practical problems. Among these are the computation of many of the
special functions of mathematical physics, particularly functions of
hypergeometric type, the computation of a large number of definite
integrals, especially integrals of generalized moment type, and the
computation of zeros of functions and eigenvalues of certain differential
operators.
This chapter is devoted to such applications. The treatment is intended
to be illustrative rather than exhaustive in either scope or detail. For
instance I do not generally present a complete account of the error
involved. My intention is to give an overview of the multifarious
applications of the algorithm and to emphasize the large number of areas in
science and engineering in which the algorithm has been found to be an
important computational resource. In fact the bibliographic sources for
most of the work in this chapter are not generally contained in
mathematics journals but rather in applied journals whose fields range the
spectrum from celestial mechanics to hydrology.
5.1 The confluent hypergeometric function ^»(a, c; x)
In the recurrence formula (C.16) let z -^ zb2 and
obvious relabeling of parameters and variables we obtain
•«>. After an
One solution is
Another is the function
<2). . x~TBn +
y (n):=
-a)y(n + 2) = 0. E.1)
E.2)
-2n;x). E.3)
.
--page0038--
62 Applications of the Miller algorithm
An elementary argument shows
m, , r(c)
T(a)
while
y<2)(n)~2c^3/2ex/2
ex
E.4)
E.5)
yA)(n) is a minimal solution and the recursion is stable for yB)(n).
Let us consider in detail the computation of yU)(n) by Miller's
algorithm.
A number of normalization relationships are possible; see Luke A969,
vol. 2, Ch. 9). A convenient one is
S=c-1=
-(c-l)k(c + 2fc-l)w(k).
(An obvious modification must be made of this formula when c = 1.) The
algorithm is not defined when any of the quantities c, c +1 - a, a are
negative integers or zero. However in such cases a meaningful
computational scheme may be obtained by a change of dependent variable in the
equation. I leave such details to the reader.
Theorem 4.2 is applicable and we find the Miller algorithm converges
with great rapidity. In fact formula D.21) furnishes an error estimate for
exact arithmetic. Let y(n) = y<2)(t)- The series ?k=o c(k)y(k) diverges so
rapidly that the last term serves as the leading term of its asymptotic
series. Using D), E) we find
= w(N+ l)c(N+ 1)[1 + CHAT1)].
Thus
F(a)
H±\\ (^l)N+lNa
i-N-2
N-
This result can be compared with the error bound produced by the Tait
theory, Corollary, Theorem 4.9. This result is applicable provided a(n)^
— 1 and b(n)<0. These conditions are certainly met when n is sufficiently
large, a, c real and x>0. The Tait error bound is rather complicated.
If not all members of the sequence {yn)(n)} are desired the Shintani
form of the algorithm should be used; see Theorem 4.1.
The confluent hypergeometric function 63
When c = 2a, the recursion reduces to the one satisfied by the modified
Bessel functions,
w(n)={-
1-1/2
e*'2
F(a
A great many other higher transcendental functions may be expressed
in terms of . For a partial list consult Erdelyi et al. A953, Section 6.9.2).
Currently of great interest are the Coulomb wave functions
\\T(l
Fb(ti, p) := Cn(r,)p\"+
1)!,
- ir,, 2n +2; 2ip),
(see Abramowitz and Stegun A964, Ch. 14)), where 17, p are real, p>0.
The function Fn is easily shown from our previous work to be the minimal
solution of the recurrence
5.2 The confluent hypergeometric function W(a, c; x)
One way of computing this function is to use the representation
A single-valued branch of this function is defined by choosing an
appropriate principal branch of x1 c. The parameter c, however, must be
nonintegral if this formula is to be used.
Another approach, not suffering from this disadvantage, is to start with
the recurrence (see Wimp A974))
(n + a)(n + a + l-c)y(n)-(n + l)[2(n + a +1) + x - c]y (n +1)
+ (n + l)(n + 2)y(n+2) = 0. E.6)
Solutions of this equation are
n!
and
We require a, a + l-c#O, -1,-2,.... However, these restrictions
may be removed by renormalization when necessary.
.
--page0039--
64 Applications of the Miller algorithm
Formulas in Erdelyi et al. (A953), vol. 1, 6.4, 6.6) show
n!
n!
, c; x)],
These, used with the results (ibid, 6.7) and A-V show the Casorati
determinant for the system is
D(n): = y A)(n)yB)(n +1) - y <2>(n)yA>(n +1)
D@)(a)B(a + l-c)B
E.7)
n!(n
and
The asymptotic analysis of the recurrence F) is given in Section B.2. We
have
+ Cln~
X(x) = -
r(a)F(a + l-c)
A result in Slater ), p. 80, 4.6.42) shows
|argx|<7r.
E.8)
r(c)ex/2xA~2c)/4
E.9)
> oo, |argx|<<77\\
For the preceding estimates to hold it is, strictly speaking, necessary that
none of the quantities a, c, a +1 - c be negative integers or zero. If this is
not the case the analysis must be modified. However, doing so is quite
straightforward, and details are left to the reader. For |arg x|<7r, yA>(n) is
a minimal solution and may be computed by the Miller algorithm. (Note
y<2>(n) is then a solution for which the equation is stable at any point in
the forward direction by Theorem 2.2.)
A normalization relation is provided by
E.10)
k-0
As in the previous section we can use formula D.21) to obtain an
estimate for the relative error of the Miller algorithm based on the
normalization A0). The necessary asymptotic estimates are easily made
The confluent hypergeometric function 65
and we have
e
-2>/(n). The result is
When the minus sign is affixed to the radical in the equation for f0 the
result provir^s an asymptotic formula for a constant multiple of y<2>(n).
These estimates, in turn, can be used to obtain a more accurate
asymptotic formula for ??(\")¦
Convergence for the simplified algorithm is considerably stronger,
being O(e-4V(nx>).
Example 5.1 Mori A980) has shown the functions
xa(a)n(a + l-c
a, c, x):
n!
are useful in evaluating infinite integrals. I shall give an example.
Let
= f
-) dt,
E.12)
where z is complex, |argz|<7r.
Let
. , . 2Au
t + 2k'
E.13)
where A is a fixed real number, Re z > A. The function A3) maps the right
half of the t-plane onto |u| |