# Deturck, Wilf - Lecture Notes on Numerical Analysis

Ïîñìîòðåòü àðõèâ öåëèêîìÏðîñìîòðåòü ôàéë â îòäåëüíîì îêíå: 14606f2d29e9a97367d71dae19dd4662.pdf

Lectures on Numerical Analysis

Dennis Deturck and Herbert S. Wilf

Department of Mathematics

University of Pennsylvania

Philadelphia, PA 19104-6395

Copyright 2002, Dennis Deturck and Herbert Wilf

April 30, 2002

2

Contents

1 Differential and Difference Equations

5

1.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2

Linear equations with constant coefficients . . . . . . . . . . . . . . . . . . .

8

1.3

Difference equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

1.4

Computing with difference equations . . . . . . . . . . . . . . . . . . . . . .

14

1.5

Stability theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

1.6

Stability theory of difference equations . . . . . . . . . . . . . . . . . . . . .

19

2 The Numerical Solution of Differential Equations

23

2.1

Euler’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

2.2

Software notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

2.3

Systems and equations of higher order . . . . . . . . . . . . . . . . . . . . .

29

2.4

How to document a program . . . . . . . . . . . . . . . . . . . . . . . . . .

34

2.5

The midpoint and trapezoidal rules . . . . . . . . . . . . . . . . . . . . . . .

38

2.6

Comparison of the methods . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

2.7

Predictor-corrector methods . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

2.8

Truncation error and step size . . . . . . . . . . . . . . . . . . . . . . . . . .

50

2.9

Controlling the step size . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

2.10 Case study: Rocket to the moon . . . . . . . . . . . . . . . . . . . . . . . .

60

2.11 Maple programs for the trapezoidal rule . . . . . . . . . . . . . . . . . . . .

65

2.11.1 Example: Computing the cosine function . . . . . . . . . . . . . . .

67

2.11.2 Example: The moon rocket in one dimension . . . . . . . . . . . . .

68

2.12 The big leagues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

2.13 Lagrange and Adams formulas . . . . . . . . . . . . . . . . . . . . . . . . .

74

4

CONTENTS

3 Numerical linear algebra

81

3.1

Vector spaces and linear mappings . . . . . . . . . . . . . . . . . . . . . . .

81

3.2

Linear systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

86

3.3

Building blocks for the linear equation solver . . . . . . . . . . . . . . . . .

92

3.4

How big is zero? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

3.5

Operation count . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

102

3.6

To unscramble the eggs . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

105

3.7

Eigenvalues and eigenvectors of matrices . . . . . . . . . . . . . . . . . . . .

108

3.8

The orthogonal matrices of Jacobi . . . . . . . . . . . . . . . . . . . . . . .

112

3.9

Convergence of the Jacobi method . . . . . . . . . . . . . . . . . . . . . . .

115

3.10 Corbat?o’s idea and the implementation of the Jacobi algorithm . . . . . . .

118

3.11 Getting it together . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

122

3.12 Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

124

Chapter 1

Differential and Difference

Equations

1.1

Introduction

In this chapter we are going to study differential equations, with particular emphasis on how

to solve them with computers. We assume that the reader has previously met differential

equations, so we’re going to review the most basic facts about them rather quickly.

A differential equation is an equation in an unknown function, say y(x), where the

equation contains various derivatives of y and various known functions of x. The problem

is to “find” the unknown function. The order of a differential equation is the order of the

highest derivative that appears in it.

Here’s an easy equation of first order:

y 0 (x) = 0.

(1.1.1)

The unknown function is y(x) = constant, so we have solved the given equation (1.1.1).

The next one is a little harder:

y 0 (x) = 2y(x).

(1.1.2)

A solution will, now doubt, arrive after a bit of thought, namely y(x) = e2x . But, if y(x)

is a solution of (1.1.2), then so is 10y(x), or 49.6y(x), or in fact cy(x) for any constant c.

Hence y = ce2x is a solution of (1.1.2). Are there any other solutions? No there aren’t,

because if y is any function that satisfies (1.1.2) then

(ye?2x )0 = e?2x (y 0 ? 2y) = 0,

(1.1.3)

and so ye?2x must be a constant, C.

In general, we can expect that if a differential equation is of the first order, then the

most general solution will involve one arbitrary constant C. This is not always the case,

6

Differential and Difference Equations

since we can write down differential equations that have no solutions at all. We would have,

for instance, a fairly hard time (why?) finding a real function y(x) for which

(y 0 )2 = ?y 2 ? 2.

(1.1.4)

There are certain special kinds of differential equations that can always be solved, and

it’s often important to be able to recognize them. Among there are the “first-order linear”

equations

y 0 (x) + a(x)y(x) = 0,

(1.1.5)

where a(x) is a given function of x.

Before we describe the solution of these equations, let’s discuss the word linear. To say

that an equation is linear is to say that if we have any two solutions y1 (x) and y2 (x) of the

equation, then c1 y1 (x) + c2 y2 (x) is also a solution of the equation, where c1 and c2 are any

two constants (in other words, the set of solutions forms a vector space).

Equation (1.1.1) is linear, in fact, y1 (x) = 7 and y2 (x) = 23 are both solutions, and so

is 7c1 + 23c2 . Less trivially, the equation

y 00 (x) + y(x) = 0

(1.1.6)

is linear. The linearity of (1.1.6) can be checked right from the equation itself, without

knowing what the solutions are (do it!). For an example, though, we might note that

y = sin x is a solution of (1.1.6), that y = cos x is another solution of (1.1.6), and finally,

by linearity, that the function y = c1 sin x + c2 cos x is a solution, whatever the constants c1

and c2 . Now let’s consider an instance of the first order linear equation (1.1.5):

y 0 (x) + xy(x) = 0.

(1.1.7)

So we’re looking for a function whose derivative is ?x times the function. Evidently y =

2

2

e?x /2 will do, and the general solution is y(x) = ce?x /2 .

If instead of (1.1.7) we had

y 0 (x) + x2 y(x) = 0,

then we would have found the general solution ce?x

As a last example, take

3 /3

.

y 0 (x) ? (cos x) y(x) = 0.

(1.1.8)

The right medicine is now y(x) = esin x . In the next paragraph we’ll give the general rule

of which the above are three examples. The reader might like to put down the book at this

point and try to formulate the rule for solving (1.1.5) before going on to read about it.

Ready? What we need is to choose some antiderivative A(x) of a(x), and then the

solution is y(x) = ce?A(x) .

Since that was so easy, next let’s put a more interesting right hand side into (1.1.5), by

considering the equation

y 0 (x) + a(x)y(x) = b(x)

(1.1.9)

1.1 Introduction

7

where now b(x) is also a given function of x (Is (1.1.9) a linear equation? Are you sure?).

To solve (1.1.9), once again choose some antiderivative A(x) of a(x), and then note that

we can rewrite (1.1.9) in the equivalent form

e?A(x)

d A(x)

e

y(x) = b(x).

dx

Now if we multiply through by eA(x) we see that

d A(x)

e

y(x) = b(x)eA(x)

dx

(1.1.10)

so , if we integrate both sides,

Z

eA(x) y(x) =

x

b(t)eA(t) dt + const. ,

(1.1.11)

where on the right side, we mean any antiderivative of the function under the integral sign.

Consequently

Z x

y(x) = e?A(x)

b(t)eA(t) dt + const. .

(1.1.12)

As an example, consider the equation

y0 +

y

= x + 1.

x

(1.1.13)

We find that A(x) = log x, then from (1.1.12) we get

Z

x

1

y(x) =

(t + 1)t dt + C

x

x2 x C

=

+ + .

3

2

x

(1.1.14)

We may be doing a disservice to the reader by beginning with this discussion of certain

types of differential equations that can be solved analytically, because it would be erroneous

to assume that most, or even many, such equations can be dealt with by these techniques.

Indeed, the reason for the importance of the numerical methods that are the main subject

of this chapter is precisely that most equations that arise in “real” problems are quite

intractable by analytical means, so the computer is the only hope.

Despite the above disclaimer, in the next section we will study yet another important

family of differential equations that can be handled analytically, namely linear equations

with constant coefficients.

Exercises 1.1

1. Find the general solution of each of the following equations:

(a) y 0 = 2 cos x

8

Differential and Difference Equations

2

y=0

x

(c) y 0 + xy = 3

1

(d) y 0 + y = x + 5

x

(e) 2yy 0 = x + 1

(b) y 0 +

2. Show that the equation (1.1.4) has no real solutions.

3. Go to your computer or terminal and familiarize yourself with the equipment, the

operating system, and the specific software you will be using. Then write a program

that will calculate and print the sum of the squares of the integers 1, 2, . . . , 100. Run

this program.

4. For each part of problem 1, find the solution for which y(1) = 1.

1.2

Linear equations with constant coefficients

One particularly pleasant, and important, type of linear differential equation is the variety

with constant coefficients, such as

y 00 + 3y 0 + 2y = 0 .

(1.2.1)

It turns out that what we have to do to solve these equations is to try a solution of a certain

form, and we will then find that all of the solutions indeed are of that form.

Let’s see if the function y(x) = e?x is a solution of (1.2.1). If we substitute in (1.2.1),

and then cancel the common factor e?x , we are left with the quadratic equation

?2 + 3? + 2 = 0

whose solutions are ? = ?2 and ? = ?1. Hence for those two values of ? our trial function

y(x) = e?x is indeed a solution of (1.2.1). In other words, e?2x is a solution, e?x is a

solution, and since the equation is linear,

y(x) = c1 e?2x + c2 e?x

(1.2.2)

is also a solution, where c1 and c2 are arbitrary constants. Finally, (1.2.2) must be the most

general solution since it has the “right” number of arbitrary constants, namely two.

Trying a solution in the form of an exponential is always the correct first step in solving

linear equations with constant coefficients. Various complications can develop, however, as

illustrated by the equation

y 00 + 4y 0 + 4y = 0 .

(1.2.3)

Again, let’s see if there is a solution of the form y = e?x . This time, substitution into

(1.2.3) and cancellation of the factor e?x leads to the quadratic equation

?2 + 4? + 4 = 0,

(1.2.4)

1.2 Linear equations with constant coefficients

9

whose two roots are identical, both being ?2. Hence e?2x is a solution, and of course so is

c1 e?2x , but we don’t yet have the general solution because there is, so far, only one arbitrary

constant. The difficulty, of course, is caused by the fact that the roots of (1.2.4) are not

distinct.

In this case, it turns out that xe?2x is another solution of the differential equation (1.2.3)

(verify this), so the general solution is (c1 + c2 x)e?2x .

Suppose that we begin with an equation of third order, and that all three roots turn

out to be the same. For instance, to solve the equation

y 000 + 3y 00 + 3y 0 + y = 0

(1.2.5)

we would try y = e?x , and we would then be facing the cubic equation

?3 + 3?2 + 3? + 1 = 0 ,

(1.2.6)

whose “three” roots are all equal to ?1. Now, not only is e?x a solution, but so are xe?x

and x2 e?x .

To see why this procedure works in general, suppose we have a linear differential equation

with constant coeficcients, say

y (n) + a1 y (n?1) + a2 y (n?2) + · · · + an y = 0

(1.2.7)

If we try to find a solution of the usual exponential form y = e?x , then after substitution into

(1.2.7) and cancellation of the common factor e?x , we would find the polynomial equation

?n + a1 ?n?1 + a2 ?n?2 + · · · + an = 0 .

(1.2.8)

The polynomial on the left side is called the characteristic polynomial of the given

differential equation. Suppose now that a certain number ? = ?? is a root of (1.2.8) of

multiplicity p. To say that ?? is a root of multiplicity p of the equation is to say that

(? ? ?? )p is a factor of the characteristic polynomial. Now look at the left side of the given

differential equation (1.2.7). We can write it in the form

(D n + a1 D n?1 + a2 D n?2 + · · · + an )y = 0 ,

(1.2.9)

in which D is the differential operator d/dx. In the parentheses in (1.2.9) we see the

polynomial ?(D), where ? is exactly the characteristic polynomial in (1.2.8).

Since ?(?) has the factor (? ? ?? )p , it follows that ?(D) has the factor (D ? ?? )p , so

the left side of (1.2.9) can be written in the form

g(D)(D ? ?? )p y = 0 ,

(1.2.10)

?

where g is a polynomial of degree n ? p. Now it’s quite easy to see that y = xk e? x satisfies

(1.2.10) (and therefore (1.2.7) also) for each k = 0, 1, . . . , p ? 1. Indeed, if we substitute this

function y into (1.2.10), we see that it is enough to show that

?

(D ? ?? )p (xk e? x ) = 0

k = 0, 1, . . . , p ? 1 .

(1.2.11)

10

Differential and Difference Equations

?

?

However, (D ? ?? )(xk e?? x ) = kxk?1 e? x , and if we apply (D ? ?? ) again,

?

?

(D ? ?? )2 (xk e?? x ) = k(k ? 1)xk?2 e? x ,

?

etc. Now since k < p it is clear that (D ? ?? )p (xk e?? x ) = 0, as claimed.

To summarize, then, if we encounter a root ?? of the characteristic equation, of multiplicity p, then corresponding to ?? we can find exactly p linearly independent solutions of

the differential equation, namely

?

?

?

?

e? x , xe? x , x2 e? x , . . . , xp?1 e? x .

(1.2.12)

Another way to state it is to say that the portion of the general solution of the given

differential equation that corresponds to a root ?? of the characteristic polynomial equation

?

is Q(x)e? x , where Q(x) is an arbitrary polynomial whose degree is one less than the

multiplicity of the root ?? .

One last mild complication may arise from roots of the characteristic equation that are

not real numbers. These don’t really require any special attention, but they do present a few

options. For instance, to solve y 00 + 4y = 0, we find the characteristic equation ?2 + 4 = 0,

and the complex roots ±2i. Hence the general solution is obtained by the usual rule as

y(x) = c1 e2ix + c2 e?2ix .

(1.2.13)

This is a perfectly acceptable form of the solution, but we could make it look a bit prettier

by using deMoivre’s theorem, which says that

e2ix = cos 2x + i sin 2x

e?2ix = cos 2x ? i sin 2x.

(1.2.14)

Then our general solution would look like

y(x) = (c1 + c2 ) cos 2x + (ic1 ? ic2 ) sin 2x.

(1.2.15)

But c1 and c2 are just arbitrary constants, hence so are c1 + c2 and ic1 ? ic2 , so we might

as well rename them c1 and c2 , in which case the solution would take the form

y(x) = c1 cos 2x + c2 sin 2x.

(1.2.16)

Here’s an example that shows the various possibilities:

y (8) ? 5y (7) + 17y (6) ? 997y (5) + 110y (4) ? 531y (3) + 765y (2) ? 567y 0 + 162y = 0. (1.2.17)

The equation was cooked up to have a characteristic polynomial that can be factored as

(? ? 2)(?2 + 9)2 (? ? 1)3 .

(1.2.18)

Hence the roots of the characteristic equation are 2 (simple), 3i (multiplicity 2), ?3i (multiplicity 2), and 1 (multiplicity 3).

1.3 Difference equations

11

Corresponding to the root 2, the general solution will contain the term c1 e2x . Corresponding to the double root at 3i we have terms (c2 + c3 x)e3ix in the solution. From the

double root at ?3i we get a contribution (c4 + c5 x)e?3ix , and finally from the triple root

at 1 we get (c6 + c7 x + c8 x2 )ex . The general solution is the sum of these eight terms.

Alternatively, we might have taken the four terms that come from 3i in the form

(c2 + c3 x) cos 3x + (c4 + c5 x) sin 3x.

(1.2.19)

Exercises 1.2

1. Obtain the general solutions of each of the following differential equations:

(a) y 00 + 5y 0 + 6y = 0

(b) y 00 ? 8y 0 + 7y = 0

(c) (D + 3)2 y = 0

(d) (D2 + 16)2 y = 0

(e) (D + 3)3 (D 2 ? 25)2 (D + 2)3 y = 0

2. Find a curve y = f (x) that passes through the origin with unit slope, and which

satisfies (D + 4)(D ? 1)y = 0.

1.3

Difference equations

Whereas a differential equation is an equation in an unknown function, a difference equation

is an equation in an unknown sequence. For example, suppose we know that a certain

sequence of numbers y0 , y1 , y2 , . . . satisfies the following conditions:

yn+2 + 5yn+1 + 6yn = 0

n = 0, 1, 2, . . .

(1.3.1)

and furthermore that y0 = 1 and y1 = 3.

Evidently, we can compute as many of the yn ’s as we need from (1.3.1), thus we would

get y2 = ?21, y3 = 87, y4 = ?309 so forth. The entire sequence of yn ’s is determined by

the difference equation (1.3.1) together with the two starting values.

Such equations are encountered when differential equations are solved on computers.

Naturally, the computer can provide the values of the unknown function only at a discrete

set of points. These values are computed by replacing the given differential equations by

a difference equation that approximates it, and then calculating successive approximate

values of the desired function from the difference equation.

Can we somehow “solve” a difference equation by obtaining a formula for the values

of the solution sequence? The answer is that we can, as long as the difference equation is

linear and has constant coefficients, as in (1.3.1). Just as in the case of differential equations

with constant coefficients, the correct strategy for solving them is to try a solution of the

12

Differential and Difference Equations

right form. In the previous section, the right form to try was y(x) = e?x . Now the winning

combination is y = ?n , where ? is a constant.

In fact, let’s substitute ?n for yn in (1.3.1) to see what happens. The left side becomes

?n+2 + 5?n+1 + 6?n = ?n (?2 + 5? + 6) = 0.

(1.3.2)

Just as we were able to cancel the common factor e?x in the differential equation case, so

here we can cancel the ?n , and we’re left with the quadratic equation

?2 + 5? + 6 = 0.

(1.3.3)

The two roots of this characteristic equation are ? = ?2 and ? = ?3. Therefore the

sequence (?2)n satisfies (1.3.1) and so does (?3)n . Since the difference equation is linear,

it follows that

yn = c1 (?2)n + c2 (?3)n

(1.3.4)

is also a solution, whatever the values of the constants c1 and c2 .

Now it is evident from (1.3.1) itself that the numbers yn are uniquely determined if we

prescribe the values of just two of them. Hence, it is very clear that when we have a solution

that contains two arbitrary constants we have the most general solution.

When we take account of the given data y0 = 1 and y1 = 3, we get the two equations

(

1 = c1 + c2

3 = (?2)c1 + (?3)c2

(1.3.5)

from which c1 = 6 and c2 = ?5. Finally, we use these values of c1 and c2 in (1.3.4) to get

yn = 6(?2)n ? 5(?3)n

n = 0, 1, 2, . . . .

(1.3.6)

Equation (1.3.6) is the desired formula that represents the unique solution of the given

difference equation together with the prescribed starting values.

Let’s step back a few paces to get a better view of the solution. Notice that the formula

(1.3.6) expresses the solution as a linear combination of nth powers of the roots of the

associated characteristic equation (1.3.3). When n is very large, is the number yn a large

number or a small one? Evidently the powers of ?3 overwhelm those of ?2, so the sequence

will behave roughly like a constant times powers of ?3. This means that we should expect

the members of the sequence to alternate in sign and to grow rapidly in magnitude.

So much for the equation (1.3.1). Now let’s look at the general case, in the form of a

linear difference equation of order p:

yn+p + a1 yn+p?1 + a2 yn+p?2 + · · · + ap yn = 0.

(1.3.7)

We try a solution of the form yn = ?n , and after substituting and canceling, we get the

characteristic equation

?p + a1 ?p?1 + a2 ?p?2 + · · · + ap = 0.

(1.3.8)

1.3 Difference equations

13

This is a polynomial equation of degree p, so it has p roots, counting multiplicities, somewhere in the complex plane.

Let ?? be one of these p roots. If ?? is simple (i.e., has multiplicity 1) then the part

of the general solution that corresponds to ?? is c(?? )n . If, however, ?? is a root of

multiplicity k > 1 then we must multiply the solution c(?? )n by an arbitrary polynomial

in n, of degree k ? 1, just as in the corresponding case for differential equations we used an

arbitrary polynomial in x of degree k ? 1.

We illustrate this, as well as the case of complex roots, by considering the following

difference equation of order five:

yn+5 ? 5yn+4 + 9yn+3 ? 9yn+2 + 8yn+1 ? 4yn = 0.

(1.3.9)

This example is rigged so that the characteristic equation can be factored as

(?2 + 1)(? ? 2)2 (? ? 1) = 0

(1.3.10)

from which the roots are obviously i, ?i, 2 (multiplicity 2), 1.

Corresponding to the roots i, ?i, the portion of the general solution is c1 in + c2 (?i)n .

Since

n?

n?

in = ein?/2 = cos

+ i sin

(1.3.11)

2

2

and similarly for (?i)n , we can also take this part of the general solution in the form

n?

c1 cos

2

n?

+ c2 sin

.

2

(1.3.12)

The double root ? = 2 contributes (c3 + c4 n)2n , and the simple root ? = 1 adds c5 to

the general solution, which in its full glory is

yn = c1 cos

n?

2

+ c2 sin

n?

2

+ (c3 + c4 n)2n + c5 .

(1.3.13)

The five constants would be determined by prescribing five initial values, say y0 , y1 , y2 , y3

and y4 , as we would expect for the equation (1.3.9).

Exercises 1.3

1. Obtain the general solution of each of the following difference equations:

(a) yn+1 = 3yn

(b) yn+1 = 3yn + 2

(c) yn+2 ? 2yn+1 + yn = 0

(d) yn+2 ? 8yn+1 + 12yn = 0

(e) yn+2 ? 6yn+1 + 9yn = 1

(f) yn+2 + yn = 0

14

Differential and Difference Equations

2. Find the solution of the given difference equation that takes the prescribed initial

values:

(a)

(b)

(c)

(d)

yn+2 = 2yn+1 + yn ; y0 = 0; y1 = 1

yn+1 = ?yn + ?; y0 = 1

yn+4 + yn = 0; y0 = 1; y1 = ?1; y2 = 1; y3 = ?1

yn+2 ? 5yn+1 + 6yn = 0; y0 = 1; y1 = 2

3. (a) For each of the difference equations in problems 1 and 2, evaluate

yn+1

lim

n?? yn

(1.3.14)

if it exists.

(b) Formulate and prove a general theorem about the existence of, and value of the

limit in part (a) for a linear difference equation with constant coefficients.

(c) Reverse the process: given a polynomial equation, find its root of largest absolute

value by computing from a certain difference equation and evaluating the ratios

of consecutive terms.

(d) Write a computer program to implement the method in part (c). Use it to

calculate the largest root of the equation

x8 = x7 + x6 + x5 + · · · + 1.

1.4

(1.3.15)

Computing with difference equations

This is, after all, a book about computing, so let’s begin with computing from difference

equations since they will give us a chance to discuss some important questions that concern

the design of computer programs. For a sample difference equation we’ll use

yn+3 = yn+2 + 5yn+1 + 3yn

(1.4.1)

together with the starting values y0 = y1 = y2 = 1. The reader might want, just for practice,

to find an explicit formula for this sequence by the methods of the previous section.

Suppose we want to compute a large number of these y’s in order to verify some property

that they have, for instance to check that

yn+1

lim

=3

(1.4.2)

n?? yn

which must be true since 3 is the root of largest absolute value of the characteristic equation.

As a first approach, we might declare y to be a linear array of some size large enough to

accommodate the expected length of the calculation. Then the rest is easy. For each n, we

would calculate the next yn+1 from (1.4.1), we would divide it by its predecessor yn to get

a new ratio. If the new ratio agrees sufficiently well with the previous ratio we announce

that the computation has terminated and print the new ratio as our answer. Otherwise, we

move the new ratio to the location of the old ratio, increase n and try again.

If we were to write this out as formal procedure (algorithm) it might look like:

1.4 Computing with difference equations

15

y0 := 1; y1 := 1; y2 := 1; n := 2;

newrat := ?10; oldrat := 1;

while |newrat ? oldrat| ? 0.000001 do

oldrat := newrat; n := n + 1;

yn := yn?1 + 5 ? yn?2 + 3 ? yn?3 ;

newrat := yn /yn?1

endwhile

print newrat; Halt.

We’ll use the symbol ‘:=’ to mean that we are to compute the quantity on the right, if

necessary, and then store it in the place named on the left. It can be read as ‘is replaced

by’ or ‘is assigned.’ Also, the block that begins with ‘while’ and ends with ‘endwhile’

represents a group of instructions that are to be executed repeatedly until the condition

that follows ‘while’ becomes false, at which point the line following ‘endwhile’ is executed.

The procedure just described is fast, but it uses lots of storage. If, for instance, such a

program needed to calculate 79 y’s before convergence occurred, then it would have used

79 locations of array storage. In fact, the problem above doesn’t need that many locations

because convergence happens a lot sooner. Suppose you wanted to find out how much

sooner, given only a programmable hand calculator with ten or twenty memory locations.

Then you might appreciate a calculation procedure that needs just four locations to hold

all necessary y’s.

That’s fairly easy to accomplish, though. At any given moment in the program, what

we need to find the next y are just the previous three y’s. So why not save only those three?

We’ll use the previous three to calculate the next one, and stow it for a moment in a fourth

location. Then we’ll compute the new ratio and compare it with the old. If they’re not

close enough, we move each one of the three newest y’s back one step into the places where

we store the latest three y’s and repeat the process. Formally, it might be:

y := 1; ym1 := 1; ym2 := 1;

newrat := ?10; oldrat := 1;

while |newrat ? oldrat| ? 0.000001 do

ym3 := ym2; ym2 := ym1; ym1 := y;

oldrat := newrat;

y := ym1 + 5 ? ym2 + 3 ? ym3;

newrat := y/ym1 endwhile;

print newrat; Halt.

The calculation can now be done in exactly six memory locations (y, ym1, ym2, ym3,

oldrat, newrat) no matter how many y’s have to be calculated, so you can undertake it

on your hand calculator with complete confidence. The price that we pay for the memory

saving is that we must move the data around a bit more.

One should not think that such programming methods are only for hand calculators. As

we progress through the numerical solution of differential equations we will see situations

in which each of the quantities that appears in the difference equation will itself be an

16

Differential and Difference Equations

array (!), and that very large numbers, perhaps thousands, of these arrays will need to be

computed. Even large computers might quake at the thought of using the first method

above, rather than the second, for doing the calculation. Fortunately, it will almost never

be necessary to save in memory all of the computed values simultaneously. Normally, they

will be computed, and then printed or plotted, and never needed except in the calculation

of their immediate successors.

Exercises 1.4

1. The Fibonacci numbers F0 , F1 , F2 , . . . are defined by the recurrence formula Fn+2 =

Fn+1 + Fn for n = 0, 1, 2, . . . together with the starting values F0 = 0, F1 = 1.

(a) Write out the first ten Fibonacci numbers.

(b) Derive an explicit formula for the nth Fibonacci number Fn .

(c) Evaluate your formula for n = 0, 1, 2, 3, 4.

(d) Prove directly from your formula that the Fibonacci numbers are integers (This is

perfectly obvious from their definition, but is not so obvious from the formula!).

(e) Evaluate

lim

n??

Fn+1

Fn

(1.4.3)

(f) Write a computer program that will compute Fibonacci numbers and print out

the limit in part (e) above, correct to six decimal places.

(g) Write a computer program that will compute the first?

40 members of the modified

Fibonacci sequence in which F0 = 1 and F1 = (1 ? 5)/2. Do these computed

numbers seem to be approaching zero? Explain carefully what you see and why

it happens.

(h) Modify the program of part (h) to run in higher (or double) precision arithmetic.

Does it change any of your answers?

2. Find the most general solution of each of the following difference equations:

(a) yn+1 ? 2yn + yn?1 = 0

(b) yn+1 = 2yn

(c) yn+2 + yn = 0

(d) yn+2 + 3yn+1 + 3yn + yn?1 = 0

1.5

Stability theory

In the study of natural phenomena it is most often true that a small change in conditions

will produce just a small change in the state of the system being studied. If, for example,

a very slight increase in atmospheric pollution could produce dramatically large changes

in populations of flora and fauna, or if tiny variations in the period of the earth’s rotation

1.5 Stability theory

17

produced huge changes in climatic conditions, the world would be a very different place to

live in, or to try to live in. In brief, we may say that most aspects of nature are stable.

When physical scientists attempt to understand some facet of nature, they often will

make a mathematical model. This model will usually not faithfully reproduce all of the

structure of the original phenomenon, but one hopes that the important features of the

system will be preserved in the model, so that predictions will be possible. One of the most

important features to preserve is that of stability.

For instance, the example of atmospheric pollution and its effect on living things referred

to above is important and very complex. Therefore considerable effort has gone into the

construction of mathematical models that will allow computer studies of the effects of

atmospheric changes. One of the first tests to which such a model should be subjected is

that of stability: does it faithfully reproduce the observed fact that small changes produce

small changes? What is true in nature need not be true in a man-made model that is a

simplification or idealization of the real world.

Now suppose that we have gotten ourselves over this hurdle, and we have constructed

a model that is indeed stable. The next step might be to go to the computer and do

calculations from the model, and to use these calculations for predicting the effects of

various proposed actions. Unfortunately, yet another layer of approximation is usually

introduced at this stage, because the model, even though it is a simplification of the real

world, might still be too complicated to solve exactly on a computer.

For instance, may models use differential equations. Models of the weather, of the motion of fluids, of the movement of astronomical objects, of spacecraft, of population growth,

of predator-prey relationships, of electric circuit transients, and so forth, all involve differential equations. Digital computers solve differential equations by approximating them by

difference equations, and then solving the difference equations. Even though the differential

equation that represents our model is indeed stable, it may be that the difference equation

that we use on the computer is no longer stable, and that small changes in initial data on

the computer, or small roundoff errors, will produce not small but very large changes in the

computed solution.

An important job of the numerical analyst is to make sure that this does not happen,

and we will find that this theme of stability recurs throughout our study of computer

approximations.

As an example of instability in differential equations, suppose that some model of a

system led us to the equation

y 00 ? y 0 ? 2y = 0

(1.5.1)

together with the initial data

y(0) = 1;

y 0 (0) = ?1.

(1.5.2)

We are thinking of the independent variable t as the time, so we will be interested in

the solution as t becomes large and positive.

The general solution of (1.5.1) is y(t) = c1 e?t + c2 e2t . The initial conditions tell us that

c1 = 1 and c2 = 0, hence the solution of our problem is y(t) = e?t , and it represents a

18

Differential and Difference Equations

function that decays rapidly to zero with increasing t. In fact, when t = 10, the solution

has the value 0.000045.

Now let’s change the initial data (1.5.2) just a bit, by asking for a solution with y 0 (0) =

?0.999. It’s easy to check that the solution is now

y(t) = (0.999666 . . .)e?t + (0.000333 . . .)e2t

(1.5.3)

instead of just y(t) = e?t . If we want the value of the solution at t = 10, we would find

that it has changed from 0.000045 to about 7.34.

At t = 20 the change is even more impressive, from 0.00000002 to 161, 720+, just from

changing the initial value of y 0 from ?1 to ?0.999. Let’s hope that there are no phenomena

in nature that behave in this way, or our lives hang by a slender thread indeed!

Now exactly what is the reason for the observed instability of the equation (1.5.1)?

The general solution of the equation contains a falling exponential term c1 e?t , and a rising

exponential term c2 e2t . By prescribing the initial data (1.5.2) we suppressed the growing

term, and picked out only the decreasing one. A small change in the initial data, however,

results in the presence of both terms in the solution.

Now it’s time for a formal

Definition: A differential equation is said to be stable if for every set of initial data (at

t = 0) the solution of the differential equation remains bounded as t approaches infinity.

A differential equation is called strongly stable if, for every set of initial data (at t = 0)

the solution not only remains bounded, but approaches zero as t approaches infinity.

What makes the equation (1.5.1) unstable, then, is the presence of a rising exponential

in its general solution. In other words, if we have a differential equation whose general

solution contains a term e?t in which ? is positive, that equation is unstable.

Let’s restrict attention now to linear differential equations with constant coefficients.

We know from section 1.2 that the general solution of such an equation is a sum of terms

of the form

(polynomial in t)e?t .

(1.5.4)

Under what circumstances does such a term remain bounded as t becomes large and positive?

Certainly if ? is negative then the term stays bounded. Likewise, if ? is a complex

number and its real part is negative, then the term remains bounded. If ? has positive real

part the term is unbounded.

This takes care of all of the possibilities except the case where ? is zero, or more generally,

the complex number ? has zero real part (is purely imaginary). In that case the question

of whether (polynomial in t)e?t remains bounded depend on whether the “polynomial in t”

is of degree zero (a constant polynomial) or of higher degree. If the polynomial is constant

then the term does indeed remain bounded for large positive t, whereas otherwise the term

will grow as t gets large, for some values of the initial conditions, thereby violating the

definition of stability.

1.6 Stability theory of difference equations

19

Now recall that the “polynomial in t” is in fact a constant if the root ? is a simple

root of the characteristic equation of the differential equation, and otherwise it is of higher

degree. This observation completes the proof of the following:

Theorem 1.5.1 A linear differential equation with constant coefficients is stable if and

only if all of the roots of its characteristic equation lie in the left half plane, and those that

lie on the imaginary axis, if any, are simple. Such an equation is strongly stable if and only

if all of the roots of its characteristic equation lie in the left half plane, and none lie on the

imaginary axis.

Exercises 1.5

1. Determine for each of the following differential equations whether it is strongly stable,

stable, or unstable.

(a) y 00 ? 5y 0 + 6y = 0

(b) y 00 + 5y 0 + 6y = 0

(c) y 00 + 3y = 0

(d) (D + 3)3 (D + 1)y = 0

(e) (D + 1)2 (D 2 + 1)2 y = 0

(f) (D4 + 1)y = 0

2. Make a list of some natural phenomena that you think are unstable. Discuss.

3. The differential equation y 00 ?y = 0 is to be solved with the initial conditions y(0) = 1,

y 0 (0) = ?1, and then solved again with y(0) = 1, y 0 (0) = ?0.99. Compare the two

solutions when x = 20.

4. For exactly which real values of the parameter ? is each of the following differential

equations stable? . . . strongly stable?

(a) y 00 + (2 + ?)y 0 + y = 0

(b) y 00 + ?y 0 + y = 0

(c) y 0 + ?y = 1

1.6

Stability theory of difference equations

In the previous section we discussed the stability of differential equations. The key ideas

were that such an equation is stable if every one of its solutions remains bounded as t

approaches infinity, and strongly stable if the solutions actually approach zero.

Similar considerations apply to difference equations, and for similar reasons. As an

example, take the equation

5

yn+1 = yn ? yn?1

2

(n ? 1)

(1.6.1)

20

Differential and Difference Equations

along with the initial equations

y0 = 1;

y1 = 0.5 .

(1.6.2)

It’s easy to see that the solution is yn = 2?n , and of course, this is a function that

rapidly approaches zero with increasing n.

Now let’s change the initial data (1.6.2), say to

y0 = 1;

y1 = 0.50000001

(1.6.3)

instead of (1.6.2).

The solution of the difference equation with these new data is

y = (0.0000000066 . . .)2n + (0.9999999933 . . .)2?n .

(1.6.4)

The point is that the coefficient of the growing term 2n is small, but 2n grows so fast that

after a while the first term in (1.6.4) will be dominant. For example, when n = 30, the

solution is y30 = 7.16, compared to the value y30 = 0.0000000009 of the solution with the

original initial data (1.6.2). A change of one part in fifty million in the initial condition

produced, thirty steps later, an answer one billion times as large.

The fault lies with the difference equation, because it has both rising and falling components to its general solution. It should be clear that it is hopeless to do extended computation with an unstable difference equation, since a small roundoff error may alter the

solution beyond recognition several steps later.

As in the case of differential equations, we’ll say that a difference equation is stable

if every solution remains bounded as n grows large, and that it is strongly stable if every

solution approaches zero as n grows large. Again, we emphasize that every solution must

be well behaved, not just the solution that is picked out by a certain set of initial data. In

other words, the stability, or lack of it, is a property of the equation and not of the starting

values.

Now consider the case where the difference equation is linear with constant coefficients.

The we know that the general solution is a sum of terms of the form

(polynomial in n)?n .

(1.6.5)

Under what circumstances will such a term remain bounded or approach zero?

Suppose |?| < 1. Then the powers of ? approach zero, and multiplication by a polynomial in n does not alter that conclusion. Suppose |?| > 1. Then the sequence of powers

grows unboundedly, and multiplication by a nonzero polynomial only speeds the parting

guest.

Finally suppose the complex number ? has absolute value 1. Then the sequence of its

powers remains bounded (in fact they all have absolute value 1), but if we multiply by a

nonconstant polynomial, the resulting expression would grow without bound.

To summarize then, the term (1.6.5), if the polynomial is not identically zero, approaches

zero with increasing n if and only if |?| < 1. It remains bounded as n increases if and only

1.6 Stability theory of difference equations

21

if either (a) |?| < 1 or (b) |?| = 1 and the polynomial is of degree zero (a constant). Now

we have proved:

Theorem 1.6.1 A linear difference equation with constant coefficients is stable if and only

if all of the roots of its characteristic equation have absolute value at most 1, and those of

absolute value 1 are simple. The equation is strongly stable if and only if all of the roots

have absolute value strictly less than 1.

Exercises 1.6

1. Determine, for each of the following difference equations whether it is strongly stable,

stable, or unstable.

(a) yn+2 ? 5yn+1 + 6yn = 0

(b) 8yn+2 + 2yn+1 ? 3yn = 0

(c) 3yn+2 + yn = 0

(d) 3yn+3 + 9yn+2 ? yn+1 ? 3yn = 0

(e) 4yn+4 + 5yn+2 + yn = 0

2. The difference equation 2yn+2 + 3yn+1 ? 2yn = 0 is to be solved with the initial

conditions y0 = 2, y1 = 1, and then solved again with y0 = 2, y1 = 0.99. Compare y20

for the two solutions.

3. For exactly which real values of the parameter ? is each of the following difference

equations stable? . . . strongly stable?

(a) yn+2 + ?yn+1 + yn = 0

(b) yn+1 + ?yn = 1

(c) yn+2 + yn+1 + ?yn = 0

4. (a) Consider the (constant-coefficient) difference equation

a0 yn+p + a1 yn+p?1 + a2 yn+p?2 + · · · + ap yn = 0.

(1.6.6)

Show that this difference equation cannot be stable if |ap /a0 | > 1.

(b) Give an example to show that the converse of the statement in part (a) is false.

Namely, exhibit a difference equation for which |ap /a0 | < 1 but the equation is unstable anyway.

22

Differential and Difference Equations

Chapter 2

The Numerical Solution of

Differential Equations

2.1

Euler’s method

Our study of numerical methods will begin with a very simple procedure, due to Euler.

We will state it as a method for solving a single differential equation of first order. One of

the nice features of the subject of numerical integration of differential equations is that the

techniques that are developed for just one first order differential equation will apply, with

very little change, both to systems of simultaneous first order equations and to equations of

higher order. Hence the consideration of a single equation of first order, seemingly a very

special case, turns out to be quite general.

By an initial-value problem we mean a differential equation together with enough given

values of the unknown function and its derivatives at an initial point x0 to determine the

solution uniquely.

Let’s suppose that we are given an initial-value problem of the form

y 0 = f (x, y);

y(x0 ) = y0 .

(2.1.1)

Our job is to find numerical approximate values of the unknown function y at points x

to the right of (larger than) x0 .

What we actually will find will be approximate values of the unknown function at a

discrete set of points x0 , x1 = x0 + h, x2 = x0 + 2h, x3 = x0 + 3h, etc. At each of these

points xn we will compute yn , our approximation to y(xn ).

Hence, suppose that the spacing h between consecutive points has been chosen. We

propose to start at the point x0 where the initial data are given, and move to the right,

obtaining y1 from y0 , then y2 from y1 and so forth until sufficiently many values have been

found.

Next we need to derive a method by which each value of y can be obtained from its

immediate predecessor. Consider the Taylor series expansion of the unknown function y(x)

24

The Numerical Solution of Differential Equations

about the point xn

y 00 (X)

,

(2.1.2)

2

where we have halted the expansion after the first power of h and in the remainder term,

the point X lies between xn and xn + h.

y(xn + h) = y(xn ) + hy 0 (xn ) + h2

Now equation (2.1.2) is exact, but of course it cannot be used for computation because

the point X is unknown. On the other hand, if we simply “forget” the error term, we’ll

have only an approximate relation instead of an exact one, with the consolation that we

will be able to compute from it. The approximate relation is

y(xn + h) ? y(xn ) + hy 0 (xn ).

(2.1.3)

Next define yn+1 to be the approximate value of y(xn+1 ) that we obtain by using the

right side of (2.1.3) instead of (2.1.2). Then we get

yn+1 = yn + hyn0 .

(2.1.4)

Now we have a computable formula for the approximate values of the unknown function,

because the quantity yn0 can be found from the differential equation (2.1.1) by writing

yn0 = f (xx , yn ),

(2.1.5)

and if we do so then (2.1.4) takes the form

yn+1 = yn + hf (xn , yn ).

(2.1.6)

This is Euler’s method, in a very explicit form, so that the computational procedure is

clear. Equation (2.1.6) is in fact a recurrence relation, or difference equation, whereby each

value of yn is computed from its immediate predecessor.

Let’s use Euler’s method to obtain a numerical solution of the differential equation

y 0 = 0.5y

(2.1.7)

together with the starting value y(0) = 1. The exact solution of this initial-value problem

is obviously y(x) = ex/2 .

Concerning the approximate solution by Euler’s method, we have, by comparing (2.1.7)

with (2.1.1), f (x, y) = 0.5y, so

yn

2

h

= 1+

yn .

2

yn+1 = yn + h

(2.1.8)

Therefore, in this example, each yn will be obtained from its predecessor by multiplication

by 1 + h2 . To be quite specific, let’s take h to be 0.05. Then we show below, for each value

of x = 0, 0.05, 0.10, 0.15, 0.20, . . . the approximate value of y computed from (2.1.8) and

the exact value y(xn ) = exn /2 :

2.1 Euler’s method

25

x

0.00

0.05

0.10

0.15

0.20

0.25

..

.

Euler(x)

1.00000

1.02500

1.05063

1.07689

1.10381

1.13141

..

.

Exact(x)

1.00000

1.02532

1.05127

1.07788

1.10517

1.13315

..

.

1.00

2.00

3.00

..

.

1.63862

2.68506

4.39979

..

.

1.64872

2.71828

4.48169

..

.

5.00

10.00

11.81372

139.56389

12.18249

148.41316

table 1

Considering the extreme simplicity of the approximation, it seems that we have done

pretty well by this equation. Let’s continue with this simple example by asking for a formula

for the numbers that are called Euler(x) in the above table. In other words, exactly what

function of x is Euler(x)?

To answer this, we note first that each computed value yn+1 is obtained according to

(2.1.8) by multiplying its predecessor yn by 1 + h2 . Since y0 = 1, it is clear tha we will

compute yn = (1 + h2 )n . Now we want to express this in terms of x rather than n. Since

xn = nh, we have n = x/h, and since h = 0.05 we have n = 20x. Hence the computed

approximation to y at a particular point x will be 1.02520x , or equivalently

Euler(x) = (1.638616 . . .)x .

(2.1.9)

The approximate values can now easily be compared with the true solution, since

x

1

Exact(x) = e 2 = e 2

x

= (1.648721 . . .)x .

(2.1.10)

Therefore both the exact solution of this differential equation and its computed solution

have the form (const.)x . The correct value of “const.” is e1/2 , and the value that is, in

effect, used by Euler’s method is (1 + h2 )1/h . For a fixed value of x, we see that if we use

Euler’s method with smaller and smaller values of h (neglecting the increase in roundoff

error that is sure to result), the values Euler(x) will converge to Exact(x), because

h

lim 1 +

h?0

2

Exercises 2.1

1/h

1

= e2 .

(2.1.11)

26

The Numerical Solution of Differential Equations

1. Verify the limit (2.1.11).

2. Use a calculator or a computer to integrate each of the following differential equations

forward ten steps, using a spacing h = 0.05 with Euler’s method. Also tabulate the

exact solution at each value of x that occurs.

(a) y 0 (x) = xy(x); y(0) = 1

(b) y 0 (x) = xy(x) + 2; y(0) = 1

y(x)

(c) y 0 (x) =

; y(0) = 1

1+x

(d) y 0 (x) = ?2xy(x)2 ; y(0) = 10

2.2

Software notes

One of the main themes of our study will be the preparation of programs that not only

work, but also are easily readable and useable by other people. The act of communication

that must take place before a program can be used by persons other than its author is a

difficult one to carry out, and we will return several times to the principles that have evolved

as guides to the preparation of readable software. Here are some of these guidelines.

1. Documentation

The documentation of a program is the set of written instructions to a user that inform

the user about the purpose and operation of the program. At the moment that the job

of writing and testing a program has been completed it is only natural to feel an urge to

get the whole thing over with and get on to the next job. Besides, one might think, it’s

perfectly obvious how to use this program. Some programs may be obscure, but not this

one.

It is amazing how rapidly our knowledge of our very own program fades. If we come

back to a program after a lapse of a few months’ time, it often happens that we will have no

idea what the program did or how to use it, at least not without making a large investment

of time.

For that reason it is important that when our program has been written and tested it

should be documented immediately, while our memory of it is still green. Furthermore, the

best place for documentation is in the program itself, in “comment” statements. That way

one can be sure that when the comments are needed they will be available.

The first mission of program documentation is to describe the purpose of the program.

State clearly the problem that the program solves, or the exact operation that it performs

on its input in order to get its output.

Already in this first mission, a good bit of technical skill can be brought to bear that

will be very helpful to the use, by intertwining the description of the program purpose with

the names of the communicating variables in the program.

Let’s see what that means by considering an example. Suppose we have written a

subroutine that searches through a specified row of a matrix to find the element of largest

2.2 Software notes

27

absolute value, and outputs a column in which it was found. Such a routine, in Maple for

instance, might look like this:

search:=proc(A,i)

local j, winner, jwin;

winner:=-1;

for j from 1 to coldim(A) do

if (abs(A[i,j])>winner) then

winner:=abs(A[i,j]) ; jwin:=j fi

od;

return(jwin);

end;

Now let’s try our hand at documenting this program:

“The purpose of this program is to search a given row of a matrix to find an

element of largest absolute value and return the column in which it was found.”

That is pretty good documentation, perhaps better than many programs get. But we

can make it a lot more useful by doing the intertwining that we referred to above. There

we said that the description should be related to the communicating variables. Those

variables are the ones that the user can see. They are the input and output variables of

the subroutine. In most important computer languages, the communicating variables are

announced in the first line of the coding of a procedure or subroutine. Maple follows this

convention at least for the input variables, although the output variable is usually specified

in the “return” statement.

In the first line of the little subroutine above we see the list (A, i) of its input variables

(or “arguments”). These are the ones that the user has to understand, as opposed to the

other “local” variables that live inside the subroutine but don’t communicate with the

outside world (like j, winner, jwin, which are listed on the second line of the program).

The best way to help the user to understand these variables is to relate them directly

to the description of the purpose of the program.

“The purpose of this program is to search row I of a given matrix A to find an

entry of largest absolute value, and returns the column jwin where that entry

lives.”

We’ll come back to the subject of documentation, but now let’s mention another ingredient of ease of use of programs, and that is:

2. Modularity

It is important to divide a long program into a number of smaller modules, each with a

clearly stated set of inputs and outputs, and each with its own documentation. That means

that we should get into the habit of writing lots of subroutines or procedures, because

28

The Numerical Solution of Differential Equations

the subroutine or procedure mode of expression forces one to be quite explicit about the

relationship of the block of coding to the rest of the world.

When we are writing a large program we would all write a subroutine if we found that

a certain sequence of steps was being called for repeatedly. Beyond this, however, there are

numerous inducements for breaking off subroutines even if the block of coding occurs just

once in the main program.

For one thing it’s easier to check out the program. The testing procedure would consist

of first testing each of the subroutines separately on test problems designed just for them.

Once the subroutines work, it would remain only to test their relationships to the calling

program.

For another reason, we might discover a better, faster, more elegant, or what-have-you

method of performing the task that one of these subroutines does. Then we would be able

to yank out the former subroutine and plug in the new one, while being careful only to make

sure that the new subroutine relates to the same inputs and outputs as the old one. If jobs

within a large program are not broken into subroutines it can be much harder to isolate the

block of coding that deals with a particular function and remove it without affecting the

whole works.

For another reason, if one be needed, it may well happen that even though the job that

is done by the subroutine occurs only once in the current program, it may recur in other

programs as yet undreamed of. If one is in the habit of writing small independent modules

and stringing them together to make large programs, then it doesn’t take long before one

has a library of useful subroutines, each one tested, working and documented, that will

greatly simplify the task of writing future programs.

Finally, the practice of subdividing the large jobs into the smaller jobs of which they are

composed is an extremely valuable analytical skill, one that is useful not only in programming, but in all sorts of organizational activities where smaller efforts are to be pooled in

order to produce a larger effect. It is therefore a quality of mind that provides much of its

own justification.

In this book, the major programs that are the objects of study have been broken up into

subroutines in the expectation that the reader will be able to start writing and checking out

these modules even before the main ideas of the current subject have been fully explained.

This was done in part because some of these programs are quite complex, and it would be

unreasonable to expect the whole program to be written in a short time. It was also done

to give examples of the process of subdivision that we have been talking about.

For instance, the general linear algebra program for solving systems of linear simultaneous equations in Chapter 3, has been divided into six modules, and they are described in

section 3.3. The reader might wish to look ahead at those routines and to verify that even

though their relationship to the whole job of solving equations is by no means clear now,

nonetheless, because of the fact that they are independent and self-contained, they can be

programmed and checked out at any time without waiting for the full explanation.

One more ingredient that is needed for the production of useful software is:

2.3 Systems and equations of higher order

29

3. Style

We don’t mean style in the sense of “class,” although this is as welcome in programming

as it is elsewhere. There have evolved a number of elements of good programming style,

and these will mainly be discussed as they arise. But two of them (one trivial and one quite

deep) are:

(a) Indentation: The instructions that lie within the range of a loop are indented in the

program listing further to the right than the instructions that announce that the loop

is about to begin, or that it has just terminated.

(b) Top-down structuring: When we visualize the overall logical structure of a complicated program we see a grand loop, within which there are several other loops and

branchings, within which . . . etc. According to the principles of top-down design the

looping and branching structure of the program should be visible at once in the listing. That is, we should see an announcement of the opening of the grand loop, then

indented under that perhaps a two-way branch (if-then-else), where, under the “then”

one sees all that will happen if the condition is met, and under the “else” one sees

what happens if it is not met.

When we say that we see all that will happen, we mean that there are not any “go-to”

instructions that would take our eye out of the flow of the if-then-else loop to some

other page. It all happens right there on the same page, under “then” and under

“else”.

These few words can scarcely convey the ideas of structuring, which we leave to the

numerous examples in the sequel.

2.3

Systems and equations of higher order

We have already remarked that the methods of numerical integration for a single firstorder differential equation carry over with very little change to systems of simultaneous

differential equations of first order. In this section we’ll discuss exactly how this is done,

and furthermore, how the same idea can be applied to equations of higher order than the

first. Euler’s method will be used as the example, but the same transformations will apply

to all of the methods that we will study.

In Euler’s method for one equation, the approximate value of the unknown function at

the next point xn+1 = xn + h is calculated from

yn+1 = yn + hf (xn , yn ).

(2.3.1)

Now suppose that we are trying to solve not just a single equation, but a system of N

simultaneous equations, say

yi0 (x) = fi (x, y1 , y2 , . . . yN )

i = 1, . . . , N.

(2.3.2)

30

The Numerical Solution of Differential Equations

Equation (2.3.2) represents N equations, in each of which just one derivative appears,

and whose right-hand side may depend on x, and on all of the unknown functions, but not

on their derivatives. The “fi ” indicates that, of course, each equation can have a different

right-hand side.

Now introduce the vector y(x) of unknown functions

y(x) = [y1 (x), y2 (x), y3 (x), . . . , yN (x)]

(2.3.3)

and the vector f = f (x, y(x)) of right-hand sides

f = [f1 (x, y), f2 (x, y), . . . , fN (x, y)].

(2.3.4)

In terms of these vectors, equation (2.3.2) can be rewritten as

y 0 (x) = f (x, y(x)).

(2.3.5)

We observe that equation (2.3.5) looks just like our standard form (2.1.1) for a single

equation in a single unknown function, except for the bold face type, i.e., except for the

fact that y and f now represent vector quantities.

To apply a numerical method such as that of Euler, then, all we need to do is to take

the statement of the method for a single differential equation in a single unknown function,

and replace y(x) and f (x, y(x)) by vector quantities as above. We will then have obtained

the generalization of the numerical method to systems.

To be specific, Euler’s method for a single equation is

yn+1 = yn + hf (x, yn )

(2.3.6)

so Euler’s method for a system of differential equations will be

y n+1 = y n + hf (xn , y n ).

(2.3.7)

This means that if we know the entire vector y of unknown functions at the point x = xn ,

then we can find the entire vector of unknown functions at the next point xn+1 = xn + h

by means of (2.3.7).

In detail, if y i (xn ) denotes the computed approximate value of the unknown function yi

at the point xn , then what we must calculate are

y i (xn+1 ) = y i (xn ) + hfi (xn , y 1 (xn ), y 2 (xn ), . . . , y N (xn ))

(2.3.8)

for each i = 1, 2, . . . , N .

As an example, take the pair of differential equations

(

y10 = x + y1 + y2

y20 = y1 y2 + 1

together with the initial values y1 (0) = 0, y2 (0) = 1.

(2.3.9)

2.3 Systems and equations of higher order

31

Now the vector of unknown functions is y = [y1 , y2 ], and the vector of right-hand sides

is f = [x + y1 + y2 , y1 y2 + 1]. Initially, the vector of unknowns is y = [0, 1]. Let’s choose a

step size h = 0.05. Then we calculate

"

and

"

y 1 (0.10)

y 2 (0.10)

y1 (0.05)

y2 (0.05)

#

"

=

0.05

1.05

#

"

=

0

1

#

#

+ 0.05

"

+ 0.05

"

0+0+1

0?1+1

0.05 + 0.05 + 1.05

0.05 ? 1.05 + 1

#

"

=

#

"

=

0.05

1.05

#

0.1075

1.102625

(2.3.10)

#

(2.3.11)

and so forth. At each step we compute the vector of approximate values of the two unknown

functions from the corresponding vector at the immediately preceding step.

Let’s consider the preparation of a computer program that will carry out the solution,

by Euler’s method, of a system of N simultaneous equations of the form (2.3.2), which we

will rewrite just slightly, in the form

yi0 = fi (x, y)

i = 1, . . . , N.

(2.3.12)

Note that on the left is just one of the unknown functions, and on the right there may

appear all N of them in each equation.

Evidently we will need an array Y of length N to hold the values of the N unknown

functions at the current point x. Suppose we have computed the array Y at a point x, and

we want to get the new array Y at the point x + h. Exactly what do we do?

Some care is necessary in answering this question because there is a bit of a snare in

the underbrush. The new values of the N unknown functions are calculated from (2.3.8) or

(2.3.12) in a certain order. For instance, we might calculate y 1 (x + h), then y 2 (x + h), then

y 3 (x + h),. . . , then yN (x + h).

The question is this: when we compute y 1 (x+h), where shall we put it? If we put it into

Y[1], the first position of the Y array in storage, then the previous contents of Y[1] are

lost, i.e., the value of y1 (x) is lost. But we aren’t finished with y 1 (x) yet; it’s still needed to

compute y 2 (x+h), y 3 (x+h), etc. This is because the new value y i (x+h) depends (or might

depend), according to (2.3.8), on the old values of all of the unknown functions, including

those whose new values have already been computed before we begin the computation of

y i (x + h).

If the point still is murky, go back to (2.3.11) and notice how, in the calculation of

y 2 (0.10) we needed to know y 1 (0.05) even though y1 (0.10) had already been computed.

Hence if we had put y 1 (0.10) into an array to replace the old value y 1 (0.05) we would not

have been able to obtain y 2 (0.10).

The conclusion is that we need at least two arrays, say YIN and YOUT, each of length N .

The array YIN holds the unknown functions evaluated at x, and YOUT will hold their values

at x + h. Initially YIN holds the given data at x0 . Then we compute all of the unknowns at

x0 + h, and store them in YOUT as we find them. When all have been done, we print them

if desired, move all entries of YOUT back to YIN, increase x by h and repeat.

32

The Numerical Solution of Differential Equations

The principal building block in this structure would be a subroutine that would advance

the solution exactly one step. The main program would initialize the arrays, call this

subroutine, increase x, move date from the output array YOUT back to the input array YIN,

print, etc. The single-step subroutine is shown below. We will use it later on to help us

get a solution started when we use methods that need information at more than one point

before they can start the integration.

Eulerstep:=proc(xin,yin,h,n) local i,yout;

# This program numerically integrates the system

# y’=f(x,y) one step forward by Euler’s method using step size

# h. Enter with values at xin in yin. Exit with values at xin+h

# in yout. Supply f as a function subprogram.

yout:=[seq(evalf(yin[i]+h*f(xin,yin,i)),i=1..n)];

return(yout);

end;

A few remarks about the program are in order. One structured data type in Maple

is a list of things, in this case, a list of floating point numbers. The seq command (for

“sequence”) creates such a list, in this case a list of length n since i goes from 1 to n in the

seq command. The brackets [ and ] convert the list into a vector. The evalf command

ensures that the results of the computation of the components of yout are floating point

numbers.

Our next remark about the program concerns the function subprogram f, which calculates the right hand sides of the differential equation. This subprogram, of course, must be

supplied by the user. Here is a sample of such a program, namely the one that describes

the system (2.3.9). In that case we have f1 (x, y) = x + y1 + y2 and f2 (x, y) = y1 y2 + 1. This

translates into the following:

f:=proc(x,y,i);

# Calculates the right-hand sides of the system of differential

# equations.

if i=1 then return(x+y[1]+y[2]) else return(y[1]*y[2]+1) fi;

end;

Our last comment about the program to solve systems is that it is perfectly possible to

use it in such a way that we would not have to move the contents of the vector YOUT back

to the vector YIN at each step. In other words, we could save N move operations, where N

is the number of equations. Such savings might be significant in an extended calculation.

To achieve this saving, we write two blocks of programming in the main program. One

block takes the contents of YIN as input, advances the solution one step by Euler’s method

and prints, leaving the new vector in YOUT. Then, without moving anything, another block

of programming takes the contents of YOUT as input, advances the solution one step, leaves

the new vector in YIN, and prints. The two blocks call Euler alternately as the integration

proceeds to the right. The reader might enjoy writing this program, and thinking about

how to generalize the idea to the situation where the new value of y is computed from two

2.3 Systems and equations of higher order

33

previously computed values, rather than from just one (then three blocks of programming

would be needed).

Now we’ve discussed the numerical solution of a single differential equation of first order,

and of a system of simultaneous differential equations of first order, and there remains the

treatment of equations of higher order than the first. Fortunately, this case is very easily

reduced to the varieties that we have already studied.

For example, suppose we want to solve a single equation of the second order, say

y 00 + xy 0 + (x + 1) cos y = 2.

(2.3.13)

The strategy is to transform the single second-order equation into a pair of simultaneous

first order equations that can then be handled as before. To do this, choose two unknown

functions u and v. The function u is to be the unknown function y in (2.3.13), and the

function v is to be the derivative of u. Then u and v satisfy two simultaneous first-order

differential equations:

(

u0 = v

(2.3.14)

v 0 = ?xv ? (x + 1) cos u + 2

and these are exactly of the form (2.3.5) that we have already discussed!

The same trick works on a general differential equation of N th order

y (N ) + G(x, y, y 0 , y 00 , . . . , y (N ?1) ) = 0.

(2.3.15)

We introduce N unknown functions u0 , u1 , . . . , uN ?1 , and let them be the solutions of the

system of N simultaneous first order equations

u00 = u1

u01 = u2

...

u0N ?2

u0N ?1

(2.3.16)

= uN ?1

= ?G(x, u0 , u1 , . . . , uN ?2 ).

The system can now be dealt with as before.

Exercises 2.3

1. Write each of the following as a system of simultaneous first-order initial-value problems in the standard form (2.3.2):

(a) y 00 + x2 y = 0; y(0) = 1; y 0 (0) = 0

(b) u0 + xv = 2; v 0 + euv = 0; u(0) = 0; v(0) = 0

(c) u0 + xv 0 = 0; v 0 + x2 u = 1; u(1) = 1; v(1) = 0

(d) y iv + 3xy 000 + x2 y 00 + 2y 0 + y = 0; y(0) = y 0 (0) = y 00 (0) = y 000 (0) = 1

(e) x000 (t) + t3 x(t) + y(t) = 0; y 00 (t) + x(t)2 = t3 ; x(0) = 2; x0 (0) = 1; x00 (0) = 0;

y(0) = 1; y 0 (0) = 0

34

The Numerical Solution of Differential Equations

2. For each of the parts of problem 1, write the function subprogram that will compute

the right-hand sides, as required by the Eulerstep subroutine.

3. For each of the parts of problem 1, assemble nad run on the computer the Euler

program, together with the relevant function subprogram of problem 2, to print out

the solutions for fifty steps of integration, each of size h = 0.03. Begin with x = x0 ,

the point at which the initial data was given.

4. Reprogram the Eulerstep subroutine, as discussed in the text, to avoid the movement

of YOUT back to YIN.

5. Modify your program as necessary (in Maple, take advantage of the plot command)

to produce graphical output (graph all of the unknown functions on the same axes).

Test your program by running it with Euler as it solves y 00 + y = 0 with y(0) = 0,

y 0 (0) = 1, and h = ?/60 for 150 steps.

6. Write a program that will compute successive values yp , yp+1 , . . . from a difference

equation of order p. Do this by storing the y’s as they are computed in a circular list,

so that it is never necessary to move back the last p computed values before finding

the next one. Write your program so that it will work with vectors, so you can solve

systems of difference equations as well as single ones.

2.4

How to document a program

One of the main themes of our study will be the preparation of programs that not only

work, but also are easily readable and useable by other people. The act of communication

that must take place before a program can be used by persons other than its author is a

difficult one to carry out, and we will return several times to the principles that serve as

guides to the preparation of readable software.

In this section we discuss further the all-important question of program documentation,

already touched upon in section 2.2. Some very nontrivial skills are called for in the creation

of good user-oriented program descriptions. One of these is the ability to enter the head of

another person, the user, and to relate to the program that you have just written through

the user’s eyes.

It’s hard to enter someone else’s head. One of the skills that make one person a better

teacher than another person is of the same kind: the ability to see the subject matter that

is being taught through the eyes of another person, the student. If one can do that, or

even make a good try at it, then obviously one will be able much better to deal with the

questions that are really concerning the audience. Relatively few actually do this to any

great extent not, I think, because it’s an ability that one either has or doesn’t have, but

because few efforts are made to train this skill and to develop it.

We’ll try to make our little contribution here.

2.4 How to document a program

35

(A) What does it do?

The first task should be to describe the precise purpose of the program. Put yourself in

the place of a potential user who is looking at a particular numerical instance of the problem

that needs solving. That user is now thumbing through a book full of program descriptions

in the library of a computer center 10,000 miles away in order to find a program that will

do the problem in question. Your description must answer that question.

Let’s now assume that your program has been written in the form of a subroutine or

procedure, rather than as a main program. Then the list of global, or communicating

variables is plainly in view, in the opening statement of the subroutine.

As we noted in section 2.2, you should state the purpose of your program using the

global variables of the subroutine in the same sentence. For one example of the difference

that makes, see section 2.2. For another, a linear equation solver might be described by

saying

“This program solves a system of simultaneous equations. To use it, put the

right-hand sides into the vector B, put the coefficients into the matrix A and call

the routine. The answer will be returned in X.”

We are, however, urging the reader to do it this way:

“This program solves the equations AX=B, where A is an N-by-N matrix and B

is an N-vector.”

Observe that the second description is shorter, only about half as long, and yet more

informative. We have found out not only what the program does, but how that function

relates to the global variables of the subroutine. This was done by using a judicious sprinkling of symbols in the documentation, along with the words. Don’t use only symbols, or

only words, but weave them together for maximum information.

Notice also that the ambiguous term “right-hand side” that appeared in the first form

has been done away with in the second form. The phrase was ambiguous because exactly

what ends up on the right-hand side and what on the left is an accident of how we happen

to write the equations, and your audience may not do it the same way you do.

(B) How is it done?

This is usually the easy part of program documentation because it is not the purpose of

this documentation to give a course in mathematics or algorithms or anything else. Hence

most of the time a reference to the literature is enough, or perhaps if the method is a

standard one, just give its name. Often though, variations on the standard method have

been chosen, and the user must be informed about those:

“. . . is solved by Gaussian elimination, using complete positioning for size. . . ”

“. . . the input array A is sorted by the Quicksort method (see D.E. Knuth, The

Art of Computer Programming, volume 3). . . ”

36

The Numerical Solution of Differential Equations

“. . . the eigenvalues and vectors are found by the Jacobi method, using Corbat?o’s method of avoiding the search for the largest off-diagonal element (see,

for instance, the description in D.R. Wilson, A First Course in Mathematical

Software).”

“. . . is found by the Simplex method, except that Charnes’ selection rule (see

F.A. Ficken, The Simplex Method. . . ) is not programmed, and so. . . ”

(C) Describe the global variables

Now it gets hard again. The global variables are the ones through which the subroutine

communicates with the user. Generally speaking, the user doesn’t care about variables

that are entirely local to your subroutine, but is vitally concerned with the communicating

variables.

First the user has to know exactly ho each of the global variables is related to the

problem that is being solved. This calls for a brief verbal description of the variable, and

what it has to do with the functioning of the program.

“A[i] is the ith element of the input list that is to be sorted, i=1..N”

“WHY is set by the subroutine to TRUE unles the return is because of overflow,

and then it will be set to FALSE.”

“B[i,j] is the coefficient of X[j] in the ith one of the input equations BX=C.”

“option is set by the calling program on input. Set it to 0 if the output is to

be rounded to the nearest integer, else set it to m if the output is to be rounded

to m decimal places (m ? 12).”

It is extremely important that each and every global variable of the subroutine should

get such a description. Just march through the parentheses in the subroutine or procedure

heading, and describe each variable in turn.

Next, the user will need more information about each of the global variables than just its

description as above. Also required is the “type” of the variable. Some computer languages

force each program to declare the types of their variables right in the opening statement.

Others declare types by observing various default rules with exceptions stated. In any case,

a little redundancy never hurts, and the program documentation should declare the type of

each and every global variable.

It’s easy to declare along with the types of the variables, their dimensions if they are

array variables. For instance we may have a

solver:=proc(A,X,n,ndim,b);

2.4 How to document a program

37

in which the communicating variables have the following types:

A

X

n

ndim

b

ndim-by-ndim array of floating point numbers

vector of floating point numbers of length n

integer

integer

vector of floating point numbers of length n

The best way to announce all of these types and dimensions of global variables to the

user is simply to list them, as above, in a table.

Now surely we’ve finished taking the pulse, blood pressure, etc. of the global variables,

haven’t we? Well, no, we haven’t. There’s still more vital data that a user will need to

know about these variables. There isn’t any standard name like “type” to apply to this

information, so we’ll call it the “role” of the variable.

First, for some of the global variables of the subroutine, it may be true that their values

at the time the subroutine is called are quite irrelevant to the operation of the subroutine.

This would be the case for the output variables, and in certain other situations. For some

other variables, the values at input time would be crucial. The user needs to know which are

which. Just for one example, if the value at input time is irrelevant, then the user can feel

free to use the same storage for other temporary purposes between calls to the subroutine.

Second, it may happen that certain variables are returned by the subroutine with their

values unchanged. This is particularly true for “implicitly passed” global variables, i.e.,

variables whose values are used by the subroutine but which do not appear explicitly in the

argument list. In such cases, the user may be delighted to hear the good news. In other

cases, the action of a subroutine may change an input variable, so if the user needs to use

those quantities again it will be necessary to save them somewhere else before calling the

subroutine. In either case, the user needs to be informed.

Third, it may be that the computation of the value of a certain variable is one of the

main purposes of the subroutine. Such variables are the outputs of the program, and the

user needs to know which these are (whether they are explicit in heading or the return

statement, or are “implicit”).

Although some high-level computer languages require type declarations immediately

in the opening instruction of a subroutine, none require the descriptions of the roles of

the variables (well, Pascal requires the VAR declaration, and Maple separates the input

variables from the output ones, but both languages allow implicit passing and changing

of global variables). These are, however, important for the user to know, so let’s invent

a shorthand for describing them in the documentation of the programs that occur in this

book.

First, if the value at input time is important, let’s say that the role of the variable is I,

otherwise it is I’.

Second, if the value of the variable is changed by the action of the subroutine, we’ll say

that its role is C, else C’.

Finally, if the computation of this variable is one of the main purposes of the subroutine,

38

The Numerical Solution of Differential Equations

it’s role is O (as in output), else O’.

In the description of each communicating variable, all three of these should be specified,.

Thus, a variable X might have role IC’O’, or a variable why might be of role I’CO, etc.

To sum up, the essential features of program documentation are a description of that

the program does, phrased in terms of the global variables, a statement of how it gets the

job done, and a list of all of the global variables, showing for each one its name, type,

dimension (or structure) if any, its role, and a brief verbal description.

Refer back to the short program in section 2.2, that searches for the largest element in

a row of a matrix. Here is the table of information about its global variables:

Name

A

i

jwin

Type

floating point matrix

integer

integer

Role

IC’O’

IC’O’

I’CO

Description

The input matrix

Which row to search

Column containing largest element

Exercises 2.4

Write programs that perform each of the jobs stated below. In each case, after testing

the program, document it with comments. Give a complete table of information about the

global variables in each case.

(a) Find and print all of the prime numbers between M and N .

(b) Find the elements of largest and smallest absolute values in a given linear array

(vector), and their positions in the array.

(c) Sort the elements of a given linear array into ascending order of size.

(d) Deal out four bridge hands (13 cards each from a standard 52-card deck – This one

is not so easy!).

(e) Solve a quadratic equation (any quadratic equation!).

2.5

The midpoint and trapezoidal rules

Euler’s formula is doubtless the simplest numerical integration procedure for differential

equations, but the accuracy that can be obtained with it is insufficient for most applications.

In this section and those that follow, we want to introduce a while family of methods for

the solution of differential equations, called the linear multistep methods, in which the user

can choose the degree of precision that will suffice for the job, and then select a member of

the family that will achieve it.

Before describing the family in all of its generality, we will produce two more of its

members, which illustrate the different sorts of creatures that inhabit the family in question.

2.5 The midpoint and trapezoidal rules

39

Recall that we derived Euler’s method by chopping off the Taylor series expansion of the

solution after the linear term. To get a more accurate method we could, of course, keep the

quadratic term, too. However, that term involves a second derivative, and we want to avoid

the calculation of higher derivatives because our differential equations will always be written

as first-order systems, so that only the first derivative will be conveniently computable.

We can have greater accuracy without having to calculate higher derivatives if we’re

willing to allow our numerical integration procedure to involve values of the unknown function and its derivative at more than one point. In other words, in Euler’s method, the next

value of the unknown function, at x + h, is gotten from the values of y and y 0 at just one

backwards point x. In the more accurate formulas that we will discuss next, the new value

of y depends in y and y 0 at more than one point, for instance, at x and x ? h, or at several

points.

As a primitive example of this kind, we will now discuss the midpoint rule. We begin

once again with the Taylor expansion of the unknown function y(x) about the point xn :

y(xn + h) = y(xn ) + hy 0 (xn ) + h2

y 00 (xn )

y 000 (xn )

+ h3

+ ···.

2

6

(2.5.1)

Now we rewrite equation (2.5.1) with h replaced by ?h to get

y(xn ? h) = y(xn ) ? hy 0 (xn ) + h2

y 00 (xn )

y 000 (xn )

? h3

+ ···

2

6

(2.5.2)

and then subtract these equations, obtaining

y(xn + h) ? y(xn ? h) = 2hy 0 (xn ) + 2h3

y 000 (xn )

+ ···.

6

(2.5.3)

Now, just as we did in the derivation of Euler’s method, we will truncate the right side

of (2.5.3) after the first term, ignoring the terms that involve h3 , h5 , etc. Further, let’s

use yn to denote the computed approximate value of y(xn ) (and yn+1 for the approximate

y(xn+1 ), etc.). Then we have

yn+1 ? yn?1 = 2hyn0 .

(2.5.4)

If, as usual, we are solving the differential equation y 0 = f (x, y), then finally (2.5.4) takes

the form

yn+1 = yn?1 + 2hf (xn , yn )

(2.5.5)

and this is the midpoint rule. The name arises from the fact that the first derivative yn0 is

being approximated by the slope of the chord that joins the two points (xn?1 , yn?1 ) and

(xn+1 , yn+1 ), instead of the chord joining (xn , yn ) and (xn+1 , yn+1 ) as in Euler’s method.

At first sight it seems that (2.5.5) can be used just like Euler’s method, because it is a

recurrence formula in which we compute the next value yn+1 from the two previous values

yn and yn?1 . Indeed the rules are quite similar, except for the fact that we can’t get started

with the midpoint rule until we know two consecutive values y0 , y1 of the unknown function

at two consecutive points x0 , x1 . Normally a differential equation is given together with

just one value of the unknown function, so if we are to use the midpoint rule we’ll need to

manufacture one more value of y(x) by some other means.

40

The Numerical Solution of Differential Equations

This kind of situation will come up again and again as we look at more accurate methods,

because to obtain greater precision without computing higher derivatives we will get the

next approximate value of y from a recurrence formula that may involve not just one or two,

but several of its predecessors. To get such a formula started we will have to find several

starting values in addition to the one that is given in the statement of the initial-value

problem.

To get back to the midpoint rule, we can get it started most easily by calculating y1 ,

the approximation to y(x0 + h), from Euler’s method, and then switching to the midpoint

rule to carry out the rest of the calculation.

Let’s do this, for example with the same differential equation (2.1.7) that we used to

illustrate Euler’s rule, so we can compare the two methods. The problem consists of the

equation y 0 = 0.5y and the initial value y(0) = 1. We’ll use the same step size h = 0.05 as

before.

Now to start the midpoint rule we need two consecutive values of y, in this case at x = 0

and x = 0.05. At 0.05 we use the value that Euler’s method gives us, namely y1 = 1.025

(see Table 1). It’s easy to continue the calculation now from (2.5.5).

For instance

y2 = y0 + 2h(0.5y1 )

= 1 + 0.1(0.5 ? 1.025)

(2.5.6)

= 1.05125

and

y3 = y1 + 2h(0.5y2 )

= 1.025 + 0.1(0.5 ? 1.05125)

(2.5.7)

= 1.0775625 .

In the table below we show for each x the value computed from the midpoint rule, from

Euler’s method, and from the exact solution y(x) = ex/2 . The superior accuracy of the

midpoint rule is apparent.

x

0.00

0.05

0.10

0.15

0.20

0.25

..

.

Midpoint(x)

1.00000

1.02500

1.05125

1.07756

1.10513

1.13282

..

.

Euler(x)

1.00000

1.02500

1.05063

1.07689

1.10381

1.13141

..

.

Exact(x)

1.00000

1.02532

1.05127

1.07788

1.10517

1.13315

..

.

1.00

2.00

3.00

..

.

1.64847

2.71763

4.48032

..

.

1.63862

2.68506

4.39979

..

.

1.64872

2.71828

4.48169

..

.

5.00

10.00

12.17743

148.31274

11.81372

139.56389

12.18249

148.41316

2.5 The midpoint and trapezoidal rules

41

y = f (x)

a

b

Figure 2.1: The trapezoidal rule

table 2

Next, we introduce a third method of numerical integration, the trapezoidal rule. The

best way to obtain it is to convert the differential equation that we’re trying to solve into

an integral equation, and then use the trapezoidal approximation for the integral.

We begin with the differential equation y 0 = f (x, y(x)), and we integrate both sides

from x to x + h, getting

Z

x+h

y(x + h) = y(x) +

f (t, y(t)) dt.

(2.5.8)

x

Now if we approximate the right-hand side in any way by a weighted sum of values of

the integrand at various points we will have found an approximate method for solving our

differential equation.

The trapezoidal rule states that for an approximate value of an integral

Z

b

f (t) dt

(2.5.9)

a

we can use, instead of the area under the curve between x = a and x = b, the area of the

trapezoid whose sides are the x axis, the lines x = a and x = b, and the line through the

points (a, f (a)) and (b, f (b)), as shown in Figure 2.1. That area is 12 (f (a) + f (b))(b ? a).

If we apply the trapezoidal rule to the integral that appears in (2.5.8), we obtain

y(xn + h) ? y(xn ) +

h

(f (xn , y(xn )) + f (xn + h, y(xn + h)))

2

(2.5.10)

in which we have used the “?” sign rather than the “=” because the right hand side is not

exactly equal to the integral that really belongs there, but is only approximately so.

If we use our usual abbreviation yn for the computed approximate value of y(xn ), then

(2.5.10) becomes

h

yn+1 = yn + (f (xn , yn ) + f (xn+1 , yn+1 )).

(2.5.11)

2

42

The Numerical Solution of Differential Equations

This is the trapezoidal rule in the form that is useful for differential equations.

At first sight, (2.5.11) looks like a recurrence formula from which the next approximate

value, yn+1 , of the unknown function, can immediately be computed from the previous

value, yn . However, this is not the case.

Upon closer examination one observes that the next value yn+1 appears not only on the

left-hand side, but also on the right (it’s hiding in the second f on the right side).

In order to find the value yn+1 it appears that we need to carry out an iterative process.

First we would guess yn+1 (guessing yn+1 to be equal to yn wouldn’t be all that bad, but

we can do better). If we use this guess value on the right side of (2.5.11) then we would

be able to calculate the entire right-hand side, and then we could use that value as a new

“improved” value of yn+1 .

Now if the new value agrees with the old sufficiently well the iteration would halt, and

we would have found the desired value of yn+1 . Otherwise we would use the improved value

on the right side just as we previously used the first guess. Then we would have a “more

improved” guess, etc.

Fortunately, in actual use, it turns out that one does not actually have to iterate to

convergence. If a good enough guess is available for the unknown value, then just one

refinement by a single application of the trapezoidal formula is sufficient. This is not the

case if a high quality guess is unavailable. We will discuss this point in more detail in

section 2.9. The pair of formulas, one of which supplies a very good guess to the next

value of y, and the other of which refines it to a better guess, is called a predictor-corrector

pair, and such pairs form the basis of many of the highly accurate schemes that are used in

practice.

As a numerical example, take the differential equation

y 0 = 2xey + 1

(2.5.12)

with the initial value y(0) = 1. If we use h = 0.05, then our first task is to calculate y1 , the

approximate value of y(0.05). The trapezoidal rule asserts that

y1 = 1 + 0.025(2 + 0.1ey1 )

(2.5.13)

and sure enough, the unknown number y1 appears on both sides.

Let’s guess y1 = 1. Since this is not a very inspired guess, we will have to iterate the

trapezoidal rule to convergence. Hence, we use this guess on the right side of (2.5.13),

compute the right side, and obtain y1 = 1.056796. If we use this new guess the same way,

the result is y1 = 1.057193. Then we get 1.057196, and since this is in sufficiently close

agreement with the previous result, we declare that the iteration has converged. Then we

take y1 = 1.057196 for the computed value of the unknown function at x = 0.05, and we go

next to x = 0.1 to repeat the same sort of thing to get y2 , the computed approximation to

y(0.1).

In Table 3 we show the results of using the trapezoidal rule (where we have iterated

until two successive guesses are within 10?4 ) on our test equation y 0 = 0.5y, y(0) = 1 as

the column Trap(x). For comparison, we show Midpoint(x) and Exact(x).

2.6 Comparison of the methods

43

x

0.00

0.05

0.10

0.15

0.20

0.25

..

.

Trap(x)

1.00000

1.02532

1.05127

1.07789

1.10518

1.13316

..

.

Midpoint(x)

1.00000

1.02500

1.05125

1.07756

1.10513

1.13282

..

.

Exact(x)

1.00000

1.02532

1.05127

1.07788

1.10517

1.13315

..

.

1.00

2.00

3.00

..

.

1.64876

2.71842

4.48203

..

.

1.64847

2.71763

4.48032

..

.

1.64872

2.71828

4.48169

..

.

5.00

10.00

12.18402

148.45089

12.17743

148.31274

12.18249

148.41316

table 3

2.6

Comparison of the methods

We are now in possession of three methods for the numerical solution of differential equations. They are Euler’s method

yn+1 = yn + hyn0 ,

(2.6.1)

the trapezoidal rule

yn+1 = yn +

h 0

0

),

(y + yn+1

2 n

(2.6.2)

and the midpoint rule

yn+1 = yn?1 + 2hyn0 .

(2.6.3)

In order to compare the performance of the three techniques it will be helpful to have

a standard differential equation on which to test them. The most natural candidate for

such an equation is y 0 = Ay, where A is constant. The reasons for this choice are first

that the equation is easy to solve exactly, second that the difference approximations are

also relatively easy to solve exactly, so comparison is readily done, third that by varying the

sign of A we can study behavior of either growing or shrinking (stable or unstable) solutions,

and finally that many problems in nature have solutions that are indeed exponential, at

least over the short term, so this is an important class of differential equations.

We will, however, write the test equation in a slightly different form for expository

reasons, namely as

y

y0 = ? ;

y(0) = 1 ;

L>0

(2.6.4)

L

where L is a constant. The most interesting and revealing case is where the true solution

is a decaying exponential, so we will assume that L > 0. Further, we will assume that

y(0) = 1 is the given initial value.

44

The Numerical Solution of Differential Equations

The exact solution is of course

Exact(x) = e?x/L .

(2.6.5)

Notice that if x increases by L, the solution changes by a factor of e. Hence L, called the

relaxation length of the problem, can be conveniently visualized as the distance over which

the solution falls by a factor of e.

Now we would like to know how well each of the methods (2.6.1)–(2.6.3) handles the

problem (2.6.4).

Suppose first that we ask Euler’s method to solve the problem. If we substitute y 0 =

f (x, y) = ?y/L into (2.6.1), we get

yn+1

yn

= yn + h ? ?

L

h

= 1?

yn .

L

(2.6.6)

Before we solve this recurrence, let’s comment on the ratio h/L that appears in it. Now

L is the distance over which the solution changes by a factor of e, and h is the step size

that we are going to use in the numerical integration. Instinctively, one feels that if the

solution is changing rapidly in a certain region, then h will have to be kept small there if

good accuracy is to be retained, while if the solution changes only slowly, then h can be

larger without sacrificing too much accuracy. The ratio h/L measures the step size of the

integration in relation to the distance over which the solution changes appreciably. Hence,

h/L is exactly the thing that one feels should be kept small for a successful numerical

solution.

Since h/L occurs frequently below, we will denote it with the symbol ? .

Now the solution of the recurrence equation (2.6.6), with the starting value y0 = 1, is

obviously

yn = (1 ? ? )n

n = 0, 1, 2, . . . .

(2.6.7)

Next we study the trapezoidal approximation to the same equation (2.6.4). We substitute y 0 = f (x, y) = ?y/L in (2.6.2) and get

yn+1

h

yn yn+1

= yn +

? ?

.

2

L

L

(2.6.8)

The unknown yn+1 appears, as usual with this method, on both sides of the equation.

However, for the particularly simple equation that we are now studying, there is no difficulty

in solving (2.6.8) for yn+1 (without any need for an iterative process) and obtaining

yn+1 =

1?

1+

?

2

?

2

yn .

(2.6.9)

Together with the initial value y0 = 1, this implies that

yn =

1?

1+

?

2

?

2

!n

n = 0, 1, 2, . . . .

(2.6.10)

2.6 Comparison of the methods

45

Before we deal with the midpoint rule, let’s pause to examine the two methods whose

solutions we have just found. Note that for a given value of h, all three of (a) the exact solution, (b) Euler’s solution and (c) the trapezoidal solution are of the form yn = (constant)n ,

in which the three values of “constant” are

(a) e??

(b) 1 ? ?

(c)

1? ?2

1+ ?2

(2.6.11)

.

It follows that to compare the two approximate methods with the “truth,” all we have

to do is see how close the constants (b) and (c) above are to the true constant (a). If we

remember that ? is being thought of as small compared to 1, then we have the power series

expansion of e??

?2 ?3

e?? = 1 ? ? +

?

+ ···

(2.6.12)

2

6

to compare with 1 ? ? and with the power series expansion of

1?

1+

?

2

?

2

=1?? +

?2 ?3

?

+ ···.

2

4

(2.6.13)

The comparison is now clear. Both the Euler and the trapezoidal methods yield approximate solutions of the form (constant)n , where “constant” is near e?? . The trapezoidal

rule does a better job of being near e?? because its constant agrees with the power series

expansion of e?? through the quadratic term, whereas that of the Euler method agrees only

up to the linear term.

Finally we study the nature of the approximation that is provided by the midpoint

rule. We will find that a new and important phenomenon rears its head in this case. The

analysis begins just as it did in the previous two cases: We substitute the right-hand side

f (x, y) = ?y/L for y 0 in (2.6.3) to get

yn+1 = yn?1 + 2h ? ?

yn

.

L

(2.6.14)

One important feature is already apparent. Instead of facing a first-order difference

equation as we did in (2.6.6) for Euler’s method and in (2.6.9) for the trapezoidal rule, we

have now to contend with a second-order difference equation.

Since the equation is linear with constant coefficients, we know to try a solution of the

form yn = r n . This leads to the quadratic equation

r 2 + 2? r ? 1 = 0.

(2.6.15)

Evidently the discriminant of this equation is positive, so its roots are distinct. If we denote

these two roots by r+ (? ) and r? (? ), then the general solution of the difference equation

(2.6.14) is

yn = c (r+ (? ))n + d (r? (? ))n ,

(2.6.16)

46

The Numerical Solution of Differential Equations

where c and d are constants whose values are determined by the initial data y0 and y1 .

The Euler and trapezoidal approximations were each of the form (constant)n . This one

is a sum of two terms of that kind. We will see that r+ (? ) is a very good approximation

to e?? . The other term, (r? (? ))n is, so to speak, the price that we pay for getting such a

good approximation in r+ (? ). We hope that the other term will stay small relative to the

first term, so as not to disturb the closeness of the approximation. We will see, however,

that it need not be so obliging, and that in fact it might do whatever it can to spoil things.

The two roots of the quadratic equation are

?

r+ (? ) = ?? + ?1 + ? 2

r? (? ) = ?? ? 1 + ? 2 .

(2.6.17)

When ? = 0 the first of these is +1, so when ? is small r+ (? ) is near +1, and it is the root

that is trying to approximate the exact constant e?? as well as possible. In fact it does

pretty well, because the power series expansion of r+ (? ) is

r+ (? ) = 1 ? ? +

?2 ?4

?

+ ···

2

8

(2.6.18)

so it agrees with e?? through the quadratic terms.

What about r? (? )? Its Taylor series is

r? (? ) = ?1 ? ? ?

?2

+ ···.

2

(2.6.19)

The bad news is now before us: When ? is a small positive number, the root r? (? ) is

larger than 1 in absolute value. This means that the stability criterion of Theorem 1.6.1 is

violated, so we say that the midpoint rule is unstable.

In practical terms, we observe that r+ (? ) is close to e?? , so the first term on the right

of (2.6.16) is very close to the exact solution. The second term of (2.6.16), the so-called

parasitic solution, is small compared to the first term when n is small, because the constant

d will be small compared with c. However, as we move to the right, n increases, and the

second term will eventually dominate the first, because the first term is shrinking to zero

as n increases, because that’s what the exact solution does, while the second term increases

steadily in size. In fact, since r? (? ) is negative and larger than 1 in absolute value, the

second term will alternate in sign as n increases, and grow without bound in magnitude.

In Table 4 below we show the result of integrating the problem y 0 = ?y, y(0) = 1 with

each of the three methods that we have discussed, using a step size of h = 0.05. To get the

midpoint method started, we used the exact value of y(0.05) (i.e., we cheated), and in the

trapezoidal rule we iterated to convergence with = 10?4 . The instability of the midpoint

rule is quite apparent.

2.6 Comparison of the methods

x

0.0

1.0

2.0

3.0

4.0

5.0

10.0

14.55

15.8

Euler(x)

1.00000

0.35849

0.12851

0.04607

0.01652

0.00592

0.00004

3.3 ? 10?7

9.1 ? 10?8

47

Trap(x)

1.00000

0.36780

0.13527

0.04975

0.01830

0.00673

0.00005

4.8 ? 10?7

1.4 ? 10?7

Midpoint(x)

1.00000

0.36806

0.13552

0.05005

0.01888

0.00822

0.21688

-20.48

71.45

Exact(x)

1.00000

0.36788

0.13534

0.04979

0.01832

0.00674

0.00005

4.8 ? 10?7

1.4 ? 10?7

table 4

In addition to the above discussion of accuracy, we summarize here there additional

properties of integration methods as they relate to the examples that we have already

studied.

First, a numerical integration method might be iterative or noniterative. A method is

noniterative if it expresses the next value of the unknown function quite explicitly in terms

of values of the function and its derivatives at preceding points. In an iterative method,

at each step of the solution process the next value of the unknown is defined implicitly by

an equation, which must be solved to obtain the next value of the unknown function. In

practice, we may either solve this equation completely by an iteration or

Dennis Deturck and Herbert S. Wilf

Department of Mathematics

University of Pennsylvania

Philadelphia, PA 19104-6395

Copyright 2002, Dennis Deturck and Herbert Wilf

April 30, 2002

2

Contents

1 Differential and Difference Equations

5

1.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2

Linear equations with constant coefficients . . . . . . . . . . . . . . . . . . .

8

1.3

Difference equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

1.4

Computing with difference equations . . . . . . . . . . . . . . . . . . . . . .

14

1.5

Stability theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

1.6

Stability theory of difference equations . . . . . . . . . . . . . . . . . . . . .

19

2 The Numerical Solution of Differential Equations

23

2.1

Euler’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

2.2

Software notes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

2.3

Systems and equations of higher order . . . . . . . . . . . . . . . . . . . . .

29

2.4

How to document a program . . . . . . . . . . . . . . . . . . . . . . . . . .

34

2.5

The midpoint and trapezoidal rules . . . . . . . . . . . . . . . . . . . . . . .

38

2.6

Comparison of the methods . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

2.7

Predictor-corrector methods . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

2.8

Truncation error and step size . . . . . . . . . . . . . . . . . . . . . . . . . .

50

2.9

Controlling the step size . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

2.10 Case study: Rocket to the moon . . . . . . . . . . . . . . . . . . . . . . . .

60

2.11 Maple programs for the trapezoidal rule . . . . . . . . . . . . . . . . . . . .

65

2.11.1 Example: Computing the cosine function . . . . . . . . . . . . . . .

67

2.11.2 Example: The moon rocket in one dimension . . . . . . . . . . . . .

68

2.12 The big leagues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

2.13 Lagrange and Adams formulas . . . . . . . . . . . . . . . . . . . . . . . . .

74

4

CONTENTS

3 Numerical linear algebra

81

3.1

Vector spaces and linear mappings . . . . . . . . . . . . . . . . . . . . . . .

81

3.2

Linear systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

86

3.3

Building blocks for the linear equation solver . . . . . . . . . . . . . . . . .

92

3.4

How big is zero? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

3.5

Operation count . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

102

3.6

To unscramble the eggs . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

105

3.7

Eigenvalues and eigenvectors of matrices . . . . . . . . . . . . . . . . . . . .

108

3.8

The orthogonal matrices of Jacobi . . . . . . . . . . . . . . . . . . . . . . .

112

3.9

Convergence of the Jacobi method . . . . . . . . . . . . . . . . . . . . . . .

115

3.10 Corbat?o’s idea and the implementation of the Jacobi algorithm . . . . . . .

118

3.11 Getting it together . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

122

3.12 Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

124

Chapter 1

Differential and Difference

Equations

1.1

Introduction

In this chapter we are going to study differential equations, with particular emphasis on how

to solve them with computers. We assume that the reader has previously met differential

equations, so we’re going to review the most basic facts about them rather quickly.

A differential equation is an equation in an unknown function, say y(x), where the

equation contains various derivatives of y and various known functions of x. The problem

is to “find” the unknown function. The order of a differential equation is the order of the

highest derivative that appears in it.

Here’s an easy equation of first order:

y 0 (x) = 0.

(1.1.1)

The unknown function is y(x) = constant, so we have solved the given equation (1.1.1).

The next one is a little harder:

y 0 (x) = 2y(x).

(1.1.2)

A solution will, now doubt, arrive after a bit of thought, namely y(x) = e2x . But, if y(x)

is a solution of (1.1.2), then so is 10y(x), or 49.6y(x), or in fact cy(x) for any constant c.

Hence y = ce2x is a solution of (1.1.2). Are there any other solutions? No there aren’t,

because if y is any function that satisfies (1.1.2) then

(ye?2x )0 = e?2x (y 0 ? 2y) = 0,

(1.1.3)

and so ye?2x must be a constant, C.

In general, we can expect that if a differential equation is of the first order, then the

most general solution will involve one arbitrary constant C. This is not always the case,

6

Differential and Difference Equations

since we can write down differential equations that have no solutions at all. We would have,

for instance, a fairly hard time (why?) finding a real function y(x) for which

(y 0 )2 = ?y 2 ? 2.

(1.1.4)

There are certain special kinds of differential equations that can always be solved, and

it’s often important to be able to recognize them. Among there are the “first-order linear”

equations

y 0 (x) + a(x)y(x) = 0,

(1.1.5)

where a(x) is a given function of x.

Before we describe the solution of these equations, let’s discuss the word linear. To say

that an equation is linear is to say that if we have any two solutions y1 (x) and y2 (x) of the

equation, then c1 y1 (x) + c2 y2 (x) is also a solution of the equation, where c1 and c2 are any

two constants (in other words, the set of solutions forms a vector space).

Equation (1.1.1) is linear, in fact, y1 (x) = 7 and y2 (x) = 23 are both solutions, and so

is 7c1 + 23c2 . Less trivially, the equation

y 00 (x) + y(x) = 0

(1.1.6)

is linear. The linearity of (1.1.6) can be checked right from the equation itself, without

knowing what the solutions are (do it!). For an example, though, we might note that

y = sin x is a solution of (1.1.6), that y = cos x is another solution of (1.1.6), and finally,

by linearity, that the function y = c1 sin x + c2 cos x is a solution, whatever the constants c1

and c2 . Now let’s consider an instance of the first order linear equation (1.1.5):

y 0 (x) + xy(x) = 0.

(1.1.7)

So we’re looking for a function whose derivative is ?x times the function. Evidently y =

2

2

e?x /2 will do, and the general solution is y(x) = ce?x /2 .

If instead of (1.1.7) we had

y 0 (x) + x2 y(x) = 0,

then we would have found the general solution ce?x

As a last example, take

3 /3

.

y 0 (x) ? (cos x) y(x) = 0.

(1.1.8)

The right medicine is now y(x) = esin x . In the next paragraph we’ll give the general rule

of which the above are three examples. The reader might like to put down the book at this

point and try to formulate the rule for solving (1.1.5) before going on to read about it.

Ready? What we need is to choose some antiderivative A(x) of a(x), and then the

solution is y(x) = ce?A(x) .

Since that was so easy, next let’s put a more interesting right hand side into (1.1.5), by

considering the equation

y 0 (x) + a(x)y(x) = b(x)

(1.1.9)

1.1 Introduction

7

where now b(x) is also a given function of x (Is (1.1.9) a linear equation? Are you sure?).

To solve (1.1.9), once again choose some antiderivative A(x) of a(x), and then note that

we can rewrite (1.1.9) in the equivalent form

e?A(x)

d A(x)

e

y(x) = b(x).

dx

Now if we multiply through by eA(x) we see that

d A(x)

e

y(x) = b(x)eA(x)

dx

(1.1.10)

so , if we integrate both sides,

Z

eA(x) y(x) =

x

b(t)eA(t) dt + const. ,

(1.1.11)

where on the right side, we mean any antiderivative of the function under the integral sign.

Consequently

Z x

y(x) = e?A(x)

b(t)eA(t) dt + const. .

(1.1.12)

As an example, consider the equation

y0 +

y

= x + 1.

x

(1.1.13)

We find that A(x) = log x, then from (1.1.12) we get

Z

x

1

y(x) =

(t + 1)t dt + C

x

x2 x C

=

+ + .

3

2

x

(1.1.14)

We may be doing a disservice to the reader by beginning with this discussion of certain

types of differential equations that can be solved analytically, because it would be erroneous

to assume that most, or even many, such equations can be dealt with by these techniques.

Indeed, the reason for the importance of the numerical methods that are the main subject

of this chapter is precisely that most equations that arise in “real” problems are quite

intractable by analytical means, so the computer is the only hope.

Despite the above disclaimer, in the next section we will study yet another important

family of differential equations that can be handled analytically, namely linear equations

with constant coefficients.

Exercises 1.1

1. Find the general solution of each of the following equations:

(a) y 0 = 2 cos x

8

Differential and Difference Equations

2

y=0

x

(c) y 0 + xy = 3

1

(d) y 0 + y = x + 5

x

(e) 2yy 0 = x + 1

(b) y 0 +

2. Show that the equation (1.1.4) has no real solutions.

3. Go to your computer or terminal and familiarize yourself with the equipment, the

operating system, and the specific software you will be using. Then write a program

that will calculate and print the sum of the squares of the integers 1, 2, . . . , 100. Run

this program.

4. For each part of problem 1, find the solution for which y(1) = 1.

1.2

Linear equations with constant coefficients

One particularly pleasant, and important, type of linear differential equation is the variety

with constant coefficients, such as

y 00 + 3y 0 + 2y = 0 .

(1.2.1)

It turns out that what we have to do to solve these equations is to try a solution of a certain

form, and we will then find that all of the solutions indeed are of that form.

Let’s see if the function y(x) = e?x is a solution of (1.2.1). If we substitute in (1.2.1),

and then cancel the common factor e?x , we are left with the quadratic equation

?2 + 3? + 2 = 0

whose solutions are ? = ?2 and ? = ?1. Hence for those two values of ? our trial function

y(x) = e?x is indeed a solution of (1.2.1). In other words, e?2x is a solution, e?x is a

solution, and since the equation is linear,

y(x) = c1 e?2x + c2 e?x

(1.2.2)

is also a solution, where c1 and c2 are arbitrary constants. Finally, (1.2.2) must be the most

general solution since it has the “right” number of arbitrary constants, namely two.

Trying a solution in the form of an exponential is always the correct first step in solving

linear equations with constant coefficients. Various complications can develop, however, as

illustrated by the equation

y 00 + 4y 0 + 4y = 0 .

(1.2.3)

Again, let’s see if there is a solution of the form y = e?x . This time, substitution into

(1.2.3) and cancellation of the factor e?x leads to the quadratic equation

?2 + 4? + 4 = 0,

(1.2.4)

1.2 Linear equations with constant coefficients

9

whose two roots are identical, both being ?2. Hence e?2x is a solution, and of course so is

c1 e?2x , but we don’t yet have the general solution because there is, so far, only one arbitrary

constant. The difficulty, of course, is caused by the fact that the roots of (1.2.4) are not

distinct.

In this case, it turns out that xe?2x is another solution of the differential equation (1.2.3)

(verify this), so the general solution is (c1 + c2 x)e?2x .

Suppose that we begin with an equation of third order, and that all three roots turn

out to be the same. For instance, to solve the equation

y 000 + 3y 00 + 3y 0 + y = 0

(1.2.5)

we would try y = e?x , and we would then be facing the cubic equation

?3 + 3?2 + 3? + 1 = 0 ,

(1.2.6)

whose “three” roots are all equal to ?1. Now, not only is e?x a solution, but so are xe?x

and x2 e?x .

To see why this procedure works in general, suppose we have a linear differential equation

with constant coeficcients, say

y (n) + a1 y (n?1) + a2 y (n?2) + · · · + an y = 0

(1.2.7)

If we try to find a solution of the usual exponential form y = e?x , then after substitution into

(1.2.7) and cancellation of the common factor e?x , we would find the polynomial equation

?n + a1 ?n?1 + a2 ?n?2 + · · · + an = 0 .

(1.2.8)

The polynomial on the left side is called the characteristic polynomial of the given

differential equation. Suppose now that a certain number ? = ?? is a root of (1.2.8) of

multiplicity p. To say that ?? is a root of multiplicity p of the equation is to say that

(? ? ?? )p is a factor of the characteristic polynomial. Now look at the left side of the given

differential equation (1.2.7). We can write it in the form

(D n + a1 D n?1 + a2 D n?2 + · · · + an )y = 0 ,

(1.2.9)

in which D is the differential operator d/dx. In the parentheses in (1.2.9) we see the

polynomial ?(D), where ? is exactly the characteristic polynomial in (1.2.8).

Since ?(?) has the factor (? ? ?? )p , it follows that ?(D) has the factor (D ? ?? )p , so

the left side of (1.2.9) can be written in the form

g(D)(D ? ?? )p y = 0 ,

(1.2.10)

?

where g is a polynomial of degree n ? p. Now it’s quite easy to see that y = xk e? x satisfies

(1.2.10) (and therefore (1.2.7) also) for each k = 0, 1, . . . , p ? 1. Indeed, if we substitute this

function y into (1.2.10), we see that it is enough to show that

?

(D ? ?? )p (xk e? x ) = 0

k = 0, 1, . . . , p ? 1 .

(1.2.11)

10

Differential and Difference Equations

?

?

However, (D ? ?? )(xk e?? x ) = kxk?1 e? x , and if we apply (D ? ?? ) again,

?

?

(D ? ?? )2 (xk e?? x ) = k(k ? 1)xk?2 e? x ,

?

etc. Now since k < p it is clear that (D ? ?? )p (xk e?? x ) = 0, as claimed.

To summarize, then, if we encounter a root ?? of the characteristic equation, of multiplicity p, then corresponding to ?? we can find exactly p linearly independent solutions of

the differential equation, namely

?

?

?

?

e? x , xe? x , x2 e? x , . . . , xp?1 e? x .

(1.2.12)

Another way to state it is to say that the portion of the general solution of the given

differential equation that corresponds to a root ?? of the characteristic polynomial equation

?

is Q(x)e? x , where Q(x) is an arbitrary polynomial whose degree is one less than the

multiplicity of the root ?? .

One last mild complication may arise from roots of the characteristic equation that are

not real numbers. These don’t really require any special attention, but they do present a few

options. For instance, to solve y 00 + 4y = 0, we find the characteristic equation ?2 + 4 = 0,

and the complex roots ±2i. Hence the general solution is obtained by the usual rule as

y(x) = c1 e2ix + c2 e?2ix .

(1.2.13)

This is a perfectly acceptable form of the solution, but we could make it look a bit prettier

by using deMoivre’s theorem, which says that

e2ix = cos 2x + i sin 2x

e?2ix = cos 2x ? i sin 2x.

(1.2.14)

Then our general solution would look like

y(x) = (c1 + c2 ) cos 2x + (ic1 ? ic2 ) sin 2x.

(1.2.15)

But c1 and c2 are just arbitrary constants, hence so are c1 + c2 and ic1 ? ic2 , so we might

as well rename them c1 and c2 , in which case the solution would take the form

y(x) = c1 cos 2x + c2 sin 2x.

(1.2.16)

Here’s an example that shows the various possibilities:

y (8) ? 5y (7) + 17y (6) ? 997y (5) + 110y (4) ? 531y (3) + 765y (2) ? 567y 0 + 162y = 0. (1.2.17)

The equation was cooked up to have a characteristic polynomial that can be factored as

(? ? 2)(?2 + 9)2 (? ? 1)3 .

(1.2.18)

Hence the roots of the characteristic equation are 2 (simple), 3i (multiplicity 2), ?3i (multiplicity 2), and 1 (multiplicity 3).

1.3 Difference equations

11

Corresponding to the root 2, the general solution will contain the term c1 e2x . Corresponding to the double root at 3i we have terms (c2 + c3 x)e3ix in the solution. From the

double root at ?3i we get a contribution (c4 + c5 x)e?3ix , and finally from the triple root

at 1 we get (c6 + c7 x + c8 x2 )ex . The general solution is the sum of these eight terms.

Alternatively, we might have taken the four terms that come from 3i in the form

(c2 + c3 x) cos 3x + (c4 + c5 x) sin 3x.

(1.2.19)

Exercises 1.2

1. Obtain the general solutions of each of the following differential equations:

(a) y 00 + 5y 0 + 6y = 0

(b) y 00 ? 8y 0 + 7y = 0

(c) (D + 3)2 y = 0

(d) (D2 + 16)2 y = 0

(e) (D + 3)3 (D 2 ? 25)2 (D + 2)3 y = 0

2. Find a curve y = f (x) that passes through the origin with unit slope, and which

satisfies (D + 4)(D ? 1)y = 0.

1.3

Difference equations

Whereas a differential equation is an equation in an unknown function, a difference equation

is an equation in an unknown sequence. For example, suppose we know that a certain

sequence of numbers y0 , y1 , y2 , . . . satisfies the following conditions:

yn+2 + 5yn+1 + 6yn = 0

n = 0, 1, 2, . . .

(1.3.1)

and furthermore that y0 = 1 and y1 = 3.

Evidently, we can compute as many of the yn ’s as we need from (1.3.1), thus we would

get y2 = ?21, y3 = 87, y4 = ?309 so forth. The entire sequence of yn ’s is determined by

the difference equation (1.3.1) together with the two starting values.

Such equations are encountered when differential equations are solved on computers.

Naturally, the computer can provide the values of the unknown function only at a discrete

set of points. These values are computed by replacing the given differential equations by

a difference equation that approximates it, and then calculating successive approximate

values of the desired function from the difference equation.

Can we somehow “solve” a difference equation by obtaining a formula for the values

of the solution sequence? The answer is that we can, as long as the difference equation is

linear and has constant coefficients, as in (1.3.1). Just as in the case of differential equations

with constant coefficients, the correct strategy for solving them is to try a solution of the

12

Differential and Difference Equations

right form. In the previous section, the right form to try was y(x) = e?x . Now the winning

combination is y = ?n , where ? is a constant.

In fact, let’s substitute ?n for yn in (1.3.1) to see what happens. The left side becomes

?n+2 + 5?n+1 + 6?n = ?n (?2 + 5? + 6) = 0.

(1.3.2)

Just as we were able to cancel the common factor e?x in the differential equation case, so

here we can cancel the ?n , and we’re left with the quadratic equation

?2 + 5? + 6 = 0.

(1.3.3)

The two roots of this characteristic equation are ? = ?2 and ? = ?3. Therefore the

sequence (?2)n satisfies (1.3.1) and so does (?3)n . Since the difference equation is linear,

it follows that

yn = c1 (?2)n + c2 (?3)n

(1.3.4)

is also a solution, whatever the values of the constants c1 and c2 .

Now it is evident from (1.3.1) itself that the numbers yn are uniquely determined if we

prescribe the values of just two of them. Hence, it is very clear that when we have a solution

that contains two arbitrary constants we have the most general solution.

When we take account of the given data y0 = 1 and y1 = 3, we get the two equations

(

1 = c1 + c2

3 = (?2)c1 + (?3)c2

(1.3.5)

from which c1 = 6 and c2 = ?5. Finally, we use these values of c1 and c2 in (1.3.4) to get

yn = 6(?2)n ? 5(?3)n

n = 0, 1, 2, . . . .

(1.3.6)

Equation (1.3.6) is the desired formula that represents the unique solution of the given

difference equation together with the prescribed starting values.

Let’s step back a few paces to get a better view of the solution. Notice that the formula

(1.3.6) expresses the solution as a linear combination of nth powers of the roots of the

associated characteristic equation (1.3.3). When n is very large, is the number yn a large

number or a small one? Evidently the powers of ?3 overwhelm those of ?2, so the sequence

will behave roughly like a constant times powers of ?3. This means that we should expect

the members of the sequence to alternate in sign and to grow rapidly in magnitude.

So much for the equation (1.3.1). Now let’s look at the general case, in the form of a

linear difference equation of order p:

yn+p + a1 yn+p?1 + a2 yn+p?2 + · · · + ap yn = 0.

(1.3.7)

We try a solution of the form yn = ?n , and after substituting and canceling, we get the

characteristic equation

?p + a1 ?p?1 + a2 ?p?2 + · · · + ap = 0.

(1.3.8)

1.3 Difference equations

13

This is a polynomial equation of degree p, so it has p roots, counting multiplicities, somewhere in the complex plane.

Let ?? be one of these p roots. If ?? is simple (i.e., has multiplicity 1) then the part

of the general solution that corresponds to ?? is c(?? )n . If, however, ?? is a root of

multiplicity k > 1 then we must multiply the solution c(?? )n by an arbitrary polynomial

in n, of degree k ? 1, just as in the corresponding case for differential equations we used an

arbitrary polynomial in x of degree k ? 1.

We illustrate this, as well as the case of complex roots, by considering the following

difference equation of order five:

yn+5 ? 5yn+4 + 9yn+3 ? 9yn+2 + 8yn+1 ? 4yn = 0.

(1.3.9)

This example is rigged so that the characteristic equation can be factored as

(?2 + 1)(? ? 2)2 (? ? 1) = 0

(1.3.10)

from which the roots are obviously i, ?i, 2 (multiplicity 2), 1.

Corresponding to the roots i, ?i, the portion of the general solution is c1 in + c2 (?i)n .

Since

n?

n?

in = ein?/2 = cos

+ i sin

(1.3.11)

2

2

and similarly for (?i)n , we can also take this part of the general solution in the form

n?

c1 cos

2

n?

+ c2 sin

.

2

(1.3.12)

The double root ? = 2 contributes (c3 + c4 n)2n , and the simple root ? = 1 adds c5 to

the general solution, which in its full glory is

yn = c1 cos

n?

2

+ c2 sin

n?

2

+ (c3 + c4 n)2n + c5 .

(1.3.13)

The five constants would be determined by prescribing five initial values, say y0 , y1 , y2 , y3

and y4 , as we would expect for the equation (1.3.9).

Exercises 1.3

1. Obtain the general solution of each of the following difference equations:

(a) yn+1 = 3yn

(b) yn+1 = 3yn + 2

(c) yn+2 ? 2yn+1 + yn = 0

(d) yn+2 ? 8yn+1 + 12yn = 0

(e) yn+2 ? 6yn+1 + 9yn = 1

(f) yn+2 + yn = 0

14

Differential and Difference Equations

2. Find the solution of the given difference equation that takes the prescribed initial

values:

(a)

(b)

(c)

(d)

yn+2 = 2yn+1 + yn ; y0 = 0; y1 = 1

yn+1 = ?yn + ?; y0 = 1

yn+4 + yn = 0; y0 = 1; y1 = ?1; y2 = 1; y3 = ?1

yn+2 ? 5yn+1 + 6yn = 0; y0 = 1; y1 = 2

3. (a) For each of the difference equations in problems 1 and 2, evaluate

yn+1

lim

n?? yn

(1.3.14)

if it exists.

(b) Formulate and prove a general theorem about the existence of, and value of the

limit in part (a) for a linear difference equation with constant coefficients.

(c) Reverse the process: given a polynomial equation, find its root of largest absolute

value by computing from a certain difference equation and evaluating the ratios

of consecutive terms.

(d) Write a computer program to implement the method in part (c). Use it to

calculate the largest root of the equation

x8 = x7 + x6 + x5 + · · · + 1.

1.4

(1.3.15)

Computing with difference equations

This is, after all, a book about computing, so let’s begin with computing from difference

equations since they will give us a chance to discuss some important questions that concern

the design of computer programs. For a sample difference equation we’ll use

yn+3 = yn+2 + 5yn+1 + 3yn

(1.4.1)

together with the starting values y0 = y1 = y2 = 1. The reader might want, just for practice,

to find an explicit formula for this sequence by the methods of the previous section.

Suppose we want to compute a large number of these y’s in order to verify some property

that they have, for instance to check that

yn+1

lim

=3

(1.4.2)

n?? yn

which must be true since 3 is the root of largest absolute value of the characteristic equation.

As a first approach, we might declare y to be a linear array of some size large enough to

accommodate the expected length of the calculation. Then the rest is easy. For each n, we

would calculate the next yn+1 from (1.4.1), we would divide it by its predecessor yn to get

a new ratio. If the new ratio agrees sufficiently well with the previous ratio we announce

that the computation has terminated and print the new ratio as our answer. Otherwise, we

move the new ratio to the location of the old ratio, increase n and try again.

If we were to write this out as formal procedure (algorithm) it might look like:

1.4 Computing with difference equations

15

y0 := 1; y1 := 1; y2 := 1; n := 2;

newrat := ?10; oldrat := 1;

while |newrat ? oldrat| ? 0.000001 do

oldrat := newrat; n := n + 1;

yn := yn?1 + 5 ? yn?2 + 3 ? yn?3 ;

newrat := yn /yn?1

endwhile

print newrat; Halt.

We’ll use the symbol ‘:=’ to mean that we are to compute the quantity on the right, if

necessary, and then store it in the place named on the left. It can be read as ‘is replaced

by’ or ‘is assigned.’ Also, the block that begins with ‘while’ and ends with ‘endwhile’

represents a group of instructions that are to be executed repeatedly until the condition

that follows ‘while’ becomes false, at which point the line following ‘endwhile’ is executed.

The procedure just described is fast, but it uses lots of storage. If, for instance, such a

program needed to calculate 79 y’s before convergence occurred, then it would have used

79 locations of array storage. In fact, the problem above doesn’t need that many locations

because convergence happens a lot sooner. Suppose you wanted to find out how much

sooner, given only a programmable hand calculator with ten or twenty memory locations.

Then you might appreciate a calculation procedure that needs just four locations to hold

all necessary y’s.

That’s fairly easy to accomplish, though. At any given moment in the program, what

we need to find the next y are just the previous three y’s. So why not save only those three?

We’ll use the previous three to calculate the next one, and stow it for a moment in a fourth

location. Then we’ll compute the new ratio and compare it with the old. If they’re not

close enough, we move each one of the three newest y’s back one step into the places where

we store the latest three y’s and repeat the process. Formally, it might be:

y := 1; ym1 := 1; ym2 := 1;

newrat := ?10; oldrat := 1;

while |newrat ? oldrat| ? 0.000001 do

ym3 := ym2; ym2 := ym1; ym1 := y;

oldrat := newrat;

y := ym1 + 5 ? ym2 + 3 ? ym3;

newrat := y/ym1 endwhile;

print newrat; Halt.

The calculation can now be done in exactly six memory locations (y, ym1, ym2, ym3,

oldrat, newrat) no matter how many y’s have to be calculated, so you can undertake it

on your hand calculator with complete confidence. The price that we pay for the memory

saving is that we must move the data around a bit more.

One should not think that such programming methods are only for hand calculators. As

we progress through the numerical solution of differential equations we will see situations

in which each of the quantities that appears in the difference equation will itself be an

16

Differential and Difference Equations

array (!), and that very large numbers, perhaps thousands, of these arrays will need to be

computed. Even large computers might quake at the thought of using the first method

above, rather than the second, for doing the calculation. Fortunately, it will almost never

be necessary to save in memory all of the computed values simultaneously. Normally, they

will be computed, and then printed or plotted, and never needed except in the calculation

of their immediate successors.

Exercises 1.4

1. The Fibonacci numbers F0 , F1 , F2 , . . . are defined by the recurrence formula Fn+2 =

Fn+1 + Fn for n = 0, 1, 2, . . . together with the starting values F0 = 0, F1 = 1.

(a) Write out the first ten Fibonacci numbers.

(b) Derive an explicit formula for the nth Fibonacci number Fn .

(c) Evaluate your formula for n = 0, 1, 2, 3, 4.

(d) Prove directly from your formula that the Fibonacci numbers are integers (This is

perfectly obvious from their definition, but is not so obvious from the formula!).

(e) Evaluate

lim

n??

Fn+1

Fn

(1.4.3)

(f) Write a computer program that will compute Fibonacci numbers and print out

the limit in part (e) above, correct to six decimal places.

(g) Write a computer program that will compute the first?

40 members of the modified

Fibonacci sequence in which F0 = 1 and F1 = (1 ? 5)/2. Do these computed

numbers seem to be approaching zero? Explain carefully what you see and why

it happens.

(h) Modify the program of part (h) to run in higher (or double) precision arithmetic.

Does it change any of your answers?

2. Find the most general solution of each of the following difference equations:

(a) yn+1 ? 2yn + yn?1 = 0

(b) yn+1 = 2yn

(c) yn+2 + yn = 0

(d) yn+2 + 3yn+1 + 3yn + yn?1 = 0

1.5

Stability theory

In the study of natural phenomena it is most often true that a small change in conditions

will produce just a small change in the state of the system being studied. If, for example,

a very slight increase in atmospheric pollution could produce dramatically large changes

in populations of flora and fauna, or if tiny variations in the period of the earth’s rotation

1.5 Stability theory

17

produced huge changes in climatic conditions, the world would be a very different place to

live in, or to try to live in. In brief, we may say that most aspects of nature are stable.

When physical scientists attempt to understand some facet of nature, they often will

make a mathematical model. This model will usually not faithfully reproduce all of the

structure of the original phenomenon, but one hopes that the important features of the

system will be preserved in the model, so that predictions will be possible. One of the most

important features to preserve is that of stability.

For instance, the example of atmospheric pollution and its effect on living things referred

to above is important and very complex. Therefore considerable effort has gone into the

construction of mathematical models that will allow computer studies of the effects of

atmospheric changes. One of the first tests to which such a model should be subjected is

that of stability: does it faithfully reproduce the observed fact that small changes produce

small changes? What is true in nature need not be true in a man-made model that is a

simplification or idealization of the real world.

Now suppose that we have gotten ourselves over this hurdle, and we have constructed

a model that is indeed stable. The next step might be to go to the computer and do

calculations from the model, and to use these calculations for predicting the effects of

various proposed actions. Unfortunately, yet another layer of approximation is usually

introduced at this stage, because the model, even though it is a simplification of the real

world, might still be too complicated to solve exactly on a computer.

For instance, may models use differential equations. Models of the weather, of the motion of fluids, of the movement of astronomical objects, of spacecraft, of population growth,

of predator-prey relationships, of electric circuit transients, and so forth, all involve differential equations. Digital computers solve differential equations by approximating them by

difference equations, and then solving the difference equations. Even though the differential

equation that represents our model is indeed stable, it may be that the difference equation

that we use on the computer is no longer stable, and that small changes in initial data on

the computer, or small roundoff errors, will produce not small but very large changes in the

computed solution.

An important job of the numerical analyst is to make sure that this does not happen,

and we will find that this theme of stability recurs throughout our study of computer

approximations.

As an example of instability in differential equations, suppose that some model of a

system led us to the equation

y 00 ? y 0 ? 2y = 0

(1.5.1)

together with the initial data

y(0) = 1;

y 0 (0) = ?1.

(1.5.2)

We are thinking of the independent variable t as the time, so we will be interested in

the solution as t becomes large and positive.

The general solution of (1.5.1) is y(t) = c1 e?t + c2 e2t . The initial conditions tell us that

c1 = 1 and c2 = 0, hence the solution of our problem is y(t) = e?t , and it represents a

18

Differential and Difference Equations

function that decays rapidly to zero with increasing t. In fact, when t = 10, the solution

has the value 0.000045.

Now let’s change the initial data (1.5.2) just a bit, by asking for a solution with y 0 (0) =

?0.999. It’s easy to check that the solution is now

y(t) = (0.999666 . . .)e?t + (0.000333 . . .)e2t

(1.5.3)

instead of just y(t) = e?t . If we want the value of the solution at t = 10, we would find

that it has changed from 0.000045 to about 7.34.

At t = 20 the change is even more impressive, from 0.00000002 to 161, 720+, just from

changing the initial value of y 0 from ?1 to ?0.999. Let’s hope that there are no phenomena

in nature that behave in this way, or our lives hang by a slender thread indeed!

Now exactly what is the reason for the observed instability of the equation (1.5.1)?

The general solution of the equation contains a falling exponential term c1 e?t , and a rising

exponential term c2 e2t . By prescribing the initial data (1.5.2) we suppressed the growing

term, and picked out only the decreasing one. A small change in the initial data, however,

results in the presence of both terms in the solution.

Now it’s time for a formal

Definition: A differential equation is said to be stable if for every set of initial data (at

t = 0) the solution of the differential equation remains bounded as t approaches infinity.

A differential equation is called strongly stable if, for every set of initial data (at t = 0)

the solution not only remains bounded, but approaches zero as t approaches infinity.

What makes the equation (1.5.1) unstable, then, is the presence of a rising exponential

in its general solution. In other words, if we have a differential equation whose general

solution contains a term e?t in which ? is positive, that equation is unstable.

Let’s restrict attention now to linear differential equations with constant coefficients.

We know from section 1.2 that the general solution of such an equation is a sum of terms

of the form

(polynomial in t)e?t .

(1.5.4)

Under what circumstances does such a term remain bounded as t becomes large and positive?

Certainly if ? is negative then the term stays bounded. Likewise, if ? is a complex

number and its real part is negative, then the term remains bounded. If ? has positive real

part the term is unbounded.

This takes care of all of the possibilities except the case where ? is zero, or more generally,

the complex number ? has zero real part (is purely imaginary). In that case the question

of whether (polynomial in t)e?t remains bounded depend on whether the “polynomial in t”

is of degree zero (a constant polynomial) or of higher degree. If the polynomial is constant

then the term does indeed remain bounded for large positive t, whereas otherwise the term

will grow as t gets large, for some values of the initial conditions, thereby violating the

definition of stability.

1.6 Stability theory of difference equations

19

Now recall that the “polynomial in t” is in fact a constant if the root ? is a simple

root of the characteristic equation of the differential equation, and otherwise it is of higher

degree. This observation completes the proof of the following:

Theorem 1.5.1 A linear differential equation with constant coefficients is stable if and

only if all of the roots of its characteristic equation lie in the left half plane, and those that

lie on the imaginary axis, if any, are simple. Such an equation is strongly stable if and only

if all of the roots of its characteristic equation lie in the left half plane, and none lie on the

imaginary axis.

Exercises 1.5

1. Determine for each of the following differential equations whether it is strongly stable,

stable, or unstable.

(a) y 00 ? 5y 0 + 6y = 0

(b) y 00 + 5y 0 + 6y = 0

(c) y 00 + 3y = 0

(d) (D + 3)3 (D + 1)y = 0

(e) (D + 1)2 (D 2 + 1)2 y = 0

(f) (D4 + 1)y = 0

2. Make a list of some natural phenomena that you think are unstable. Discuss.

3. The differential equation y 00 ?y = 0 is to be solved with the initial conditions y(0) = 1,

y 0 (0) = ?1, and then solved again with y(0) = 1, y 0 (0) = ?0.99. Compare the two

solutions when x = 20.

4. For exactly which real values of the parameter ? is each of the following differential

equations stable? . . . strongly stable?

(a) y 00 + (2 + ?)y 0 + y = 0

(b) y 00 + ?y 0 + y = 0

(c) y 0 + ?y = 1

1.6

Stability theory of difference equations

In the previous section we discussed the stability of differential equations. The key ideas

were that such an equation is stable if every one of its solutions remains bounded as t

approaches infinity, and strongly stable if the solutions actually approach zero.

Similar considerations apply to difference equations, and for similar reasons. As an

example, take the equation

5

yn+1 = yn ? yn?1

2

(n ? 1)

(1.6.1)

20

Differential and Difference Equations

along with the initial equations

y0 = 1;

y1 = 0.5 .

(1.6.2)

It’s easy to see that the solution is yn = 2?n , and of course, this is a function that

rapidly approaches zero with increasing n.

Now let’s change the initial data (1.6.2), say to

y0 = 1;

y1 = 0.50000001

(1.6.3)

instead of (1.6.2).

The solution of the difference equation with these new data is

y = (0.0000000066 . . .)2n + (0.9999999933 . . .)2?n .

(1.6.4)

The point is that the coefficient of the growing term 2n is small, but 2n grows so fast that

after a while the first term in (1.6.4) will be dominant. For example, when n = 30, the

solution is y30 = 7.16, compared to the value y30 = 0.0000000009 of the solution with the

original initial data (1.6.2). A change of one part in fifty million in the initial condition

produced, thirty steps later, an answer one billion times as large.

The fault lies with the difference equation, because it has both rising and falling components to its general solution. It should be clear that it is hopeless to do extended computation with an unstable difference equation, since a small roundoff error may alter the

solution beyond recognition several steps later.

As in the case of differential equations, we’ll say that a difference equation is stable

if every solution remains bounded as n grows large, and that it is strongly stable if every

solution approaches zero as n grows large. Again, we emphasize that every solution must

be well behaved, not just the solution that is picked out by a certain set of initial data. In

other words, the stability, or lack of it, is a property of the equation and not of the starting

values.

Now consider the case where the difference equation is linear with constant coefficients.

The we know that the general solution is a sum of terms of the form

(polynomial in n)?n .

(1.6.5)

Under what circumstances will such a term remain bounded or approach zero?

Suppose |?| < 1. Then the powers of ? approach zero, and multiplication by a polynomial in n does not alter that conclusion. Suppose |?| > 1. Then the sequence of powers

grows unboundedly, and multiplication by a nonzero polynomial only speeds the parting

guest.

Finally suppose the complex number ? has absolute value 1. Then the sequence of its

powers remains bounded (in fact they all have absolute value 1), but if we multiply by a

nonconstant polynomial, the resulting expression would grow without bound.

To summarize then, the term (1.6.5), if the polynomial is not identically zero, approaches

zero with increasing n if and only if |?| < 1. It remains bounded as n increases if and only

1.6 Stability theory of difference equations

21

if either (a) |?| < 1 or (b) |?| = 1 and the polynomial is of degree zero (a constant). Now

we have proved:

Theorem 1.6.1 A linear difference equation with constant coefficients is stable if and only

if all of the roots of its characteristic equation have absolute value at most 1, and those of

absolute value 1 are simple. The equation is strongly stable if and only if all of the roots

have absolute value strictly less than 1.

Exercises 1.6

1. Determine, for each of the following difference equations whether it is strongly stable,

stable, or unstable.

(a) yn+2 ? 5yn+1 + 6yn = 0

(b) 8yn+2 + 2yn+1 ? 3yn = 0

(c) 3yn+2 + yn = 0

(d) 3yn+3 + 9yn+2 ? yn+1 ? 3yn = 0

(e) 4yn+4 + 5yn+2 + yn = 0

2. The difference equation 2yn+2 + 3yn+1 ? 2yn = 0 is to be solved with the initial

conditions y0 = 2, y1 = 1, and then solved again with y0 = 2, y1 = 0.99. Compare y20

for the two solutions.

3. For exactly which real values of the parameter ? is each of the following difference

equations stable? . . . strongly stable?

(a) yn+2 + ?yn+1 + yn = 0

(b) yn+1 + ?yn = 1

(c) yn+2 + yn+1 + ?yn = 0

4. (a) Consider the (constant-coefficient) difference equation

a0 yn+p + a1 yn+p?1 + a2 yn+p?2 + · · · + ap yn = 0.

(1.6.6)

Show that this difference equation cannot be stable if |ap /a0 | > 1.

(b) Give an example to show that the converse of the statement in part (a) is false.

Namely, exhibit a difference equation for which |ap /a0 | < 1 but the equation is unstable anyway.

22

Differential and Difference Equations

Chapter 2

The Numerical Solution of

Differential Equations

2.1

Euler’s method

Our study of numerical methods will begin with a very simple procedure, due to Euler.

We will state it as a method for solving a single differential equation of first order. One of

the nice features of the subject of numerical integration of differential equations is that the

techniques that are developed for just one first order differential equation will apply, with

very little change, both to systems of simultaneous first order equations and to equations of

higher order. Hence the consideration of a single equation of first order, seemingly a very

special case, turns out to be quite general.

By an initial-value problem we mean a differential equation together with enough given

values of the unknown function and its derivatives at an initial point x0 to determine the

solution uniquely.

Let’s suppose that we are given an initial-value problem of the form

y 0 = f (x, y);

y(x0 ) = y0 .

(2.1.1)

Our job is to find numerical approximate values of the unknown function y at points x

to the right of (larger than) x0 .

What we actually will find will be approximate values of the unknown function at a

discrete set of points x0 , x1 = x0 + h, x2 = x0 + 2h, x3 = x0 + 3h, etc. At each of these

points xn we will compute yn , our approximation to y(xn ).

Hence, suppose that the spacing h between consecutive points has been chosen. We

propose to start at the point x0 where the initial data are given, and move to the right,

obtaining y1 from y0 , then y2 from y1 and so forth until sufficiently many values have been

found.

Next we need to derive a method by which each value of y can be obtained from its

immediate predecessor. Consider the Taylor series expansion of the unknown function y(x)

24

The Numerical Solution of Differential Equations

about the point xn

y 00 (X)

,

(2.1.2)

2

where we have halted the expansion after the first power of h and in the remainder term,

the point X lies between xn and xn + h.

y(xn + h) = y(xn ) + hy 0 (xn ) + h2

Now equation (2.1.2) is exact, but of course it cannot be used for computation because

the point X is unknown. On the other hand, if we simply “forget” the error term, we’ll

have only an approximate relation instead of an exact one, with the consolation that we

will be able to compute from it. The approximate relation is

y(xn + h) ? y(xn ) + hy 0 (xn ).

(2.1.3)

Next define yn+1 to be the approximate value of y(xn+1 ) that we obtain by using the

right side of (2.1.3) instead of (2.1.2). Then we get

yn+1 = yn + hyn0 .

(2.1.4)

Now we have a computable formula for the approximate values of the unknown function,

because the quantity yn0 can be found from the differential equation (2.1.1) by writing

yn0 = f (xx , yn ),

(2.1.5)

and if we do so then (2.1.4) takes the form

yn+1 = yn + hf (xn , yn ).

(2.1.6)

This is Euler’s method, in a very explicit form, so that the computational procedure is

clear. Equation (2.1.6) is in fact a recurrence relation, or difference equation, whereby each

value of yn is computed from its immediate predecessor.

Let’s use Euler’s method to obtain a numerical solution of the differential equation

y 0 = 0.5y

(2.1.7)

together with the starting value y(0) = 1. The exact solution of this initial-value problem

is obviously y(x) = ex/2 .

Concerning the approximate solution by Euler’s method, we have, by comparing (2.1.7)

with (2.1.1), f (x, y) = 0.5y, so

yn

2

h

= 1+

yn .

2

yn+1 = yn + h

(2.1.8)

Therefore, in this example, each yn will be obtained from its predecessor by multiplication

by 1 + h2 . To be quite specific, let’s take h to be 0.05. Then we show below, for each value

of x = 0, 0.05, 0.10, 0.15, 0.20, . . . the approximate value of y computed from (2.1.8) and

the exact value y(xn ) = exn /2 :

2.1 Euler’s method

25

x

0.00

0.05

0.10

0.15

0.20

0.25

..

.

Euler(x)

1.00000

1.02500

1.05063

1.07689

1.10381

1.13141

..

.

Exact(x)

1.00000

1.02532

1.05127

1.07788

1.10517

1.13315

..

.

1.00

2.00

3.00

..

.

1.63862

2.68506

4.39979

..

.

1.64872

2.71828

4.48169

..

.

5.00

10.00

11.81372

139.56389

12.18249

148.41316

table 1

Considering the extreme simplicity of the approximation, it seems that we have done

pretty well by this equation. Let’s continue with this simple example by asking for a formula

for the numbers that are called Euler(x) in the above table. In other words, exactly what

function of x is Euler(x)?

To answer this, we note first that each computed value yn+1 is obtained according to

(2.1.8) by multiplying its predecessor yn by 1 + h2 . Since y0 = 1, it is clear tha we will

compute yn = (1 + h2 )n . Now we want to express this in terms of x rather than n. Since

xn = nh, we have n = x/h, and since h = 0.05 we have n = 20x. Hence the computed

approximation to y at a particular point x will be 1.02520x , or equivalently

Euler(x) = (1.638616 . . .)x .

(2.1.9)

The approximate values can now easily be compared with the true solution, since

x

1

Exact(x) = e 2 = e 2

x

= (1.648721 . . .)x .

(2.1.10)

Therefore both the exact solution of this differential equation and its computed solution

have the form (const.)x . The correct value of “const.” is e1/2 , and the value that is, in

effect, used by Euler’s method is (1 + h2 )1/h . For a fixed value of x, we see that if we use

Euler’s method with smaller and smaller values of h (neglecting the increase in roundoff

error that is sure to result), the values Euler(x) will converge to Exact(x), because

h

lim 1 +

h?0

2

Exercises 2.1

1/h

1

= e2 .

(2.1.11)

26

The Numerical Solution of Differential Equations

1. Verify the limit (2.1.11).

2. Use a calculator or a computer to integrate each of the following differential equations

forward ten steps, using a spacing h = 0.05 with Euler’s method. Also tabulate the

exact solution at each value of x that occurs.

(a) y 0 (x) = xy(x); y(0) = 1

(b) y 0 (x) = xy(x) + 2; y(0) = 1

y(x)

(c) y 0 (x) =

; y(0) = 1

1+x

(d) y 0 (x) = ?2xy(x)2 ; y(0) = 10

2.2

Software notes

One of the main themes of our study will be the preparation of programs that not only

work, but also are easily readable and useable by other people. The act of communication

that must take place before a program can be used by persons other than its author is a

difficult one to carry out, and we will return several times to the principles that have evolved

as guides to the preparation of readable software. Here are some of these guidelines.

1. Documentation

The documentation of a program is the set of written instructions to a user that inform

the user about the purpose and operation of the program. At the moment that the job

of writing and testing a program has been completed it is only natural to feel an urge to

get the whole thing over with and get on to the next job. Besides, one might think, it’s

perfectly obvious how to use this program. Some programs may be obscure, but not this

one.

It is amazing how rapidly our knowledge of our very own program fades. If we come

back to a program after a lapse of a few months’ time, it often happens that we will have no

idea what the program did or how to use it, at least not without making a large investment

of time.

For that reason it is important that when our program has been written and tested it

should be documented immediately, while our memory of it is still green. Furthermore, the

best place for documentation is in the program itself, in “comment” statements. That way

one can be sure that when the comments are needed they will be available.

The first mission of program documentation is to describe the purpose of the program.

State clearly the problem that the program solves, or the exact operation that it performs

on its input in order to get its output.

Already in this first mission, a good bit of technical skill can be brought to bear that

will be very helpful to the use, by intertwining the description of the program purpose with

the names of the communicating variables in the program.

Let’s see what that means by considering an example. Suppose we have written a

subroutine that searches through a specified row of a matrix to find the element of largest

2.2 Software notes

27

absolute value, and outputs a column in which it was found. Such a routine, in Maple for

instance, might look like this:

search:=proc(A,i)

local j, winner, jwin;

winner:=-1;

for j from 1 to coldim(A) do

if (abs(A[i,j])>winner) then

winner:=abs(A[i,j]) ; jwin:=j fi

od;

return(jwin);

end;

Now let’s try our hand at documenting this program:

“The purpose of this program is to search a given row of a matrix to find an

element of largest absolute value and return the column in which it was found.”

That is pretty good documentation, perhaps better than many programs get. But we

can make it a lot more useful by doing the intertwining that we referred to above. There

we said that the description should be related to the communicating variables. Those

variables are the ones that the user can see. They are the input and output variables of

the subroutine. In most important computer languages, the communicating variables are

announced in the first line of the coding of a procedure or subroutine. Maple follows this

convention at least for the input variables, although the output variable is usually specified

in the “return” statement.

In the first line of the little subroutine above we see the list (A, i) of its input variables

(or “arguments”). These are the ones that the user has to understand, as opposed to the

other “local” variables that live inside the subroutine but don’t communicate with the

outside world (like j, winner, jwin, which are listed on the second line of the program).

The best way to help the user to understand these variables is to relate them directly

to the description of the purpose of the program.

“The purpose of this program is to search row I of a given matrix A to find an

entry of largest absolute value, and returns the column jwin where that entry

lives.”

We’ll come back to the subject of documentation, but now let’s mention another ingredient of ease of use of programs, and that is:

2. Modularity

It is important to divide a long program into a number of smaller modules, each with a

clearly stated set of inputs and outputs, and each with its own documentation. That means

that we should get into the habit of writing lots of subroutines or procedures, because

28

The Numerical Solution of Differential Equations

the subroutine or procedure mode of expression forces one to be quite explicit about the

relationship of the block of coding to the rest of the world.

When we are writing a large program we would all write a subroutine if we found that

a certain sequence of steps was being called for repeatedly. Beyond this, however, there are

numerous inducements for breaking off subroutines even if the block of coding occurs just

once in the main program.

For one thing it’s easier to check out the program. The testing procedure would consist

of first testing each of the subroutines separately on test problems designed just for them.

Once the subroutines work, it would remain only to test their relationships to the calling

program.

For another reason, we might discover a better, faster, more elegant, or what-have-you

method of performing the task that one of these subroutines does. Then we would be able

to yank out the former subroutine and plug in the new one, while being careful only to make

sure that the new subroutine relates to the same inputs and outputs as the old one. If jobs

within a large program are not broken into subroutines it can be much harder to isolate the

block of coding that deals with a particular function and remove it without affecting the

whole works.

For another reason, if one be needed, it may well happen that even though the job that

is done by the subroutine occurs only once in the current program, it may recur in other

programs as yet undreamed of. If one is in the habit of writing small independent modules

and stringing them together to make large programs, then it doesn’t take long before one

has a library of useful subroutines, each one tested, working and documented, that will

greatly simplify the task of writing future programs.

Finally, the practice of subdividing the large jobs into the smaller jobs of which they are

composed is an extremely valuable analytical skill, one that is useful not only in programming, but in all sorts of organizational activities where smaller efforts are to be pooled in

order to produce a larger effect. It is therefore a quality of mind that provides much of its

own justification.

In this book, the major programs that are the objects of study have been broken up into

subroutines in the expectation that the reader will be able to start writing and checking out

these modules even before the main ideas of the current subject have been fully explained.

This was done in part because some of these programs are quite complex, and it would be

unreasonable to expect the whole program to be written in a short time. It was also done

to give examples of the process of subdivision that we have been talking about.

For instance, the general linear algebra program for solving systems of linear simultaneous equations in Chapter 3, has been divided into six modules, and they are described in

section 3.3. The reader might wish to look ahead at those routines and to verify that even

though their relationship to the whole job of solving equations is by no means clear now,

nonetheless, because of the fact that they are independent and self-contained, they can be

programmed and checked out at any time without waiting for the full explanation.

One more ingredient that is needed for the production of useful software is:

2.3 Systems and equations of higher order

29

3. Style

We don’t mean style in the sense of “class,” although this is as welcome in programming

as it is elsewhere. There have evolved a number of elements of good programming style,

and these will mainly be discussed as they arise. But two of them (one trivial and one quite

deep) are:

(a) Indentation: The instructions that lie within the range of a loop are indented in the

program listing further to the right than the instructions that announce that the loop

is about to begin, or that it has just terminated.

(b) Top-down structuring: When we visualize the overall logical structure of a complicated program we see a grand loop, within which there are several other loops and

branchings, within which . . . etc. According to the principles of top-down design the

looping and branching structure of the program should be visible at once in the listing. That is, we should see an announcement of the opening of the grand loop, then

indented under that perhaps a two-way branch (if-then-else), where, under the “then”

one sees all that will happen if the condition is met, and under the “else” one sees

what happens if it is not met.

When we say that we see all that will happen, we mean that there are not any “go-to”

instructions that would take our eye out of the flow of the if-then-else loop to some

other page. It all happens right there on the same page, under “then” and under

“else”.

These few words can scarcely convey the ideas of structuring, which we leave to the

numerous examples in the sequel.

2.3

Systems and equations of higher order

We have already remarked that the methods of numerical integration for a single firstorder differential equation carry over with very little change to systems of simultaneous

differential equations of first order. In this section we’ll discuss exactly how this is done,

and furthermore, how the same idea can be applied to equations of higher order than the

first. Euler’s method will be used as the example, but the same transformations will apply

to all of the methods that we will study.

In Euler’s method for one equation, the approximate value of the unknown function at

the next point xn+1 = xn + h is calculated from

yn+1 = yn + hf (xn , yn ).

(2.3.1)

Now suppose that we are trying to solve not just a single equation, but a system of N

simultaneous equations, say

yi0 (x) = fi (x, y1 , y2 , . . . yN )

i = 1, . . . , N.

(2.3.2)

30

The Numerical Solution of Differential Equations

Equation (2.3.2) represents N equations, in each of which just one derivative appears,

and whose right-hand side may depend on x, and on all of the unknown functions, but not

on their derivatives. The “fi ” indicates that, of course, each equation can have a different

right-hand side.

Now introduce the vector y(x) of unknown functions

y(x) = [y1 (x), y2 (x), y3 (x), . . . , yN (x)]

(2.3.3)

and the vector f = f (x, y(x)) of right-hand sides

f = [f1 (x, y), f2 (x, y), . . . , fN (x, y)].

(2.3.4)

In terms of these vectors, equation (2.3.2) can be rewritten as

y 0 (x) = f (x, y(x)).

(2.3.5)

We observe that equation (2.3.5) looks just like our standard form (2.1.1) for a single

equation in a single unknown function, except for the bold face type, i.e., except for the

fact that y and f now represent vector quantities.

To apply a numerical method such as that of Euler, then, all we need to do is to take

the statement of the method for a single differential equation in a single unknown function,

and replace y(x) and f (x, y(x)) by vector quantities as above. We will then have obtained

the generalization of the numerical method to systems.

To be specific, Euler’s method for a single equation is

yn+1 = yn + hf (x, yn )

(2.3.6)

so Euler’s method for a system of differential equations will be

y n+1 = y n + hf (xn , y n ).

(2.3.7)

This means that if we know the entire vector y of unknown functions at the point x = xn ,

then we can find the entire vector of unknown functions at the next point xn+1 = xn + h

by means of (2.3.7).

In detail, if y i (xn ) denotes the computed approximate value of the unknown function yi

at the point xn , then what we must calculate are

y i (xn+1 ) = y i (xn ) + hfi (xn , y 1 (xn ), y 2 (xn ), . . . , y N (xn ))

(2.3.8)

for each i = 1, 2, . . . , N .

As an example, take the pair of differential equations

(

y10 = x + y1 + y2

y20 = y1 y2 + 1

together with the initial values y1 (0) = 0, y2 (0) = 1.

(2.3.9)

2.3 Systems and equations of higher order

31

Now the vector of unknown functions is y = [y1 , y2 ], and the vector of right-hand sides

is f = [x + y1 + y2 , y1 y2 + 1]. Initially, the vector of unknowns is y = [0, 1]. Let’s choose a

step size h = 0.05. Then we calculate

"

and

"

y 1 (0.10)

y 2 (0.10)

y1 (0.05)

y2 (0.05)

#

"

=

0.05

1.05

#

"

=

0

1

#

#

+ 0.05

"

+ 0.05

"

0+0+1

0?1+1

0.05 + 0.05 + 1.05

0.05 ? 1.05 + 1

#

"

=

#

"

=

0.05

1.05

#

0.1075

1.102625

(2.3.10)

#

(2.3.11)

and so forth. At each step we compute the vector of approximate values of the two unknown

functions from the corresponding vector at the immediately preceding step.

Let’s consider the preparation of a computer program that will carry out the solution,

by Euler’s method, of a system of N simultaneous equations of the form (2.3.2), which we

will rewrite just slightly, in the form

yi0 = fi (x, y)

i = 1, . . . , N.

(2.3.12)

Note that on the left is just one of the unknown functions, and on the right there may

appear all N of them in each equation.

Evidently we will need an array Y of length N to hold the values of the N unknown

functions at the current point x. Suppose we have computed the array Y at a point x, and

we want to get the new array Y at the point x + h. Exactly what do we do?

Some care is necessary in answering this question because there is a bit of a snare in

the underbrush. The new values of the N unknown functions are calculated from (2.3.8) or

(2.3.12) in a certain order. For instance, we might calculate y 1 (x + h), then y 2 (x + h), then

y 3 (x + h),. . . , then yN (x + h).

The question is this: when we compute y 1 (x+h), where shall we put it? If we put it into

Y[1], the first position of the Y array in storage, then the previous contents of Y[1] are

lost, i.e., the value of y1 (x) is lost. But we aren’t finished with y 1 (x) yet; it’s still needed to

compute y 2 (x+h), y 3 (x+h), etc. This is because the new value y i (x+h) depends (or might

depend), according to (2.3.8), on the old values of all of the unknown functions, including

those whose new values have already been computed before we begin the computation of

y i (x + h).

If the point still is murky, go back to (2.3.11) and notice how, in the calculation of

y 2 (0.10) we needed to know y 1 (0.05) even though y1 (0.10) had already been computed.

Hence if we had put y 1 (0.10) into an array to replace the old value y 1 (0.05) we would not

have been able to obtain y 2 (0.10).

The conclusion is that we need at least two arrays, say YIN and YOUT, each of length N .

The array YIN holds the unknown functions evaluated at x, and YOUT will hold their values

at x + h. Initially YIN holds the given data at x0 . Then we compute all of the unknowns at

x0 + h, and store them in YOUT as we find them. When all have been done, we print them

if desired, move all entries of YOUT back to YIN, increase x by h and repeat.

32

The Numerical Solution of Differential Equations

The principal building block in this structure would be a subroutine that would advance

the solution exactly one step. The main program would initialize the arrays, call this

subroutine, increase x, move date from the output array YOUT back to the input array YIN,

print, etc. The single-step subroutine is shown below. We will use it later on to help us

get a solution started when we use methods that need information at more than one point

before they can start the integration.

Eulerstep:=proc(xin,yin,h,n) local i,yout;

# This program numerically integrates the system

# y’=f(x,y) one step forward by Euler’s method using step size

# h. Enter with values at xin in yin. Exit with values at xin+h

# in yout. Supply f as a function subprogram.

yout:=[seq(evalf(yin[i]+h*f(xin,yin,i)),i=1..n)];

return(yout);

end;

A few remarks about the program are in order. One structured data type in Maple

is a list of things, in this case, a list of floating point numbers. The seq command (for

“sequence”) creates such a list, in this case a list of length n since i goes from 1 to n in the

seq command. The brackets [ and ] convert the list into a vector. The evalf command

ensures that the results of the computation of the components of yout are floating point

numbers.

Our next remark about the program concerns the function subprogram f, which calculates the right hand sides of the differential equation. This subprogram, of course, must be

supplied by the user. Here is a sample of such a program, namely the one that describes

the system (2.3.9). In that case we have f1 (x, y) = x + y1 + y2 and f2 (x, y) = y1 y2 + 1. This

translates into the following:

f:=proc(x,y,i);

# Calculates the right-hand sides of the system of differential

# equations.

if i=1 then return(x+y[1]+y[2]) else return(y[1]*y[2]+1) fi;

end;

Our last comment about the program to solve systems is that it is perfectly possible to

use it in such a way that we would not have to move the contents of the vector YOUT back

to the vector YIN at each step. In other words, we could save N move operations, where N

is the number of equations. Such savings might be significant in an extended calculation.

To achieve this saving, we write two blocks of programming in the main program. One

block takes the contents of YIN as input, advances the solution one step by Euler’s method

and prints, leaving the new vector in YOUT. Then, without moving anything, another block

of programming takes the contents of YOUT as input, advances the solution one step, leaves

the new vector in YIN, and prints. The two blocks call Euler alternately as the integration

proceeds to the right. The reader might enjoy writing this program, and thinking about

how to generalize the idea to the situation where the new value of y is computed from two

2.3 Systems and equations of higher order

33

previously computed values, rather than from just one (then three blocks of programming

would be needed).

Now we’ve discussed the numerical solution of a single differential equation of first order,

and of a system of simultaneous differential equations of first order, and there remains the

treatment of equations of higher order than the first. Fortunately, this case is very easily

reduced to the varieties that we have already studied.

For example, suppose we want to solve a single equation of the second order, say

y 00 + xy 0 + (x + 1) cos y = 2.

(2.3.13)

The strategy is to transform the single second-order equation into a pair of simultaneous

first order equations that can then be handled as before. To do this, choose two unknown

functions u and v. The function u is to be the unknown function y in (2.3.13), and the

function v is to be the derivative of u. Then u and v satisfy two simultaneous first-order

differential equations:

(

u0 = v

(2.3.14)

v 0 = ?xv ? (x + 1) cos u + 2

and these are exactly of the form (2.3.5) that we have already discussed!

The same trick works on a general differential equation of N th order

y (N ) + G(x, y, y 0 , y 00 , . . . , y (N ?1) ) = 0.

(2.3.15)

We introduce N unknown functions u0 , u1 , . . . , uN ?1 , and let them be the solutions of the

system of N simultaneous first order equations

u00 = u1

u01 = u2

...

u0N ?2

u0N ?1

(2.3.16)

= uN ?1

= ?G(x, u0 , u1 , . . . , uN ?2 ).

The system can now be dealt with as before.

Exercises 2.3

1. Write each of the following as a system of simultaneous first-order initial-value problems in the standard form (2.3.2):

(a) y 00 + x2 y = 0; y(0) = 1; y 0 (0) = 0

(b) u0 + xv = 2; v 0 + euv = 0; u(0) = 0; v(0) = 0

(c) u0 + xv 0 = 0; v 0 + x2 u = 1; u(1) = 1; v(1) = 0

(d) y iv + 3xy 000 + x2 y 00 + 2y 0 + y = 0; y(0) = y 0 (0) = y 00 (0) = y 000 (0) = 1

(e) x000 (t) + t3 x(t) + y(t) = 0; y 00 (t) + x(t)2 = t3 ; x(0) = 2; x0 (0) = 1; x00 (0) = 0;

y(0) = 1; y 0 (0) = 0

34

The Numerical Solution of Differential Equations

2. For each of the parts of problem 1, write the function subprogram that will compute

the right-hand sides, as required by the Eulerstep subroutine.

3. For each of the parts of problem 1, assemble nad run on the computer the Euler

program, together with the relevant function subprogram of problem 2, to print out

the solutions for fifty steps of integration, each of size h = 0.03. Begin with x = x0 ,

the point at which the initial data was given.

4. Reprogram the Eulerstep subroutine, as discussed in the text, to avoid the movement

of YOUT back to YIN.

5. Modify your program as necessary (in Maple, take advantage of the plot command)

to produce graphical output (graph all of the unknown functions on the same axes).

Test your program by running it with Euler as it solves y 00 + y = 0 with y(0) = 0,

y 0 (0) = 1, and h = ?/60 for 150 steps.

6. Write a program that will compute successive values yp , yp+1 , . . . from a difference

equation of order p. Do this by storing the y’s as they are computed in a circular list,

so that it is never necessary to move back the last p computed values before finding

the next one. Write your program so that it will work with vectors, so you can solve

systems of difference equations as well as single ones.

2.4

How to document a program

One of the main themes of our study will be the preparation of programs that not only

work, but also are easily readable and useable by other people. The act of communication

that must take place before a program can be used by persons other than its author is a

difficult one to carry out, and we will return several times to the principles that serve as

guides to the preparation of readable software.

In this section we discuss further the all-important question of program documentation,

already touched upon in section 2.2. Some very nontrivial skills are called for in the creation

of good user-oriented program descriptions. One of these is the ability to enter the head of

another person, the user, and to relate to the program that you have just written through

the user’s eyes.

It’s hard to enter someone else’s head. One of the skills that make one person a better

teacher than another person is of the same kind: the ability to see the subject matter that

is being taught through the eyes of another person, the student. If one can do that, or

even make a good try at it, then obviously one will be able much better to deal with the

questions that are really concerning the audience. Relatively few actually do this to any

great extent not, I think, because it’s an ability that one either has or doesn’t have, but

because few efforts are made to train this skill and to develop it.

We’ll try to make our little contribution here.

2.4 How to document a program

35

(A) What does it do?

The first task should be to describe the precise purpose of the program. Put yourself in

the place of a potential user who is looking at a particular numerical instance of the problem

that needs solving. That user is now thumbing through a book full of program descriptions

in the library of a computer center 10,000 miles away in order to find a program that will

do the problem in question. Your description must answer that question.

Let’s now assume that your program has been written in the form of a subroutine or

procedure, rather than as a main program. Then the list of global, or communicating

variables is plainly in view, in the opening statement of the subroutine.

As we noted in section 2.2, you should state the purpose of your program using the

global variables of the subroutine in the same sentence. For one example of the difference

that makes, see section 2.2. For another, a linear equation solver might be described by

saying

“This program solves a system of simultaneous equations. To use it, put the

right-hand sides into the vector B, put the coefficients into the matrix A and call

the routine. The answer will be returned in X.”

We are, however, urging the reader to do it this way:

“This program solves the equations AX=B, where A is an N-by-N matrix and B

is an N-vector.”

Observe that the second description is shorter, only about half as long, and yet more

informative. We have found out not only what the program does, but how that function

relates to the global variables of the subroutine. This was done by using a judicious sprinkling of symbols in the documentation, along with the words. Don’t use only symbols, or

only words, but weave them together for maximum information.

Notice also that the ambiguous term “right-hand side” that appeared in the first form

has been done away with in the second form. The phrase was ambiguous because exactly

what ends up on the right-hand side and what on the left is an accident of how we happen

to write the equations, and your audience may not do it the same way you do.

(B) How is it done?

This is usually the easy part of program documentation because it is not the purpose of

this documentation to give a course in mathematics or algorithms or anything else. Hence

most of the time a reference to the literature is enough, or perhaps if the method is a

standard one, just give its name. Often though, variations on the standard method have

been chosen, and the user must be informed about those:

“. . . is solved by Gaussian elimination, using complete positioning for size. . . ”

“. . . the input array A is sorted by the Quicksort method (see D.E. Knuth, The

Art of Computer Programming, volume 3). . . ”

36

The Numerical Solution of Differential Equations

“. . . the eigenvalues and vectors are found by the Jacobi method, using Corbat?o’s method of avoiding the search for the largest off-diagonal element (see,

for instance, the description in D.R. Wilson, A First Course in Mathematical

Software).”

“. . . is found by the Simplex method, except that Charnes’ selection rule (see

F.A. Ficken, The Simplex Method. . . ) is not programmed, and so. . . ”

(C) Describe the global variables

Now it gets hard again. The global variables are the ones through which the subroutine

communicates with the user. Generally speaking, the user doesn’t care about variables

that are entirely local to your subroutine, but is vitally concerned with the communicating

variables.

First the user has to know exactly ho each of the global variables is related to the

problem that is being solved. This calls for a brief verbal description of the variable, and

what it has to do with the functioning of the program.

“A[i] is the ith element of the input list that is to be sorted, i=1..N”

“WHY is set by the subroutine to TRUE unles the return is because of overflow,

and then it will be set to FALSE.”

“B[i,j] is the coefficient of X[j] in the ith one of the input equations BX=C.”

“option is set by the calling program on input. Set it to 0 if the output is to

be rounded to the nearest integer, else set it to m if the output is to be rounded

to m decimal places (m ? 12).”

It is extremely important that each and every global variable of the subroutine should

get such a description. Just march through the parentheses in the subroutine or procedure

heading, and describe each variable in turn.

Next, the user will need more information about each of the global variables than just its

description as above. Also required is the “type” of the variable. Some computer languages

force each program to declare the types of their variables right in the opening statement.

Others declare types by observing various default rules with exceptions stated. In any case,

a little redundancy never hurts, and the program documentation should declare the type of

each and every global variable.

It’s easy to declare along with the types of the variables, their dimensions if they are

array variables. For instance we may have a

solver:=proc(A,X,n,ndim,b);

2.4 How to document a program

37

in which the communicating variables have the following types:

A

X

n

ndim

b

ndim-by-ndim array of floating point numbers

vector of floating point numbers of length n

integer

integer

vector of floating point numbers of length n

The best way to announce all of these types and dimensions of global variables to the

user is simply to list them, as above, in a table.

Now surely we’ve finished taking the pulse, blood pressure, etc. of the global variables,

haven’t we? Well, no, we haven’t. There’s still more vital data that a user will need to

know about these variables. There isn’t any standard name like “type” to apply to this

information, so we’ll call it the “role” of the variable.

First, for some of the global variables of the subroutine, it may be true that their values

at the time the subroutine is called are quite irrelevant to the operation of the subroutine.

This would be the case for the output variables, and in certain other situations. For some

other variables, the values at input time would be crucial. The user needs to know which are

which. Just for one example, if the value at input time is irrelevant, then the user can feel

free to use the same storage for other temporary purposes between calls to the subroutine.

Second, it may happen that certain variables are returned by the subroutine with their

values unchanged. This is particularly true for “implicitly passed” global variables, i.e.,

variables whose values are used by the subroutine but which do not appear explicitly in the

argument list. In such cases, the user may be delighted to hear the good news. In other

cases, the action of a subroutine may change an input variable, so if the user needs to use

those quantities again it will be necessary to save them somewhere else before calling the

subroutine. In either case, the user needs to be informed.

Third, it may be that the computation of the value of a certain variable is one of the

main purposes of the subroutine. Such variables are the outputs of the program, and the

user needs to know which these are (whether they are explicit in heading or the return

statement, or are “implicit”).

Although some high-level computer languages require type declarations immediately

in the opening instruction of a subroutine, none require the descriptions of the roles of

the variables (well, Pascal requires the VAR declaration, and Maple separates the input

variables from the output ones, but both languages allow implicit passing and changing

of global variables). These are, however, important for the user to know, so let’s invent

a shorthand for describing them in the documentation of the programs that occur in this

book.

First, if the value at input time is important, let’s say that the role of the variable is I,

otherwise it is I’.

Second, if the value of the variable is changed by the action of the subroutine, we’ll say

that its role is C, else C’.

Finally, if the computation of this variable is one of the main purposes of the subroutine,

38

The Numerical Solution of Differential Equations

it’s role is O (as in output), else O’.

In the description of each communicating variable, all three of these should be specified,.

Thus, a variable X might have role IC’O’, or a variable why might be of role I’CO, etc.

To sum up, the essential features of program documentation are a description of that

the program does, phrased in terms of the global variables, a statement of how it gets the

job done, and a list of all of the global variables, showing for each one its name, type,

dimension (or structure) if any, its role, and a brief verbal description.

Refer back to the short program in section 2.2, that searches for the largest element in

a row of a matrix. Here is the table of information about its global variables:

Name

A

i

jwin

Type

floating point matrix

integer

integer

Role

IC’O’

IC’O’

I’CO

Description

The input matrix

Which row to search

Column containing largest element

Exercises 2.4

Write programs that perform each of the jobs stated below. In each case, after testing

the program, document it with comments. Give a complete table of information about the

global variables in each case.

(a) Find and print all of the prime numbers between M and N .

(b) Find the elements of largest and smallest absolute values in a given linear array

(vector), and their positions in the array.

(c) Sort the elements of a given linear array into ascending order of size.

(d) Deal out four bridge hands (13 cards each from a standard 52-card deck – This one

is not so easy!).

(e) Solve a quadratic equation (any quadratic equation!).

2.5

The midpoint and trapezoidal rules

Euler’s formula is doubtless the simplest numerical integration procedure for differential

equations, but the accuracy that can be obtained with it is insufficient for most applications.

In this section and those that follow, we want to introduce a while family of methods for

the solution of differential equations, called the linear multistep methods, in which the user

can choose the degree of precision that will suffice for the job, and then select a member of

the family that will achieve it.

Before describing the family in all of its generality, we will produce two more of its

members, which illustrate the different sorts of creatures that inhabit the family in question.

2.5 The midpoint and trapezoidal rules

39

Recall that we derived Euler’s method by chopping off the Taylor series expansion of the

solution after the linear term. To get a more accurate method we could, of course, keep the

quadratic term, too. However, that term involves a second derivative, and we want to avoid

the calculation of higher derivatives because our differential equations will always be written

as first-order systems, so that only the first derivative will be conveniently computable.

We can have greater accuracy without having to calculate higher derivatives if we’re

willing to allow our numerical integration procedure to involve values of the unknown function and its derivative at more than one point. In other words, in Euler’s method, the next

value of the unknown function, at x + h, is gotten from the values of y and y 0 at just one

backwards point x. In the more accurate formulas that we will discuss next, the new value

of y depends in y and y 0 at more than one point, for instance, at x and x ? h, or at several

points.

As a primitive example of this kind, we will now discuss the midpoint rule. We begin

once again with the Taylor expansion of the unknown function y(x) about the point xn :

y(xn + h) = y(xn ) + hy 0 (xn ) + h2

y 00 (xn )

y 000 (xn )

+ h3

+ ···.

2

6

(2.5.1)

Now we rewrite equation (2.5.1) with h replaced by ?h to get

y(xn ? h) = y(xn ) ? hy 0 (xn ) + h2

y 00 (xn )

y 000 (xn )

? h3

+ ···

2

6

(2.5.2)

and then subtract these equations, obtaining

y(xn + h) ? y(xn ? h) = 2hy 0 (xn ) + 2h3

y 000 (xn )

+ ···.

6

(2.5.3)

Now, just as we did in the derivation of Euler’s method, we will truncate the right side

of (2.5.3) after the first term, ignoring the terms that involve h3 , h5 , etc. Further, let’s

use yn to denote the computed approximate value of y(xn ) (and yn+1 for the approximate

y(xn+1 ), etc.). Then we have

yn+1 ? yn?1 = 2hyn0 .

(2.5.4)

If, as usual, we are solving the differential equation y 0 = f (x, y), then finally (2.5.4) takes

the form

yn+1 = yn?1 + 2hf (xn , yn )

(2.5.5)

and this is the midpoint rule. The name arises from the fact that the first derivative yn0 is

being approximated by the slope of the chord that joins the two points (xn?1 , yn?1 ) and

(xn+1 , yn+1 ), instead of the chord joining (xn , yn ) and (xn+1 , yn+1 ) as in Euler’s method.

At first sight it seems that (2.5.5) can be used just like Euler’s method, because it is a

recurrence formula in which we compute the next value yn+1 from the two previous values

yn and yn?1 . Indeed the rules are quite similar, except for the fact that we can’t get started

with the midpoint rule until we know two consecutive values y0 , y1 of the unknown function

at two consecutive points x0 , x1 . Normally a differential equation is given together with

just one value of the unknown function, so if we are to use the midpoint rule we’ll need to

manufacture one more value of y(x) by some other means.

40

The Numerical Solution of Differential Equations

This kind of situation will come up again and again as we look at more accurate methods,

because to obtain greater precision without computing higher derivatives we will get the

next approximate value of y from a recurrence formula that may involve not just one or two,

but several of its predecessors. To get such a formula started we will have to find several

starting values in addition to the one that is given in the statement of the initial-value

problem.

To get back to the midpoint rule, we can get it started most easily by calculating y1 ,

the approximation to y(x0 + h), from Euler’s method, and then switching to the midpoint

rule to carry out the rest of the calculation.

Let’s do this, for example with the same differential equation (2.1.7) that we used to

illustrate Euler’s rule, so we can compare the two methods. The problem consists of the

equation y 0 = 0.5y and the initial value y(0) = 1. We’ll use the same step size h = 0.05 as

before.

Now to start the midpoint rule we need two consecutive values of y, in this case at x = 0

and x = 0.05. At 0.05 we use the value that Euler’s method gives us, namely y1 = 1.025

(see Table 1). It’s easy to continue the calculation now from (2.5.5).

For instance

y2 = y0 + 2h(0.5y1 )

= 1 + 0.1(0.5 ? 1.025)

(2.5.6)

= 1.05125

and

y3 = y1 + 2h(0.5y2 )

= 1.025 + 0.1(0.5 ? 1.05125)

(2.5.7)

= 1.0775625 .

In the table below we show for each x the value computed from the midpoint rule, from

Euler’s method, and from the exact solution y(x) = ex/2 . The superior accuracy of the

midpoint rule is apparent.

x

0.00

0.05

0.10

0.15

0.20

0.25

..

.

Midpoint(x)

1.00000

1.02500

1.05125

1.07756

1.10513

1.13282

..

.

Euler(x)

1.00000

1.02500

1.05063

1.07689

1.10381

1.13141

..

.

Exact(x)

1.00000

1.02532

1.05127

1.07788

1.10517

1.13315

..

.

1.00

2.00

3.00

..

.

1.64847

2.71763

4.48032

..

.

1.63862

2.68506

4.39979

..

.

1.64872

2.71828

4.48169

..

.

5.00

10.00

12.17743

148.31274

11.81372

139.56389

12.18249

148.41316

2.5 The midpoint and trapezoidal rules

41

y = f (x)

a

b

Figure 2.1: The trapezoidal rule

table 2

Next, we introduce a third method of numerical integration, the trapezoidal rule. The

best way to obtain it is to convert the differential equation that we’re trying to solve into

an integral equation, and then use the trapezoidal approximation for the integral.

We begin with the differential equation y 0 = f (x, y(x)), and we integrate both sides

from x to x + h, getting

Z

x+h

y(x + h) = y(x) +

f (t, y(t)) dt.

(2.5.8)

x

Now if we approximate the right-hand side in any way by a weighted sum of values of

the integrand at various points we will have found an approximate method for solving our

differential equation.

The trapezoidal rule states that for an approximate value of an integral

Z

b

f (t) dt

(2.5.9)

a

we can use, instead of the area under the curve between x = a and x = b, the area of the

trapezoid whose sides are the x axis, the lines x = a and x = b, and the line through the

points (a, f (a)) and (b, f (b)), as shown in Figure 2.1. That area is 12 (f (a) + f (b))(b ? a).

If we apply the trapezoidal rule to the integral that appears in (2.5.8), we obtain

y(xn + h) ? y(xn ) +

h

(f (xn , y(xn )) + f (xn + h, y(xn + h)))

2

(2.5.10)

in which we have used the “?” sign rather than the “=” because the right hand side is not

exactly equal to the integral that really belongs there, but is only approximately so.

If we use our usual abbreviation yn for the computed approximate value of y(xn ), then

(2.5.10) becomes

h

yn+1 = yn + (f (xn , yn ) + f (xn+1 , yn+1 )).

(2.5.11)

2

42

The Numerical Solution of Differential Equations

This is the trapezoidal rule in the form that is useful for differential equations.

At first sight, (2.5.11) looks like a recurrence formula from which the next approximate

value, yn+1 , of the unknown function, can immediately be computed from the previous

value, yn . However, this is not the case.

Upon closer examination one observes that the next value yn+1 appears not only on the

left-hand side, but also on the right (it’s hiding in the second f on the right side).

In order to find the value yn+1 it appears that we need to carry out an iterative process.

First we would guess yn+1 (guessing yn+1 to be equal to yn wouldn’t be all that bad, but

we can do better). If we use this guess value on the right side of (2.5.11) then we would

be able to calculate the entire right-hand side, and then we could use that value as a new

“improved” value of yn+1 .

Now if the new value agrees with the old sufficiently well the iteration would halt, and

we would have found the desired value of yn+1 . Otherwise we would use the improved value

on the right side just as we previously used the first guess. Then we would have a “more

improved” guess, etc.

Fortunately, in actual use, it turns out that one does not actually have to iterate to

convergence. If a good enough guess is available for the unknown value, then just one

refinement by a single application of the trapezoidal formula is sufficient. This is not the

case if a high quality guess is unavailable. We will discuss this point in more detail in

section 2.9. The pair of formulas, one of which supplies a very good guess to the next

value of y, and the other of which refines it to a better guess, is called a predictor-corrector

pair, and such pairs form the basis of many of the highly accurate schemes that are used in

practice.

As a numerical example, take the differential equation

y 0 = 2xey + 1

(2.5.12)

with the initial value y(0) = 1. If we use h = 0.05, then our first task is to calculate y1 , the

approximate value of y(0.05). The trapezoidal rule asserts that

y1 = 1 + 0.025(2 + 0.1ey1 )

(2.5.13)

and sure enough, the unknown number y1 appears on both sides.

Let’s guess y1 = 1. Since this is not a very inspired guess, we will have to iterate the

trapezoidal rule to convergence. Hence, we use this guess on the right side of (2.5.13),

compute the right side, and obtain y1 = 1.056796. If we use this new guess the same way,

the result is y1 = 1.057193. Then we get 1.057196, and since this is in sufficiently close

agreement with the previous result, we declare that the iteration has converged. Then we

take y1 = 1.057196 for the computed value of the unknown function at x = 0.05, and we go

next to x = 0.1 to repeat the same sort of thing to get y2 , the computed approximation to

y(0.1).

In Table 3 we show the results of using the trapezoidal rule (where we have iterated

until two successive guesses are within 10?4 ) on our test equation y 0 = 0.5y, y(0) = 1 as

the column Trap(x). For comparison, we show Midpoint(x) and Exact(x).

2.6 Comparison of the methods

43

x

0.00

0.05

0.10

0.15

0.20

0.25

..

.

Trap(x)

1.00000

1.02532

1.05127

1.07789

1.10518

1.13316

..

.

Midpoint(x)

1.00000

1.02500

1.05125

1.07756

1.10513

1.13282

..

.

Exact(x)

1.00000

1.02532

1.05127

1.07788

1.10517

1.13315

..

.

1.00

2.00

3.00

..

.

1.64876

2.71842

4.48203

..

.

1.64847

2.71763

4.48032

..

.

1.64872

2.71828

4.48169

..

.

5.00

10.00

12.18402

148.45089

12.17743

148.31274

12.18249

148.41316

table 3

2.6

Comparison of the methods

We are now in possession of three methods for the numerical solution of differential equations. They are Euler’s method

yn+1 = yn + hyn0 ,

(2.6.1)

the trapezoidal rule

yn+1 = yn +

h 0

0

),

(y + yn+1

2 n

(2.6.2)

and the midpoint rule

yn+1 = yn?1 + 2hyn0 .

(2.6.3)

In order to compare the performance of the three techniques it will be helpful to have

a standard differential equation on which to test them. The most natural candidate for

such an equation is y 0 = Ay, where A is constant. The reasons for this choice are first

that the equation is easy to solve exactly, second that the difference approximations are

also relatively easy to solve exactly, so comparison is readily done, third that by varying the

sign of A we can study behavior of either growing or shrinking (stable or unstable) solutions,

and finally that many problems in nature have solutions that are indeed exponential, at

least over the short term, so this is an important class of differential equations.

We will, however, write the test equation in a slightly different form for expository

reasons, namely as

y

y0 = ? ;

y(0) = 1 ;

L>0

(2.6.4)

L

where L is a constant. The most interesting and revealing case is where the true solution

is a decaying exponential, so we will assume that L > 0. Further, we will assume that

y(0) = 1 is the given initial value.

44

The Numerical Solution of Differential Equations

The exact solution is of course

Exact(x) = e?x/L .

(2.6.5)

Notice that if x increases by L, the solution changes by a factor of e. Hence L, called the

relaxation length of the problem, can be conveniently visualized as the distance over which

the solution falls by a factor of e.

Now we would like to know how well each of the methods (2.6.1)–(2.6.3) handles the

problem (2.6.4).

Suppose first that we ask Euler’s method to solve the problem. If we substitute y 0 =

f (x, y) = ?y/L into (2.6.1), we get

yn+1

yn

= yn + h ? ?

L

h

= 1?

yn .

L

(2.6.6)

Before we solve this recurrence, let’s comment on the ratio h/L that appears in it. Now

L is the distance over which the solution changes by a factor of e, and h is the step size

that we are going to use in the numerical integration. Instinctively, one feels that if the

solution is changing rapidly in a certain region, then h will have to be kept small there if

good accuracy is to be retained, while if the solution changes only slowly, then h can be

larger without sacrificing too much accuracy. The ratio h/L measures the step size of the

integration in relation to the distance over which the solution changes appreciably. Hence,

h/L is exactly the thing that one feels should be kept small for a successful numerical

solution.

Since h/L occurs frequently below, we will denote it with the symbol ? .

Now the solution of the recurrence equation (2.6.6), with the starting value y0 = 1, is

obviously

yn = (1 ? ? )n

n = 0, 1, 2, . . . .

(2.6.7)

Next we study the trapezoidal approximation to the same equation (2.6.4). We substitute y 0 = f (x, y) = ?y/L in (2.6.2) and get

yn+1

h

yn yn+1

= yn +

? ?

.

2

L

L

(2.6.8)

The unknown yn+1 appears, as usual with this method, on both sides of the equation.

However, for the particularly simple equation that we are now studying, there is no difficulty

in solving (2.6.8) for yn+1 (without any need for an iterative process) and obtaining

yn+1 =

1?

1+

?

2

?

2

yn .

(2.6.9)

Together with the initial value y0 = 1, this implies that

yn =

1?

1+

?

2

?

2

!n

n = 0, 1, 2, . . . .

(2.6.10)

2.6 Comparison of the methods

45

Before we deal with the midpoint rule, let’s pause to examine the two methods whose

solutions we have just found. Note that for a given value of h, all three of (a) the exact solution, (b) Euler’s solution and (c) the trapezoidal solution are of the form yn = (constant)n ,

in which the three values of “constant” are

(a) e??

(b) 1 ? ?

(c)

1? ?2

1+ ?2

(2.6.11)

.

It follows that to compare the two approximate methods with the “truth,” all we have

to do is see how close the constants (b) and (c) above are to the true constant (a). If we

remember that ? is being thought of as small compared to 1, then we have the power series

expansion of e??

?2 ?3

e?? = 1 ? ? +

?

+ ···

(2.6.12)

2

6

to compare with 1 ? ? and with the power series expansion of

1?

1+

?

2

?

2

=1?? +

?2 ?3

?

+ ···.

2

4

(2.6.13)

The comparison is now clear. Both the Euler and the trapezoidal methods yield approximate solutions of the form (constant)n , where “constant” is near e?? . The trapezoidal

rule does a better job of being near e?? because its constant agrees with the power series

expansion of e?? through the quadratic term, whereas that of the Euler method agrees only

up to the linear term.

Finally we study the nature of the approximation that is provided by the midpoint

rule. We will find that a new and important phenomenon rears its head in this case. The

analysis begins just as it did in the previous two cases: We substitute the right-hand side

f (x, y) = ?y/L for y 0 in (2.6.3) to get

yn+1 = yn?1 + 2h ? ?

yn

.

L

(2.6.14)

One important feature is already apparent. Instead of facing a first-order difference

equation as we did in (2.6.6) for Euler’s method and in (2.6.9) for the trapezoidal rule, we

have now to contend with a second-order difference equation.

Since the equation is linear with constant coefficients, we know to try a solution of the

form yn = r n . This leads to the quadratic equation

r 2 + 2? r ? 1 = 0.

(2.6.15)

Evidently the discriminant of this equation is positive, so its roots are distinct. If we denote

these two roots by r+ (? ) and r? (? ), then the general solution of the difference equation

(2.6.14) is

yn = c (r+ (? ))n + d (r? (? ))n ,

(2.6.16)

46

The Numerical Solution of Differential Equations

where c and d are constants whose values are determined by the initial data y0 and y1 .

The Euler and trapezoidal approximations were each of the form (constant)n . This one

is a sum of two terms of that kind. We will see that r+ (? ) is a very good approximation

to e?? . The other term, (r? (? ))n is, so to speak, the price that we pay for getting such a

good approximation in r+ (? ). We hope that the other term will stay small relative to the

first term, so as not to disturb the closeness of the approximation. We will see, however,

that it need not be so obliging, and that in fact it might do whatever it can to spoil things.

The two roots of the quadratic equation are

?

r+ (? ) = ?? + ?1 + ? 2

r? (? ) = ?? ? 1 + ? 2 .

(2.6.17)

When ? = 0 the first of these is +1, so when ? is small r+ (? ) is near +1, and it is the root

that is trying to approximate the exact constant e?? as well as possible. In fact it does

pretty well, because the power series expansion of r+ (? ) is

r+ (? ) = 1 ? ? +

?2 ?4

?

+ ···

2

8

(2.6.18)

so it agrees with e?? through the quadratic terms.

What about r? (? )? Its Taylor series is

r? (? ) = ?1 ? ? ?

?2

+ ···.

2

(2.6.19)

The bad news is now before us: When ? is a small positive number, the root r? (? ) is

larger than 1 in absolute value. This means that the stability criterion of Theorem 1.6.1 is

violated, so we say that the midpoint rule is unstable.

In practical terms, we observe that r+ (? ) is close to e?? , so the first term on the right

of (2.6.16) is very close to the exact solution. The second term of (2.6.16), the so-called

parasitic solution, is small compared to the first term when n is small, because the constant

d will be small compared with c. However, as we move to the right, n increases, and the

second term will eventually dominate the first, because the first term is shrinking to zero

as n increases, because that’s what the exact solution does, while the second term increases

steadily in size. In fact, since r? (? ) is negative and larger than 1 in absolute value, the

second term will alternate in sign as n increases, and grow without bound in magnitude.

In Table 4 below we show the result of integrating the problem y 0 = ?y, y(0) = 1 with

each of the three methods that we have discussed, using a step size of h = 0.05. To get the

midpoint method started, we used the exact value of y(0.05) (i.e., we cheated), and in the

trapezoidal rule we iterated to convergence with = 10?4 . The instability of the midpoint

rule is quite apparent.

2.6 Comparison of the methods

x

0.0

1.0

2.0

3.0

4.0

5.0

10.0

14.55

15.8

Euler(x)

1.00000

0.35849

0.12851

0.04607

0.01652

0.00592

0.00004

3.3 ? 10?7

9.1 ? 10?8

47

Trap(x)

1.00000

0.36780

0.13527

0.04975

0.01830

0.00673

0.00005

4.8 ? 10?7

1.4 ? 10?7

Midpoint(x)

1.00000

0.36806

0.13552

0.05005

0.01888

0.00822

0.21688

-20.48

71.45

Exact(x)

1.00000

0.36788

0.13534

0.04979

0.01832

0.00674

0.00005

4.8 ? 10?7

1.4 ? 10?7

table 4

In addition to the above discussion of accuracy, we summarize here there additional

properties of integration methods as they relate to the examples that we have already

studied.

First, a numerical integration method might be iterative or noniterative. A method is

noniterative if it expresses the next value of the unknown function quite explicitly in terms

of values of the function and its derivatives at preceding points. In an iterative method,

at each step of the solution process the next value of the unknown is defined implicitly by

an equation, which must be solved to obtain the next value of the unknown function. In

practice, we may either solve this equation completely by an iteration or