Table Of ContentNumerical Solution
of
Differential Equations
Isaac Fried
Department ofMathematics
Boston University
Boston, Massachusetts
ACADEMIC PRESS
A SUBSIDIARY OFHARCOURT BRACE JOVANOVICH, PUBLISHERS
New York London Toronto Sydney San Francisco
COPYRIGHT © 1979, BY ACADEMIC PRESS, INC.
ALL RIGHTS RESERVED.
NO PART OF THIS PUBLICATION MAY BE REPRODUCED OR
TRANSMITTED IN ANY FORM OR BY ANY MEANS, ELECTRONIC
OR MECHANICAL, INCLUDING PHOTOCOPY, RECORDING, OR ANY
INFORMATION STORAGE AND RETRIEVAL SYSTEM, WITHOUT
PERMISSION IN WRITING FROM THE PUBLISHER.
ACADEMIC PRESS, INC.
Ill Fifth Avenue, New York, New York 10003
United Kingdom Edition published by
ACADEMIC PRESS, INC. (LONDON) LTD.
24/28 Oval Road, London NW1 7DX
Library of Congress Cataloging in Publication Data
Fried, Isaac.
Numerical solution of differential equations.
(Computer science and applied mathematics)
Bibliography: p.
1. Differential equations—Numerical-solutions—
Data processing. 2. Difference equations—Numerical
solutions—Data processing. 3. Finite element method—
Data Processing. I. Title.
QA371.F795 519.4 78-26454
ISBN 0-12-267780-3
PRINTED IN THE UNITED STATES OF AMERICA
82 83 84 98765432
T : T i*\
Preface
Occupying a unique place along the border between applied mathemat
ics and the concrete world of industry, the numerical solution of differen
tial equations, probably more than any other branch of numerical analy
sis, is in a constant state of unrest and evolution. Being so widely and
variously applied in the real world, its techniques are relentlessly put to
the ruthless test of practical success and usefulness. Nor does it evolve
solely through the cross influences of the practical necessities of engineer
ing; unusual impetus is also given to this field by the astounding advances
in computer technology, which is gathering now to miniaturize hardware
to lower the cost of the equipment, the arithmetic, the logic, the storage,
and the output that is made more comprehensively grasped by directly
presenting it to that most remarkable of the human senses—vision,
through computer graphics, shifting thereby the engineer's or program
mer's priorities in selecting the most appropriate solution algorithm.
Discretization of a differential equation, or for that matter the physical
process behind it, in preparation for the programming of its solution is
greatly responsible for the surge in interest in finite mathematics and
particularly linear algebra as the mathematical link between the continu
ous formulation of the problem and the properties of the computer chosen
to solve it. Linear algebra is now making its way from the periphery to the
center of the undergraduate curriculum in the sciences, and my underly
ing assumption in writing this book was that the reader is familiar with
computational linear algebra.
Several outstanding textbooks, some highly theoretical, others more
technically inclined, are available on the numerical solution of differential
equations; but they either concentrate on the initial value problem, which
XI
xii PREFACE
requires only a modicum of linear algebra, or they almost exclusively
deal, as is the case mostly, in partial differential equations with boundary
value problems in two or three dimensions. What, it is hoped, makes this
book different is that it includes at once the initial, boundary, and eigen
value problems. The book is also more emphatic on the relationship
among the analytic formulation of the physical event, the discretization
techniques applied to it, the algebraic properties of the discrete systems
created, and the properties of the digital computer. It is the proper
interaction between these components of the solution plan that decides its
ultimate success. I have put strong, maybe inordinate for the circum
stances, emphasis on the most-recent exciting and influential invention in
computational mechanics-finite elements. This was done partly for its
intrinsic usefulness in the numerical solution of differential equations and
its theoretical soundness, and partly in order to fully introduce this re
markable method in a setting simpler than partial differential equations
where its power is much more apparent but where the technique is
decidedly the same as for ordinary differential equations.
The book is intended as an introduction to the numerical solution of
differential equations but proceeds to embrace all that I deemed useful for
the person interested in the actual application of the techniques. Starting
from the fundamentals, it goes on to present the principal useful discreti
zation techniques and their theoretical aspects but without undue lengthy
and pedantic mathematical arguments, and it includes ample geometrical
and physical examples, mainly from continuum mechanics, to enhance the
reader's practical perspectives. Huge effort has been put into the de
velopment of the techniques of the numerical solution of differential
equations and their analysis. But even so, when it comes to applying them
in a realistic engineering situation, the most theoretical analysis can offer
is broad guidelines, insights, and general advice, rarely a concrete deci
sive answer. Much of the success of the numerical analyst rests with the
combination of theoretical knowledge and skill that comes from experi
ence.
The ultimate aim, of course, of all the techniques for the numerical
solution of differential equations is to render the described physical prob
lem programmable on a digital computer. This can be done either directly
from the discretization of the basic physical principles, skipping the stage
of analytical modeling and with disregard to the differential equation; and
there are engineers who prefer it this way, claiming that their responsibil
ity is to physical reality rather than to some simplified differential equa
tion. Or the starting point of the discretization can be the differential
equation with the appropriate boundary conditions in which differentials
are replaced by finite differences. This latter approach is more widespread
PREFACE xiii
and accepted, maybe for historical reasons mainly; but there has been
accumulated on differential equations a great wealth of information which
can be put to great use in the analysis of their discretization. Also, the
finite difference or finite element approximation permits greater sophisti
cation in the discretization than shear physical discrete modeling; I fol
lowed this formal discretization procedure of differential equations
throughout the book.
A complete list of references on the topics covered in this book could
fill by itself several hefty volumes, and I restricted it therefore to a mere
list of books only. As for the notation, no special typographical distinction
is made between scalars, vectors, and matrices; the reader is judged keen
enough to distinguish among them contextually.
Thanks are due to my graduate student Arthur R. Johnson who proof
read the manuscript.
j Finite Differences
1. Calculus to Algebra to Arithmetic
At the heart of all methods for the computer solution of differential equa
tions that we have in mind to discuss in this book lies the idea of computing
with small finite, rather than infinitesimal, quantities and the postponement
of the thrust to the limit, so typical of calculus, to the numerical stage of the
solution, instead of carrying it out at the stage of the analytical formulation.
Before proceeding to apply this idea to differential equations, we find it
instructive to scrutinize its working in a more familiar setting.
Consider, for that purpose, the problem of computing the area A between
the curve y = cos x and y = 0 for 0 < x ^ π/2. First we wish to apply cal
culus to the solution of this problem and effortlessly write
A = cosxdx = sin(7i/2) — sin(O) = 1 (1.1)
Jo
The solution seems simple to us, but the simpler Eq. (1.1) is the more as
tounded need we be at the power and conciseness of the symbol manipulation
of calculus. Once we remember that sin x is the antiderivative of cos x the
computation of A is reduced, according to the fundamental theorem of
integral calculus, to the mere evaluation of sin x at x = 0 and x = π/2.
It now occurs to us to engage the digital computer to compute the area
under cos x, or for that matter, any other given curve. But by its very nature,
the digital computer is not suited to performing this computation by way of
calculus. The automatic programming of the fundamental theorem of the
integral calculus on the digital computer is a formidable, if not impossible,
task of providing the computer with a symbolic list of functions and the
1
2 /. FINITE DIFFERENCES
corresponding antiderivatives (that is, a table of indefinite integrals), estab
lishing a procedure for enlarging upon this list through combinations and
substitutions, and having the computer decide when the integral can or
cannot be expressed in terms of the listed functions. Moreover, even if the
computer succeeds in expressing the integral in terms of these functions, and
we easily imagine that it will often fail, a considerable amount of arithmetic
still lies ahead since the integral, which can become vastly involved, need yet
be evaluated at the limits of integration. Integral calculus is most potent, and
hence its triumph in the precomputer days, in mental and hand computations
of reasonable integrals that can be expressed in terms of common elementary
functions. It is utterly unsuitable for automatic programming on the digital
computer.
The way to compute the area A on the present-day computer is by the
ancient precalculus technique of subdivision and summation. The domain
whose area A we wish to compute is subdivided into small finite portions, the
areas of which can be algebraically computed and summed together to pro
vide an approximation for A. A better approximation is gained with a finer
subdivision. The symbol manipulation of calculus, so alien to the digital
computer, is avoided, and the computational procedure is reduced to an
algebraic formulation for the sums from which a program is devised, in
structing the computer to perform all of the many arithmetical operations
that finally result in a number (or numbers) approximating A.
In our particular case, that of computing the area under y = cos x, the
interval 0 ^ x ^ π/2 is subdivided into N, say equal, segments of finite length
h, as shown in Fig. 1.1. A lower bound A and an upper bound A onA are com
x 2
puted by summing the areas of the rectangles and we have that
A = h(cosh + cos2h + - - - + cosNh), Nh = π/2
l
(1.2)
A = h(l + cos/i + · · · + cos(iV - \)h)
2
such that A ^ A ^ A . Admittedly the algebraic formulations for A and A
1 2 1 2
given in Eq. (1.2) are easy to program [assuming the computer knows cos(/z)],
and the limiting process or convergence is carried out by increasing N—the
number of subdivisions. The computer is incapable of performing the
theoretical approach to the limit N -> oo but in practice this is also not essen
tial. The area is desired only up to a certain practical accuracy that is mea-
Fig. 1.1 A lower bound A and an upper
x
bound A on the area under y = cosx, 0 ^ x ^
2
Ä2 π/2.
2. DIFFERENTIALS AND FINITE DIFFERENCE APPROXIMA TIONS 3
sured by how close A is to A and a reasonable accuracy requirement is
1 29
usually achieved on the computer quickly and at a fair price.
Equation (1.2) is specifically written for y = cosx, but its generalization to
any other curve is apparent. Moreover, this simple integration algorithm is
infallible: It has nothing to do with the question of whether the integral is
expressible in terms of some elementary functions or not. It might even be
cheaper to compute the area by this algorithm when the integral is expressible
in terms of tabulated functions, the values of which need be evaluated by
another complicated program. Of course the practical execution of this
program is feasible only with the fortunate help of the fast and accurate
programmable computer. Without it, by hand calculation, the accurate
numerical summation of the series (assuming no formula for the sum is found)
for A and A may well become a formidable undertaking.
1 2
With differential equations the situation is similar. The theoretical limiting
process in the mathematical formulation creates differential equations to
describe physical phenomena or engineering processes. A symbolic solution
to these equations in terms of some elementary functions, even if existing, is,
as a rule, very hard to come by and the programming necessary to obtain
such a symbolic solution impossible. To overcome this difficulty, numerical
techniques of great generality have been invented, based on the discretization
of the problem, on the division of the continuous flow of events or continuous
change of state into a series of discrete states formulated algebraically, with
the limiting process or convergence deferred to the numerical stage of the
solution.
2. Differentials and Finite Difference Approximations
To fix ideas, let us consider in detail the formation and discretization of the
differential equation describing the variation of atmospheric pressure with
height. To this end we isolate a small cylindrical volume of air shown in Fig.
1.2 and apply to it the relevant laws of nature required to establish its equilib
rium. As is commonly done by physicists, we idealize the situation in order to
simplify the mathematical analysis, and still retain a valuable model. We thus
\apz
Fig. 1.2 A cylindrical element of air between heights hx and /i2, with afhz-hi)/?g
hydrostatic and weight forces acting on it.
\hz
4 1. FINITE DIFFERENCES
assume that the air is stationary, of constant temperature, obeys the laws of an
ideal gas, and has properties that vary uniformly and continuously (un
interruptedly) with the height h. That we considerably simplify nature is
obvious. Let the small imaginary air cylinder in Fig. 1.2 be located between
the heights h and h and have a base area a. Only forces in the vertical direc
1 2
tion need be considered since the pressure p and density p of the air are
assumed to be uniform at any given height. These forces are pa exerted on
x
the bottom of the cylinder by the pressure ρ at that level, p a on top of the
γ 2
cylinder, and the weight a(h — h\)pg, in which p is the mean density and g
2
the gravitational acceleration. Inasmuch as the cylinder of air is at rest, the
algebraic sum of forces acting on it is zero resulting in the equation of
equilibrium
pa - p a - a(h - h)pg = 0 (1.3)
x 2 2 {
or
{Pi-Pi)l{h -h)=-pg (1.4)
2 l
Having assumed the air to be perfect gas, we have from Boyle's law that its
density is proportional to the pressure, or p = Rp, R being the constant of
proportionality, and Eq. (1.4) becomes
(P2-Pi)/(h -h)=-Rgp (1.5)
2 i
By our assumption that the pressure is a continuous function of h, we may
proceed to the limit h — h -> 0, that occasions p — p -> 0 and we get the
2 x 2 {
differential equation of atmospheric pressure
£--««. («)
to which we add the initial condition p = p at h = 0. There is nothing in
0
Eq. (1.6) to fret about; we know well that the symbolic solution to this initial
value problem is
P = Poe~R9h (1.7)
Nevertheless, this simple problem serves well to illustrate the way the physical
problem or differential equation is discretized and we retain it for that
purpose.
One way to discretize the change of pressure with height is to halt the
limiting process at Eq. (1.5) and consider the pressure only at discrete heights
h. The pressure p appearing in Eq. (1.5) is average and we replace it approxi
mately by (/?! + p )/2, by which the equilibrium Eq. (1.5) becomes approxi
2
mately
(Pi ~ Pi)/(h ~ Λι) = -iRg(Pi + p ) (1.8)
2 2
