h„n+1

y(x)

u— en+1

n+1

yn

un+1

un

tn+1

tn

FIGURE 11.1. Geometrical interpretation of the local and global truncation er-

rors at node tn+1 for the forward Euler method

As a consequence,

|en+1 | ¤ h|„n+1 (h)| + |en | + h|f (tn , yn ) ’ f (tn , un )| ¤ h„ (h) + (1 + hL)|en |,

L being the Lipschitz constant of f . By recursion on n, we ¬nd

|en+1 | ¤ [1 + (1 + hL) + . . . + (1 + hL)n ] h„ (h)

eL(tn+1 ’t0 ) ’ 1

(1 + hL)n+1 ’ 1

„ (h) ¤

= „ (h).

L L

11.3 Analysis of One-Step Methods 479

The last inequality follows from noticing that 1 + hL ¤ ehL and (n + 1)h =

tn+1 ’ t0 .

On the other hand, if y ∈ C 2 (I), the LTE for the forward Euler method

is (see Section 10.10.1)

h

y (ξ), ξ ∈ (tn , tn+1 ),

„n+1 (h) =

2

and thus, „ (h) ¤ (M/2)h, where M = maxξ∈I |y (ξ)|. In conclusion,

eL(tn+1 ’t0 ) ’ 1 M

|en+1 | ¤ h, ∀n ≥ 0, (11.22)

L 2

from which it follows that the global error tends to zero with the same

order as the local truncation error.

If also the rounding errors are accounted for, we can assume that the

solution un+1 , actually computed by the forward Euler method at time

¯

tn+1 , is such that

u0 = y0 + ζ0 , un+1 = un + hf (tn , un ) + ζn+1 ,

¯ ¯ ¯ ¯ (11.23)

having denoted the rounding error by ζj , for j ≥ 0.

Problem (11.23) is an instance of (11.16), provided that we identify ζn+1

(h)

and un with hδn+1 and zn in (11.16), respectively. Combining Theorems

¯

11.1 and 11.2 we get, instead of (11.22), the following error estimate

1 M ζ

|yn+1 ’ un+1 | ¤ eL(tn+1 ’t0 ) |ζ0 | +

¯ h+ ,

L 2 h

where ζ = max1¤j¤n+1 |ζj |. The presence of rounding errors does not allow,

therefore, to conclude that as h ’ 0, the error goes to zero. Actually,

there exists an optimal (non null) value of h, hopt , for which the error is

minimized. For h < hopt , the rounding error dominates the truncation error

and the global error increases.

11.3.3 The Absolute Stability

The property of absolute stability is in some way specular to zero-stability,

as far as the roles played by h and I are concerned. Heuristically, we say that

a numerical method is absolutely stable if, for h ¬xed, un remains bounded

as tn ’ +∞. This property, thus, deals with the asymptotic behavior of

un , as opposed to a zero-stable method for which, for a ¬xed integration

interval, un remains bounded as h ’ 0.

For a precise de¬nition, consider the linear Cauchy problem (that from now

on, we shall refer to as the test problem)

y (t) = »y(t), t > 0,

(11.24)

y(0) = 1,

480 11. Numerical Solution of Ordinary Di¬erential Equations

with » ∈ C, whose solution is y(t) = e»t . Notice that lim |y(t)| = 0 if

t’+∞

Re(») < 0.

De¬nition 11.6 A numerical method for approximating (11.24) is abso-

lutely stable if

|un | ’’ 0 tn ’’ +∞.

as (11.25)

Let h be the discretization stepsize. The numerical solution un of (11.24)

obviously depends on h and ». The region of absolute stability of the nu-

merical method is the subset of the complex plane

A = {z = h» ∈ C : (11.25) is satis¬ed } . (11.26)

Thus, A is the set of the values of the product h» for which the numerical

method furnishes solutions that decay to zero as tn tends to in¬nity.

Let us check whether the one-step methods introduced previously are ab-

solutely stable.

1. Forward Euler method: applying (11.7) to problem (11.24) yields un+1 =

un + h»un for n ≥ 0, with u0 = 1. Proceeding recursively on n we get

n ≥ 0.

un = (1 + h»)n ,

Therefore, condition (11.25) is satis¬ed i¬ |1 + h»| < 1, that is, if h» lies

within the unit circle with center at (’1, 0) (see Figure 11.3). This amounts

to requiring that

2Re(»)

h» ∈ C’ and 0 < h < ’ (11.27)

|»|2

where

C’ = {z ∈ C : Re(z) < 0} .

Example 11.1 For the Cauchy problem y (x) = ’5y(x) for x > 0 and y(0) = 1,

condition (11.27) implies 0 < h < 2/5. Figure 11.2 (left) shows the behavior of

the computed solution for two values of h which do not ful¬ll this condition, while

on the right we show the solutions for two values of h that do. Notice that in this

•

second case the oscillations, if present, damp out as t grows.

2. Backward Euler method: proceeding as before, we get this time

1

n ≥ 0.

un = ,

(1 ’ h»)n

The absolute stability property (11.25) is satis¬ed for any value of h» that

does not belong to the unit circle of center (1, 0) (see Figure 11.3, right).

11.3 Analysis of One-Step Methods 481

3

1

0.8

2

0.6

0.4

1

0.2

0 0

’0.2

’1

’0.4

’0.6

’2

’0.8

’3 ’1

0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8

FIGURE 11.2. Left: computed solutions for h = 0.41 > 2/5 (dashed line) and

h = 2/5 (solid line). Notice how, in the limiting case h = 2/5, the oscillations

remain unmodi¬ed as t grows. Right: two solutions are reported for h = 0.39

(solid line) and h = 0.15 (dashed line)

Example 11.2 The numerical solution given by the backward Euler method in

the case of Example 11.1 does not exhibit any oscillation for any value of h. On

the other hand, the same method, if applied to the problem y (t) = 5y(t) for t > 0

and with y(0) = 1, computes a solution that decays anyway to zero as t ’ ∞ if

h > 2/5, despite the fact that the exact solution of the Cauchy problem tends to

•

in¬nity.

3. Trapezoidal (or Crank-Nicolson) method: we get

n

1 1

1 + »h / 1 ’ »h n ≥ 0,

un = ,

2 2

hence (11.25) is ful¬lled for any h» ∈ C’ .

4. Heun™s method: applying (11.10) to problem (11.24) and proceeding

by recursion on n, we obtain

n

(h»)2

n ≥ 0.

un = 1 + h» + ,

2

As shown in Figure 11.3 the region of absolute stability of Heun™s method

is larger than the corresponding one of Euler™s method. However, its re-

striction to the real axis is the same.

We say that a method is A-stable if A © C’ = C’ , i.e., if for Re(») < 0,

condition (11.25) is satis¬ed for all values of h.

The backward Euler and Crank-Nicolson methods are A-stable, while

the forward Euler and Heun methods are conditionally stable.

Remark 11.2 Notice that the implicit one-step methods examined so far

are unconditionally absolutely stable, while explicit schemes are condition-

482 11. Numerical Solution of Ordinary Di¬erential Equations

Im

1.75

H

BE

FE

’1 Re

1

’1.75

FIGURE 11.3. Regions of absolute stability for the forward (FE) and backward

Euler (BE) methods and for Heun™s method (H). Notice that the region of abso-

lute stability of the BE method lies outside the unit circle of center (1, 0) (shaded

area)

ally absolutely stable. This is, however, not a general rule: in fact, there ex-

ist implicit unstable or only conditionally stable schemes. On the contrary,

there are no explicit unconditionally absolutely stable schemes [Wid67].

11.4 Di¬erence Equations

For any integer k ≥ 1, an equation of the form

un+k + ±k’1 un+k’1 + . . . + ±0 un = •n+k , n = 0, 1, . . . (11.28)

is called a linear di¬erence equation of order k. The coe¬cients ±0 = 0,

±1 , . . . , ±k’1 may or may not depend on n. If, for any n, the right side

•n+k is equal to zero, the equation is said homogeneous, while if the ±j s

are independent of n it is called linear di¬erence equation with constant

coe¬cients.

Di¬erence equations arise for instance in the discretization of ordinary dif-

ferential equations. Regarding this, we notice that all the numerical meth-

ods examined so far end up with equations like (11.28). More generally,

equations like (11.28) are encountered when quantities are de¬ned through

linear recursive relations. Another relevant application is concerned with

the discretization of boundary value problems (see Chapter 12). For fur-

ther details on the subject, we refer to Chapters 2 and 5 of [BO78] and to

Chapter 6 of [Gau97].

11.4 Di¬erence Equations 483

Any sequence {un , n = 0, 1, . . . } of values that satisfy (11.28) is called

a solution to the equation (11.28). Given k initial values u0 , . . . , uk’1 , it

is always possible to construct a solution of (11.28) by computing (sequen-

tially)

un+k = [•n+k ’ (±k’1 un+k’1 + . . . + ±0 un )], n = 0, 1, . . .

However, our interest is to ¬nd an expression of the solution un+k which

depends only on the coe¬cients and on the initial values.

We start by considering the homogeneous case with constant coe¬cients,

un+k + ±k’1 un+k’1 + . . . + ±0 un = 0, n = 0, 1, . . . (11.29)

and associate with (11.29) the characteristic polynomial Π ∈ Pk de¬ned as

Π(r) = rk + ±k’1 rk’1 + . . . + ±1 r + ±0 . (11.30)

Denoting its roots by rj , j = 0, . . . , k ’ 1, any sequence of the form

for j = 0, . . . , k ’ 1

n

rj , n = 0, 1, . . . , (11.31)

is a solution of (11.29), since

n+k n+k’1 n

rj + ±k’1 rj + . . . + ±0 rj

k’1

n k n

= rj rj + ±k’1 rj + . . . + ±0 = rj Π(rj ) = 0.

We say that the k sequences de¬ned in (11.31) are the fundamental solutions

of the homogeneous equation (11.29). Any sequence of the form

n n n

un = γ0 r0 + γ1 r1 + . . . + γk’1 rk’1 , n = 0, 1, . . . (11.32)

is still a solution to (11.29), since it is a linear equation.

The coe¬cients γ0 , . . . , γk’1 can be determined by imposing the k initial

conditions u0 , . . . , uk’1 . Moreover, it can be proved that if all the roots

of Π are simple, then all the solutions of (11.29) can be cast in the form

(11.32).

This last statement no longer holds if there are roots of Π with multi-

plicity greater than 1. If, for a certain j, the root rj has multiplicity m ≥ 2,

in order to obtain a system of fundamental solutions that generate all the

solutions of (11.29), it su¬ces to replace the corresponding fundamental

n

solution rj , n = 0, 1, . . . with the m sequences

n n

nm’1 rj , n = 0, 1, . . . .

n

rj , n = 0, 1, . . . , nrj , n = 0, 1, . . . , . . . ,

More generally, assuming that r0 , . . . , rk are distinct roots of Π, with mul-

tiplicities equal to m0 , . . . , mk , respectively, we can write the solution of

(11.29) as

mj ’1

k

γsj ns n

un = rj , n = 0, 1, . . . . (11.33)

s=0

j=0

484 11. Numerical Solution of Ordinary Di¬erential Equations

Notice that even in presence of complex conjugate roots one can still obtain

a real solution (see Exercise 3).

Example 11.3 For the di¬erence equation un+2 ’un = 0, we have Π(r) = r2 ’1,

then r0 = ’1 and r1 = 1, therefore the solution is given by un = γ00 (’1)n + γ01 .

In particular, enforcing the initial conditions u0 and u1 gives γ00 = (u0 ’ u1 )/2,

•

γ01 = (u0 + u1 )/2.

Example 11.4 For the di¬erence equation un+3 ’ 2un+2 ’ 7un+1 ’ 4un = 0,

Π(r) = r3 ’ 2r2 ’ 7r ’ 4. Its roots are r0 = ’1 (with multiplicity 2), r1 = 4 and

the solution is un = (γ00 + nγ10 )(’1)n + γ01 4n . Enforcing the initial conditions

we can compute the unknown coe¬cients as the solution of the following linear

system

±

γ00 + γ01 = u0 ,

’γ00 ’ γ10 + 4γ01 = u1 ,

γ00 + 2γ10 + 16γ01 = u2

that yields γ00 = (24u0 ’ 2u1 ’ u2 )/25, γ10 = (u2 ’ 3u1 ’ 4u0 )/5 and γ01 =

•

(2u1 + u0 + u2 )/25.

The expression (11.33) is of little practical use since it does not outline

the dependence of un on the k initial conditions. A more convenient rep-

(n)

resentation is obtained by introducing a new set ψj , n = 0, 1, . . . of

fundamental solutions that satisfy

(i)

i, j = 0, 1, . . . , k ’ 1. (11.34)

ψj = δij ,

Then, the solution of (11.29) subject to the initial conditions u0 , . . . , uk’1

is given by

k’1

(n)

un = uj ψj , n = 0, 1, . . . . (11.35)

j=0

(n)

The new fundamental solutions ψj , n = 0, 1, . . . can be represented in

terms of those in (11.31) as follows

k’1

(n)

for j = 0, . . . , k ’ 1, n = 0, 1, . . .

n

ψj = βj,m rm (11.36)

m=0

By requiring (11.34), we obtain the k linear systems

k’1

i, j = 0, . . . , k ’ 1,

i

βj,m rm = δij ,

m=0

whose matrix form is

j = 0, . . . , k ’ 1.

Rbj = ej , (11.37)

11.4 Di¬erence Equations 485

Here ej denotes the unit vector of Rk , R = (rim ) = (rm ) and bj = i

(βj,0 , . . . , βj,k’1 )T . If all rj s are simple roots of Π, the matrix R is nonsin-

gular (see Exercise 5).

The general case where Π has k + 1 distinct roots r0 , . . . , rk with multi-

plicities m0 , . . . , mk respectively, can be dealt with by replacing in (11.36)

rj , n = 0, 1, . . . with rj ns , n = 0, 1, . . . , where j = 0, . . . , k and

n n

s = 0, . . . , mj ’ 1.

Example 11.5 We consider again the di¬erence equation of Example 11.4. Here

we have {r0 , nr0 , r1 , n = 0, 1, . . . } so that the matrix R becomes

n n n

®0 ®

0

r0 0 r2 1 01

R = ° r0 r0 r2 » = ° ’1 ’1 4 » .

1 1 1

2 2 2

r0 2r0 r2 1 2 16

Solving the three systems (11.37) yields

24 4 1n

(n)

(’1)n ’ n(’1)n +

ψ0 = 4,

25 5 25

2 3 2n

(n)

= ’ (’1)n ’ n(’1)n +

ψ1 4,

25 5 25

1 1 1n

(n)

= ’ (’1)n + n(’1)n +

ψ2 4,

25 5 25

(n)

2

from which it can be checked that the solution un = uj ψ j coincides with

j=0

•

the one already found in Example 11.4.

Now we return to the case of nonconstant coe¬cients and consider the

following homogeneous equation

k