0 1 0 0 0 0 0 1

4

-1 2

1 0 0 0 0

3 3 3

18 9 2 6

2 - 11 0 0 0

11 11 11

48

- 36 16 3 12

3 - 25 0 0

25 25 25 25

300 300 200 75 12 60

4 - 137 - 137 0

137 137 137 137

360

- 450 400

- 225 72 10 60

5 - 147

147 147 147 147 147 137

TABLE 11.2. Coe¬cients of zero-stable BDF methods for p = 0, 1, . . . , 5

11.6.1 Consistency

The property of consistency of a multistep method introduced in De¬nition

11.9 can be veri¬ed by checking that the coe¬cients satisfy certain algebraic

equations, as stated in the following theorem.

Theorem 11.3 The multistep method (11.45) is consistent i¬ the following

algebraic relations among the coe¬cients are satis¬ed

p p p

aj = 1, ’ jaj + bj = 1. (11.52)

j=0 j=0 j=’1

Moreover, if y ∈ C q+1 (I) for some q ≥ 1, where y is the solution of the

Cauchy problem (11.1), then the method is of order q i¬ (11.52) holds and

the following additional conditions are satis¬ed

p p

i

(’j)i’1 bj = 1, i = 2, . . . , q.

(’j) aj + i

j=0 j=’1

Proof. Expanding y and f in a Taylor series yields, for any n ≥ p

yn’j = yn ’ jhyn + O(h2 ), fn’j = fn + O(h). (11.53)

Plugging these values back into the multistep scheme and neglecting the terms

in h of order higher than 1 gives

p p

yn+1 ’ aj yn’j ’ h bj fn’j

j=0 j=’1

p p p p p

= yn+1 ’ jaj yn ’ h bj fn ’ O(h ) aj ’

2

a j yn + h bj

j=0 j=0 j=’1 j=0 j=’1

p p p p p

= yn+1 ’ aj yn ’ hyn ’ ’ O(h2 ) aj ’

jaj + bj bj

j=0 j=0 j=’1 j=0 j=’1

where we have replaced yn by fn . From the de¬nition (11.47) we then obtain

p p p p p

h„n+1 (h) = yn+1 ’ aj yn ’ hyn ’ ’ O(h ) aj ’

2

jaj + bj bj

j=0 j=0 j’1 j=0 j=’1

494 11. Numerical Solution of Ordinary Di¬erential Equations

hence the local truncation error is

p

yn+1 ’ yn yn

1’

„n+1 (h) = + aj

h h j=0

p p p p

jaj ’ ’ O(h) aj ’

+yn bj bj .

j=0 j=’1 j=0 j=’1

Since, for any n, (yn+1 ’ yn )/h ’ yn , as h ’ 0, it follows that „n+1 (h) tends to

0 as h goes to 0 i¬ the algebraic conditions (11.52) are satis¬ed. The rest of the

proof can be carried out in a similar manner, accounting for terms of progressively

3

higher order in the expansions (11.53).

11.6.2 The Root Conditions

Let us employ the multistep method (11.45) to approximately solve the

model problem (11.24). The numerical solution satis¬es the linear di¬erence

equation

p p

un+1 = aj un’j + h» bj un’j , (11.54)

j=0 j=’1

which ¬ts the form (11.29). We can therefore apply the theory devel-

oped in Section 11.4 and look for fundamental solutions of the form uk =

[ri (h»)]k , k = 0, 1, . . . , where ri (h»), for i = 0, . . . , p, are the roots of the

polynomial Π ∈ Pp+1

Π(r) = ρ(r) ’ h»σ(r). (11.55)

We have denoted by

p p

’ p’j

bj rp’j

p+1 p+1

ρ(r) = r aj r , σ(r) = b’1 r +

j=0 j=0

the ¬rst and second characteristic polynomials of the multistep method

(11.45), respectively. The polynomial Π(r) is the characteristic polynomial

associated with the di¬erence equation (11.54), and rj (h») are its charac-

teristic roots.

The roots of ρ are ri (0), i = 0, . . . , p, and will be abbreviated henceforth

by ri . From the ¬rst condition in (11.52) it follows that if a multistep

method is consistent then 1 is a root of ρ. We shall assume that such a root

(the consistency root) is labelled as r0 (0) = r0 and call the corresponding

root r0 (h») the principal root.

De¬nition 11.10 (Root condition) The multistep method (11.45) is said

to satisfy the root condition if all roots ri are contained within the unit

11.6 Analysis of Multistep Methods 495

circle centered at the origin of the complex plane, otherwise, if they fall on

its boundary, they must be simple roots of ρ. Equivalently,

|rj | ¤ 1, j = 0, . . . , p;

(11.56)

furthermore, for those j such that |rj | = 1, then ρ (rj ) = 0.

De¬nition 11.11 (Strong root condition) The multistep method (11.45)

is said to satisfy the strong root condition if it satis¬es the root condition

and r0 = 1 is the only root lying on the boundary of the unit circle. Equiv-

alently,

|rj | < 1 j = 1, . . . , p. (11.57)

De¬nition 11.12 (Absolute root condition) The multistep method

(11.45) satis¬es the absolute root condition if there exists h0 > 0 such that

|rj (h»)| < 1 ∀h ¤ h0 .

j = 0, . . . , p,

11.6.3 Stability and Convergence Analysis for Multistep

Methods

Let us now examine the relation between root conditions and the stability

of multistep methods. Generalizing the De¬nition 11.4, we can get the

following.

De¬nition 11.13 (Zero-stability of multistep methods) The multi-

step method (11.45) is zero-stable if

∃h0 > 0, ∃C > 0 : ∀h ∈ (0, h0 ], |zn ’ u(h) | ¤ Cµ, 0 ¤ n ¤ Nh ,

(h)

n

(11.58)

(h) (h)

where Nh = max {n : tn ¤ t0 + T } and zn and un are, respectively, the

solutions of problems

± p p

(h)

z (h) (h)

n+1 = aj zn’j + h bj f (tn’j , zn’j ) + hδn+1 ,

(11.59)

(h) j=0 j=’1

(h)

zk = wk + δk , k = 0, . . . , p

496 11. Numerical Solution of Ordinary Di¬erential Equations

± p p

(h)

u (h) (h)

n+1 = aj un’j + h bj f (tn’j , un’j ),

(11.60)

(h) j=0 j=’1

(h)

uk = wk , k = 0, . . . , p

(h) (h)

for p ¤ n ¤ Nh ’ 1, where |δk | ¤ µ, 0 ¤ k ¤ Nh , w0 = y0 and wk ,

k = 1, . . . , p, are p initial values generated by using another numerical

scheme.

Theorem 11.4 (Equivalence of zero-stability and root condition)

For a consistent multistep method, the root condition is equivalent to zero-

stability.

Proof. Let us begin by proving that the root condition is necessary for the zero-

stability to hold. We proceed by contradiction and assume that the method is

zero-stable and there exists a root ri which violates the root condition.

Since the method is zero-stable, condition (11.58) must be satis¬ed for any

Cauchy problem, in particular for the problem y (t) = 0 with y(0) = 0, whose

(h)

solution is, clearly, the null function. Similarly, the solution un of (11.60) with

(h)

f = 0 and wk = 0 for k = 0, . . . , p is identically zero.

Consider ¬rst the case |ri | > 1. Then, de¬ne

±

µri if ri ∈ R,

n

δn =

µ(r + r )n if r ∈ C,

¯i

i i

(h)

for µ > 0. It is simple to check that the sequence zn = δn for n = 0, 1, . . .

(h)

is a solution of (11.59) with initial conditions zk = δk and that |δk | ¤ µ for

¯

k = 0, 1, . . . , p. Let us now choose t in (t0 , t0 + T ) and let xn be the nearest grid

(h) (h)

node to t. Clearly, n is the integral part of t/h and limh’0 |zn | = limh’0 |un ’

¯ ¯

(h) (h) (h)

zn | ’ ∞ as h ’ 0. This proves that |un ’ zn | cannot be uniformly bounded

with respect to h as h ’ 0, which contradicts the assumption that the method

is zero-stable.

A similar proof can be carried out if |ri | = 1 but has multiplicity greater than

1, provided that one takes into account the form of the solution (11.33).

Let us now prove that the root condition is su¬cient for method (11.45) to

(h) (h)

be zero-stable. Recalling (11.46) and denoting by zn+j and un+j the solutions

to (11.59) and (11.60), respectively, for j ≥ 1, it turns out that the function

(h) (h) (h)

wn+j = zn+j ’ un+j satis¬es the following di¬erence equation

p+1

(h)

n = 0, . . . , Nh ’ (p + 1),

±j wn+j = •n+p+1 , (11.61)

j=0

having set

p+1

(h) (h)

βj f (tn+j , zn+j ) ’ f (tn+j , un+j ) + hδn+p+1 .

•n+p+1 = h (11.62)

j=0

11.6 Analysis of Multistep Methods 497

(n)

Denote by ψj a sequence of fundamental solutions to the homogeneous equa-

tion associated with (11.61). Recalling (11.42), the general solution of (11.61) is

given by

p n

(h) (h) (n) (n’l+p)

wn = wj ψj + ψp •l , n = p + 1, . . .

j=0 l=p+1

The following result expresses the connection between the root condition and the

boundedness of the solution of a di¬erence equation (for the proof, see [Gau97],

Theorem 6.3.2).

Lemma 11.3 There exists a constant M > 0 for any solution {un } of the dif-

ference equation (11.28) such that

n

|un | ¤ M |uj | + |•l |

max , n = 0, 1, . . . (11.63)

j=0,... ,k’1

l=k

i¬ the root condition is satis¬ed for the polynomial (11.30), i.e., (11.56) holds for

the zeros of the polynomial (11.30).

(n)

Let us now recall that, for any j, {ψj } is solution of a homogeneous di¬erence

(i)

equation whose initial data are ψj = δij , i, j = 0, . . . , p. On the other hand, for

(n’l+p)

any l, ψp is solution of a di¬erence equation with zero initial conditions and

right-hand sides equal to zero except for the one corresponding to n = l which is

(p)

ψp = 1.

Therefore, Lemma 11.3 can be applied in both cases so we can conclude that

(n) (n’l+p)

there exists a constant M > 0 such that |ψj | ¤ M and |ψp | ¤ M,

uniformly with respect to n and l. The following estimate thus holds

±

n

(h) (h)

|wn | ¤ M (p + 1) max |wj | + |•l | , n = 0, 1, . . . , Nh . (11.64)

j=0,... ,p

l=p+1

If L denotes the Lipschitz constant of f , from (11.62) we get

p+1

(h)

|•n+p+1 | ¤ h |βj |L |wn+j | + h|δn+p+1 |.

max

j=0,... ,p+1

j=0

|βj | and ∆[q,r] = max |δj+q |, q and r being some integers

Let β = max

j=0,... ,p+1 j=q,... ,r

with q ¤ r. From (11.64), the following estimate is therefore obtained

±

p+1

n

(h) (h)

|wn | ¤ |wl’p’1+j | + Nh h∆[p+1,n]

M (p + 1)∆[0,p] + hβL

l=p+1 j=0

n

¤ |wm | + T ∆[p+1,n]

(h)

M (p + 1)∆[0,p] + hβL(p + 2) .

m=0

498 11. Numerical Solution of Ordinary Di¬erential Equations

Let Q = 2(p + 2)βLM and h0 = 1/Q, so that 1 ’ h Q ≥ if h ¤ h0 . Then

1

2 2

1 (h) (h)

|wn | ¤ |wn |(1 ’ h Q )

2

2

n’1

¤ |wm | + T ∆[p+1,n]

(h)

M (p + 1)∆[0,p] + hβL(p + 2) .

m=0

Letting R = 2M (p + 1)∆[0,p] + T ∆[p+1,n] , we ¬nally obtain

n’1

|wn | ¤ hQ |wm | + R.

(h) (h)

m=0

(h)

Applying Lemma 11.2 with the following identi¬cations: •n = |wn |, g0 = R,

ps = 0 and ks = hQ for any s = 0, . . . , n ’ 1, yields

|wn | ¤ 2M eT Q (p + 1)∆[0,p] + T ∆[p+1,n] .

(h)

(11.65)

Method (11.45) is thus zero-stable for any h ¤ h0 . 3

Theorem 11.4 allows for characterizing the stability behavior of several

families of discretization methods.

In the special case of consistent one-step methods, the polynomial ρ ad-

mits only the root r0 = 1. They thus automatically satisfy the root condition

and are zero-stable.

For the Adams methods (11.49), the polynomial ρ is always of the form

ρ(r) = rp+1 ’ rp . Thus, its roots are r0 = 1 and r1 = 0 (with multiplicity

p) so that all Adams methods are zero-stable.

Also the midpoint method (11.43) and Simpson method (11.44) are zero-

stable: for both of them, the ¬rst characteristic polynomial is ρ(r) = r2 ’ 1,

so that r0 = 1 and r1 = ’1.

Finally, the BDF methods of Section 11.5.2 are zero-stable provided that

p ¤ 5, since in such a case the root condition is satis¬ed (see [Cry73]).

We are in position to give the following convergence result.

Theorem 11.5 (Convergence) A consistent multistep method is conver-

gent i¬ it satis¬es the root condition and the error on the initial data tends

to zero as h ’ 0. Moreover, the method converges with order q if it has

order q and the error on the initial data tends to zero as O(hq ).

Proof. Suppose that the MS method is consistent and convergent. To prove that

the root condition is satis¬ed, we refer to the problem y (t) = 0 with y(0) = 0

and on the interval I = (0, T ). Convergence means that the numerical solution

{un } must tend to the exact solution y(t) = 0 for any converging set of initial

data uk , k = 0, . . . , p, i.e. max |uk | ’ 0 as h ’ 0. From this observation, the

k=0,... ,p

proof follows by contradiction along the same lines as the proof of Theorem 11.4,

where the parameter µ is now replaced by h.

11.6 Analysis of Multistep Methods 499

Let us now prove that consistency, together with the root condition, implies

convergence under the assumption that the error on the initial data tends to zero

(h)

as h ’ 0. We can apply Theorem 11.4, setting un = un (approximate solution

(h)

of the Cauchy problem) and zn = yn (exact solution), and from (11.47) it turns

out that δm = „m (h). Then, due to (11.65), for any n ≥ p + 1 we obtain

|un ’ yn | ¤ 2M eT Q (p + 1) max |uj ’ yj | + T |„j (h)| .

max

j=0,... ,p j=p+1,... ,n

Convergence holds by noticing that the right-hand side of this inequality tends

to zero as h ’ 0. 3

A remarkable consequence of the above theorem is the following equivalence