Atkinson K. An Introduction to Numerical Analysis
Подождите немного. Документ загружается.
34
ERROR:
ITS SOURCES, PROPAGATION, AND ANALYSIS
For example, when xi and
Yi
are
in
single precision, then compute the products
and sums in double precision.
[On most computers, single and double precision
are fairly close in running time; although some computers do not implement
double precision in their hardware, but only in their software,
which
is
slower.]
The resulting sum
Sm
will
satisfy
(1.5.19)
a considerable improvement on (1.5.18) or (1.5.15). This can be used in parts of a
single precision calculation, significantly improving the accuracy without having
to do the entire calculation
in
double precision. For linear algebra problems, this
may halve the storage requirements
as
compared to that needed for an entirely
double precision computation.
1.6 Stability
in
Numerical Analysis
A number of mathematical problems have solutions that are quite sensitive
to
small computational errors, for example rounding errors. To deal with this
phenomenon, we introduce the concepts of
stability and condition number. The
condition number of a problem
will
be closely related to the maximum accuracy
that can be attained
in
the· solution when using finite-length numbers and
computer arithmetic. These concepts
will
then be extended to the numerical
methods that are used
to
calculate the solution. Generally
we
will want to use
numerical methods that have no greater sensitivity to small errors than was true
of the original mathematical problem.
To
simplify the presentation, the discussion
is
limited
to
problems that have
the form
of
an equation
F(x,
y)
= 0 (1.6.1)
The variable x
is
the unknown being sought, and the variable y
is
data on which
the solution depends. This equation
may
represent many different kinds of
problems.
For
example,
(1)
F may be a real valued function of the real variable
x,
and y may
be
a vector of coefficients present
in
the definition of F; or
(2)
the
equation may be an integral or differential equation, with
x an unknown function
and
y a given function or given boundary values.
We say that the problem
(1.6.1)
is
stable if the solution x depends in a
continuous way on the variable
y. This means that if {
Yn}
is
a sequence of values
approaching
y in some sense, then the associated solution values { x n } must also
approach
x in some
way.
Equivalently, if
we
make ever smaller changes in y,
these must lead to correspondingly smaller changes in x. The sense in which the
changes are small
will
depend on the norm being used
to
measure the sizes of the
vectors
x and
y;
there are many possible choices, varying with the problem.
Stable problems are also called
well-posed problems, and
we
will
use the two
terms interchangeably.
If
a problem
is
not stable, it
is
called unstable or ill-posed.
Example (a) Consider the solution of
ax
2
+ bx + c = 0
a=I=O
STABILITY IN NUMERICAL ANALYSIS 35
Any solution x
is
a complex number. For the data
in
this case,
we
use
y
=(a,
b,
c), the vector of coefficients.
It
should be clear from the quadratic
formula
x=
-b
±
Vb
2
-
4ac
2a
that the two solutions for x
will
vary in a continuous
way
with the data
y = (a,
b,
c).
(b) Consider the integral equation problem
1
.?5x(t)dt
1
=
y(s)
0 1.25 -
cos
( 2
7T
( s +
t))
O~s~1
(1.6.2)
This
is
an unstable problem. There are perturbations on(s) =
Yn(s)-
y(s)
for
which
Max
lon(s) I
~
0
O:Ss:Sl
and the corresponding solutions xn(s) satisfy
Max I X n ( s) - X ( s) I = 1
O:Ss:Sl
Specifically, define Yn(s) =
y(s)
+ on(s)
1
as
n
~
oo
all n
~
1
on(s) =
2
n COS
(2n7TS)
O~s~1
n~1
Then it can be shown that
thus proving (1.6.4).
(1.6.3)
(1.6.4)
If
a problem (1.6.1)
is
unstable, then there are serious difficulties in attempting
to solve it.
It
is
usually not possible
to
solve such problems without first
attempting to understand more about the properties of the solution, usually
by
returning to the context in which the mathematical problem
was
formulated. This
is
currently a very active area of research in applied mathematics and numerical
analysis [see, for example, Tikhonov and Arsenio (1977) and Wahba (1980)].
For practical purposes there are many problems that are stable in the
previously given sense, but that are still very troublesome
as
far
as
numerical
computations are concerned. To deal with this difficulty,
we
introduce a measure
of stability called a
condition number.
It
shows that practical stability represents a
continuum of problems, some better behaved than others.
The condition number attempts to measure the worst possible effect on the
solution
x of (1.6.1) when the variable y
is
perturbed by a small amount. Let oy
36
ERROR.: ITS SOURCES, PROPAGA1]0N, AND ANALYSIS
be
a perturbation
of
y, and let x +
ox
be the solution of the perturbed equation
F(x +ox, y + oy) = 0
(1.6.5)
Define
(1.6.6)
We
have used
the
notation
II
·II
to denote a measure of size. Recall the definitions
(1.1.16)-(1.1.18) for vectors from
R"
and C[ a, b ].-The-example (1.6.2) used the
norm
(1.1.18) for
mea.Suring
the perturbations in both x
andy.
Commonly x
and
y may be different kinds of variables, and then different norms are appropriate.
The
supremum
in
(1.6.6) is taken over all small perturbations o y for which the
perturbed problem
(1.6.5).will still make sense. Problems that are unstable lead to
K(x)
= oo.
The
number
K (
x)
is
called the condition number for (1.6.1
).
It
is a measure of
the sensitivity
of
the solution x to small changes in the data
y.
If
K(x)
is
quite
large, then there exists small relative changes
8 y in y that lead to large relative
changes
8x
in
x. But if
K(x)
is
small, say
K(x)
~
10,
then small relative changes
in
y always lead to correspondingly small relative changes in x. Since numerical
calculations almost always involve a variety of small computational errors, we do
not
want problems with a large condition numbe:r. Such problems are called
ill-conditioned,
and
they are generally very hard to solve accurately.
Example Consider solving
Perturbing y
by
oy,
we
have
8x
X
X-
aY
= 0
a>O
ay+By-
aY
----
=a
8
Y-1
aY
For
the condition number for (1.6.7).
Restricting 8 y
to
be small,
we
have
K(x).=jy
·ln(a)l
(1.6.7)
(1.6.8)
Regardless
of
how
we
compute x in (1.6.1),
if
K(x)
is large, then small relative
changes
in
y will lead to much larger relative changes in x.
If
K(x)
=
10
4
and
if
the value
of
y being used has relative error
10-
7
due to using finite-length
computer
arithmetic and rounding, then
it
is
likely that the resulting value
of
x
will
have relative error of about
10-
3
•
This
is
a large drop in accuracy,
and
there
STABILITY
IN
NUMERICAL ANALYSIS 37
is little way to avoid
it
except perhaps by doing all computations in longer
precision computer arithmetic, provided
y can then be obtained with greater
accuracy.
Example
Consider the n X n nonsingular matrix
1
1
1
1
-
-
2
3
n
1 1
1
1
Y=
2 3 4
n + 1
(1.6.9)
1 1
1
n
n+1
2n-
1
which is called the Hilbert matrix. The problem
of
calculating the inverse of
Y,
or
equivalently
of
solvirig
YX
= I with I the identity matrix,
is
a well-posed
problem. The solution
X can be obtained in a finite number of steps using only
simple-arithmetic operations. But the problem
of
calculating X is increasingly
·ill-conditioned -as -n -increases.
The ill-conditioning of the numerical inversion of
Y will be shown in a
practical setting. Let
Y denote the result of entering the matrix Y into an
IBM
370 computer and storing the matrix entries using single precision floating-point
format. The
fraction,al elements of Y will be expanded in the hexadecimal (base
16) number system and then chopped after six hexadecimal digits (about seven
decimal digits). Since most of the .entries in
Y
do
not have finite hexadecimal
expansions, there
will be a relative error of about
10-
6
in each such element
of
Y.
Using higher precision arithmetic,
we
can calculate the exact value of
f-
1
.
The inverse
y-
1
is known analytically, and
we
can compare it with
f-
1
.
For
n =
6,
some of the elements of
f-
1
differ from the corresponding elements in
y-
1
in the first nonzero digit. For example, the entries in row
6,
column 2 are
(Y-
1
)6,2
= 73866.34
This makes the calculation of
y-
1
an ill-conditioned problem, and it becomes
increasingly
·so
as
n increases. The condition number in (1.62) will
be
at
least 10
6
as a reflection of the poor accuracy in y-
1
compared with
y-
1
•
Lest this be
thought of as an odd pathological example that
oould not occur in practice, this
particular example occurs when doing least squares approximation theory (e.g.,
see Section 4.3). The general area of ill-conditioned problems for linear systems
and matrix inverses
is
considered in greater detail in Chapter
8.
Stability of
·numerical
algorithms
A numerical method for solving a mathemati-
cal problem is considered stable if the sensitivity of the numerical answer to the
data
is
no greater than in the original mathematical problem. We will make this
more precise, again using
(1.6.1)
as
a model for the problem. A numerical
methOd
38 ERROR: ITS SOURCES, PROPAGATION, AND ANALYSIS
for solving (1.6.1)
will
generally result in a sequence of approximate problems
(1.6.10)
depending on some parameter, say
n.
The data
Yn
are to approach y
as
n
~
oo;
the function values Fn(z,
w)
are to approach F(z, w)
as
n
~
oo,
for all (z,
w)
near
(x,
y);
and
hopefully the resulting approximate solutions
xn
will approach
x as n
~
oo.
For
example, (1.6.1) may represent a differential equation initial
value problem, and
(1.6.10) may present a sequence of finite-difference approxi-
mations depending on h =
1/n,
as
in and following (1.3.5). Another case would
be
where n represents the number of digits being used in the calculations, and
we
may
be
solving
F(x,
y)
= 0 as exactly
as
possible within this finite precision
arithmetic.
For
each
of
the problems (1.6.10)
we
can define a condition number Kn(xn),
just
as in (1.6.6). Using these cbndition numbers, define
K(x)
= Limit Supremum Kk(xk) (1.6.11)
n->oo
k"?!n
We say the numerical method is stable if
K(x)
is
of about the same magnitude
as
K (
x)
from (1.6.6), for example, if
K(x)::;;
2K(x)
If
this is true, then the sensitivity of (1.6.10) to changes in the data
is
about the
same as that
of
the original problem (1.6.1).
Some problems and numerical methods may not
fit
easily within the frame-
work of (1.6.1), (1.6.6), (1.6.10), and (1.6.11), but there is a general idea of stable
problems
and·
condition numbers that can be introduced and given similar
meaning. The main use of these concepts in this text
is
in
(1) rootfinding for
polynomial equations, (2) solving differential equations, and (3) problems in
numerical linear algebra. Generally there
is
little problem with unstable numeri-
cal methods in this text. The main difficulty will be the solution of ill-conditioned
problems.
·
E:mmple Consider the evaluation of a
Bessel
function,
m
~
0
(1.6.12)
This series converges very rapidly, and the evaluation of x
is
easily shown to be a
well-conditioned problem in its dependence
on
y.
Now consider the evaluation of Jm(Y) using the triple recursion relation
m~1
(1.6.13)
assuming J
0
(y)
and J
1
(y)
are known.
We
now demonstrate numerically that this
DISCUSSION OF THE LITERATURE
39
Table
1.5 Computed values of
J,.(l)
m Computed
Jm
(1) True Jm(1)
0
.7651976866 .7651976866
1
.4400505857 .4400505857
2
.1149034848 .1149034849
3
.195633535E-
1
.1956335398E - 1
4
.2476636iE - 2 .2476638964E - 2
5
.2497361E - 3 .2497577302E - 3
6
.207248E-
4
.2093833800E - 4
7
-.10385E-
5 .1502325817E - 5
is
an
unstable numerical method for evaluating Jm(y), for even moderately large
m.
We
take
y = 1, so
that
(1.6.13) becomes
m
~
1
(1.6.14)
We
use
values for 1
0
(1)
and
1
1
(1)
that
are
accurate
to
10 significant digits.
The
subsequent
values Jm(l) are calculated from (1.6.4) using exact arithmetic,
and
the
results
are
given
in
Table
1.5.
The
true
values
are
given for comparison,
and
they
show
the
rapid
divergence
of
the
approximate
values from the
true
values:
The
only
errors
introduced were the rounding
errors
in
1
0
(1)
and
1
1
(1),
and
they
cause
an
increasingly large
perturbation
in
Jm(l)
as
m increases.
The
use
of
three-term recursion relations
m
~
1
is a
common
tool in applied mathematics
and
numerical analysis.
But
as
previously shown, they
can
lead
to
unstable
numerical
methods.
For
a general
analysis
of
triple recursion relations, see Gautschi. (1967).
In
the
case
of
(1.6.13)
and
(1.6.14), large loss
of
significance errors are occurring.
Discussion of the Literature
A
knowledge
of
computer
arithmetic is
important
for
programmers
who
are
concerned
with
numerical accuracy, particularly
when
writing
programs
that
are
to
be
widely used. Also, when writing programs
to
be
run
on
various computers,
their different floating-point characteristics
must
be
taken
into
account. Classic
treatments
of
floating-point arithmetic are given
in
Knuth
(1981, chap. 4)
and
Sterbenz
(1974).
The
topic
of
error
propagation, especially
that
due
to
rounding/chopping
error,
has
been
difficult
to
treat
in
a precise,
but
useful manner.
There
are
some
important
early papers,
but
the
current
approaches
to
the
subject are
due
in
large
40 ERROR: ITS SOURCES, PROPAGATION, AND ANALYSIS
part
to the late
J.
H. Wilkinson. Much
of
his work was in numerical linear
algebra,
but
he
made
important contributions to
many
areas of numerical
analysis.
For
a general introduction to his techniques for analyzing the propa-
gation
of
errors, with applications to several
important
problems, see Wilkinson
(1963), (1965), (1984).
Another
approach
to the control
of
error is called interval analysis. With it, we
carry along
an
interval [x
1
,
x,]
in
our
calculations,
rather
than a single number
xA,
and
the
numbers
x
1
and x., are guaranteed to
bound
the true value
Xr~
The
difficulty
with
this approach
is
that the size
of
x., - x
1
is generally much larger
than
lxr-
xAI,
mainly because the possible cancellation
of
errors of opposite
sign is often
not
consideredwhen~corp.puting
x
1
and
x,.
For
an
introduction to
this area, showing how to improve
on
these conservative bounds in particular
cases, see
Moore
(1966). More recently, this area
and
that
of
computer arithmetic
have
been
combined to give a general theoretical framework allowing the
development
of
algorithms with rigorous error bounds. As examples
of
this area,
see
the
texts
of
Alefeld
and
Herzberger (1983),
and
Kulisch
and
Miranker (1981),
the symposium proceedings of Alefeld and Grigorieff (1980), and the survey in
Moore
(1979).
The
topic
of
ill-posed problems was
just
touched
on
in
Section 1.6,
but
it
has
been
of
increasing interest in recent years. There
are
many
problems
of
indirect
physical measurement that lead to ill-posed problems,
and
in
this form they are
called
inverse problems. The book
by
Lavrentiev (1967) gives a general introduc-
tion, although
it
discusses mainly (1) analytic continuation
of
analytic functions
of
a complex variable, and (2) inverse problems for differential equations.
One
of
the
major
numerical tools used in dealing with ill-posed problems is called
regularization,
and
an extensive development
of
it is given
in
Tikhonov
and
Arsenin (1977).
As
important examples
of
the
more current literature
on
numeri-
cal
methods
for ill-posed problems, see Groetsch (1984)
and
Wahba (1980).
Two
new types
of
computers have appeared
in
the last 10 to
15
years,
and
they
are
now
having
an
increasingly important impact
on
numerical analysis. These
are
microcomputers and supercomputers. Everyone is aware
of
microcomputers;
scientists
and
engineers are using them for
an
increasing amount
of
their
numerical calculations. Initially the arithmetic design
of
microcomputers was
quite
poor, with some having errors in their basic arithmetic operations. Re-
cently,
an
excellent new standard has been produced for arithmetic
on
microcom-
puters,
and
with
it
one can write high-quality
and
efficient numerical programs.
This
standard,
the
IEEE
Standard for Binary Floating-Point Arithmetic, is de-
scribed
in
IEEE
(1985). Implementation
on
the major families
of
microprocessors
are
becoming available; for example, see Palmer
and
Morse (1984).
The
name
supercomputer refers to a variety
of
machine designs, all having
ii:t
common
the
ability to
do
very high-speed numerical computations, say greater
than
20 million floating-point operations per second. This
area
is·
developing
and
changing very rapidly,
and·so
we can only give a few references to hint
at
the
effect
of
these machines
on
the design
of
numerical algorithms. Hackney
and
Jesshope (1981), and Quinn (1987) are general texts
on
the architecture
of
supercomputers
and
the design
of
numerical algorithms
on
them; Parter (1984) is
a symposium proceedings giving some applications
of
supercomputers
in
a
variety
of
physical problems; and. Ortega
and
Voigt
(19~5)
discuss supercom-
BIBLIOGRAPHY
41
puters
as they are being used to solve partial differential _equations. These
machines
will become increasingly
important
in
all areas of computing, and their
architectures
are
likely to affect smaller mainframe computers
of
the type
that
are
now more widely used.
Symbolic mathematics is a rapidly growing area,
and
with
it
one
can
do
analytic
rather
than numerical mathematics, for example, finding antiderivatives
exactly when possible.
lbis
area has
not
significantly affected numerical analysis
to date,
but
that
appears to
be
changing.
In
many
situations, symbolic mathe-
matics are used for
part
of
a calculation, with numerical methods used for the
remainder
of
the calculation.
One
of
the most sophisticated
of
the programming
languages for carrying
out
symbolic mathematics is MACSYMA, which is
described
in
Rand
(1984).
For
a survey and historical account
of
programming
languages for this area, see Van Hulzen
and
Calmet (1983).
We
conclude by discussing the area
of
mathematical software. This area deals
with the implementation
of
numerical algorithms as computer programs, with
careful attention given to questions
of
accuracy, efficiency, flexibility, portability,
and
other
characteristics that improve the usefulness
of
the programs. A major
journal
of
the
area is the
ACM
Transactions
on
Mathematical Software.
For
an
extensive survey of the area, including the most
important
program libraries
that
have been developed
in
recent years, see Cowell (1984).
In
the ap_pendix to this
book, we give a further discussion
of
some currently available numerical analysis
computer
program packages.
Bibliography
Alefeld, G.,
and
R.
Grigorieff, eds. (1980). Fundamentals
of
Numerical Computa-
tion (Computer-oriented Numerical Analysis). Computing Supplementum 2,
Springer-Verlag, Vienna.
·
Alefeld, G.,
and
J. Herzberger (1983). Introduction
to
Interval Computations.
Academic Press, New York.
Bender, E. (1978).
An
Introduction to Mathematical Modelling. Wiley,
New
York.
Cowell, W., ed. (1984). Sources and Development
of
Mathematical Software.
Prentice-Hall, Englewood Cliffs, N.J.
Fadeeva,
V.
(1959). Computational Methods
of
Linear Algebra. Dover,
New
York.
Forsythe,
G.
(1969). What is a satisfactory quadratic equation solver?
In
B.
Dejon
and
P.
Henrici (eds.), Constructive Aspects
of
the Fundamental
Theorem
of
Algebra, pp.
53-61,
Wiley, New York.
Forsythe, G., and
C.
Moler (1967). Computer Solution
of
Linear Algebraic
Systems.
Prentice-Hall, Englewood Cliffs,
NJ.
Forsythe, G., M. Malcolm,
and
C. Moler (1977). Computer Methods for Mathe-
matical Computations. Prentice-Hall, Englewood Cliffs,
NJ.
Gautschi,
W.
(1967). Computational aspects
of
three term recurrence relations.
SIAM
Rev., 9, 24-82.
42 ERROR: ITS SOURCES, PROPAGATION, AND ANALYSIS
Groetsch,
C.
(1984). The Theory
of
Tikhonov Regularization for Fredholm Equa-
tions
of
the First Kind. Pitman, Boston .
. Henrici,
P.
(1962). Discrete Variable Methods
in
Ordinary Differential Equations.
Wiley, New York.
Hockney, R.,
and
C.
Jesshope (1981). Parallel Computers: Architecture, Program-
ming, and Algorithms.
Adam Hilger, Bristol, England.
Hogg, R., and
A.
Craig
(1978).
Introduction
to
Mathematical Statistics, 4th
ed.
Macmillan, New York.
Institute
of
Electrical and Electronics Engineers
(1985).
Proposed standard for
binary floating-point arithmetic. (IEEE Task
P754),
Draft 10.0. IEEE
Society, New
York.
Knuth, D. (1981). The
Art
of
Computer Programming,
vol.
2,
Seminumerica/
Algorithms,
2nd ed. Addison-Wesley, Reading, Mass.
Kulisch,
U.,
and
W.
Miranker (1981). Computer Arithmetic in Theory and Prac-
tice.
Academic Press, New York.
Lavrentiev, M. (1967). Some Improperly Posed Problems
of
Mathematical Physics.
Springer-Verlag, New York.
Lin, C., and L. Segal (1974). Mathematics Applied
to
Deterministic Problems in the
Natural Sciences.
Macmillan, New York.
Maki, D., and
M.
Thompson (1973). Mathematical Models and Applications.
Prentice-Hall, Englewood Cliffs, N.J.
Moore, R. (1966).
Interval Analysis. Prentice-Hall, Englewood Cliffs, N.J.
Moore,
R.
(1979). Methods and Applications
of
Interval Analysis. Society for
Industrial and Applied Mathematics, Philadelphia.
Ortega, J., and
R.
Voigt
(1985).
Solution of partial differential equations on
vector and parallel computers.
SIAM
Rev.,
27,
149-240.
Parter,
S.,
ed. (1984). Large Scale Scientific Computation. Academic Press, New
York.
Pahner, J., and
S.
Morse (1984). The 8087 Primer.
Wiley,
New York.
Quinn,
M.
(1987). Designing Efficient Algorithms for Parallel Computers.
McGraw-Hill, New York.
Rand, R. (1984). Computer Algebra
in
Applied Mathematics:
An
Introduction to
MA
CSYMA.
Pitman, Boston.
Rubinow,
S.
(1975). Introduction
to
Mathematical Biology. Wiley, New York.
Sterbenz,
P.
(1974). Floating-Point Computation. Prentice-Hall, Englewood Cliffs,
N.J.
Tikhonov,
A., and
V.
Arsenin (1977). Solutions
of
Ill-posed Problems. Wiley, New
York.
Van
Hulzen, J., and J. Calmet (1983). Computer algebra systems. In
B.
Buch-
berger, G. Collins,
R.
Loos (eds.), Computer Algebra: Symbolic and Alge-
braic Computation,
2nd ed. Springer-Verlag, Vienna.
Wahba, G. (1980). Ill-posed problems: Numerical and statistical methods for
mildly, moderately, and severely ill-posed problems with noisy data. Tech.
PROBLEMS
43
Rep. # 595, Statistics Department, Univ. of Wisconsin, Madison. Prepared
for the
Proc. Int. Symp.
on
Ill-Posed Problems, Newark, Del., 1979.
Wahba, G. (1984). Cross-validated spline methods for the estimation of multi-
variate functions from data on functionals, in
Statistics:
An
Appraisal,
H.
A. David and
H.
T.
David (eds.), Iowa State Univ. Press, Ames,
pp.
205-235.
Wilkinson, J. (1963). Rounding Errors
in
Algebraic Processes. Prentice-Hall,
Englewood
Cliffs,
N.J.
Wilkinson,
J.
(1965). The Algebraic Eigenvalue Problem. Oxford, England.
Wilkinson,
J. (1984). The perfidious polynomial. In Golub (ed.), Studies in
Num.erical Analysis,
pp.
1-28,
Mathematical Association of America.
Problems
1.
(a) Assume
f(x)
is
continuous on
as
x s
b,
and consider the average
1 n
S = - L
f(xj)
n
j=l
with all points
xj
in the interval [a, b
].
Show
that
S=f(r)
for some r in [a, b]. Hint: Use the intermediate value theorem
and
consider the range of values of
f(x)
and thus of S.
(b)
Generaliz~
part
(a)
to the sum
n
S=
L:wJ(x)
j-1
with all xj in
[a,
b)
and all
wj
~
0.
2. Derive the following inequalities:
for
all x, z
.s
0.
'IT
'IT
(b)
lx-
zi
s
ltan(x)-
tan(z)l
--<x
z<-
2 ' 2
0 s y s
x,
p
~
1.
3. (a) Bound the error in the approximation
sin(x)
= x
lxl s 8