j=1

The goal is to transform it into an ODE by means of a function F , called

the generating function of the equation (11.38). F depends on the real

variable t and is derived as follows. We require that the n-th coe¬cient

of the Taylor series of F around t = 0 can be written as γn un , for some

unknown constant γn , so that

∞

γn un tn .

F (t) = (11.39)

n=0

The coe¬cients {γn } are unknown and must be determined in such a way

that

®

∞

k k

°un+k + ±k’j (n)un+k’j » tn ,

(k’j)

cj F (t) = (11.40)

n=0

j=0 j=1

486 11. Numerical Solution of Ordinary Di¬erential Equations

where cj are suitable unknown constants not depending on n. Note that

owing to (11.39) we obtain the ODE

k

cj F (k’j) (t) = 0

j=0

to which we must add the initial conditions F (j) (0) = γj uj for j = 0, . . . , k’

1. Once F is available, it is simple to recover un through the de¬nition of

F itself.

Example 11.6 Consider the di¬erence equation

(n + 2)(n + 1)un+2 ’ 2(n + 1)un+1 ’ 3un = 0, n = 0, 1, . . . (11.41)

with the initial conditions u0 = u1 = 2. We look for a generating function of the

form (11.39). By term-to-term derivation of the series, we get

∞ ∞

γn n(n ’ 1)un tn’2 ,

n’1

F (t) = γn nun t , F (t) =

n=0 n=0

and, after some algebra, we ¬nd

∞ ∞

n’1

γn+1 (n + 1)un+1 tn ,

F (t) = γn nun t =

n=0 n=0

∞ ∞

γn n(n ’ 1)un t n’2

γn+2 (n + 2)(n + 1)un+2 tn .

F (t) = =

n=0 n=0

As a consequence, (11.40) becomes

∞ ∞ ∞

(n + 1)(n + 2)un+2 t ’ 2 (n + 1)un+1 t ’ 3

n n

un tn

n=0 n=0 n=0

∞ ∞ ∞

n n

γn un tn ,

= c0 γn+2 (n + 2)(n + 1)un+2 t + c1 γn+1 (n + 1)un+1 t + c2

n=0 n=0 n=0

so that, equating both sides, we ¬nd

γn = 1 ∀n ≥ 0, c0 = 1, c1 = ’2, c2 = ’3.

We have thus associated with the di¬erence equation the following ODE with

constant coe¬cients

F (t) ’ 2F (t) ’ 3F (t) = 0,

with the initial condition F (0) = F (0) = 2. The n-th coe¬cient of the solution

F (t) = e3t + e’t is

1 (n) 1

[(’1)n + 3n ] ,

F (0) =

n! n!

•

so that un = (1/n!) [(’1)n + 3n ] is the solution of (11.41).

11.5 Multistep Methods 487

The nonhomogeneous case (11.28) can be tackled by searching for solutions

of the form

un = u(0) + u(•) ,

n n

(0) (•)

where un is the solution of the associated homogeneous equation and un

is a particular solution of the nonhomogeneous equation. Once the solution

of the homogeneous equation is available, a general technique to obtain

the solution of the nonhomogeneous equation is based on the method of

variation of parameters, combined with a reduction of the order of the

di¬erence equation (see [BO78]).

In the special case of di¬erence equations with constant coe¬cients, with

•n of the form cn Q(n), where c is a constant and Q is a polynomial of degree

p with respect to the variable n, a possible approach is that of undetermined

coe¬cients, where one looks for a particular solution that depends on some

undetermined constants and has a known form for some classes of right

sides •n . It su¬ces to look for a particular solution of the form

u(•) = cn (bp np + bp’1 np’1 + . . . + b0 ),

n

(•)

where bp , . . . , b0 are constants to be determined in such a way that un is

actually a solution of (11.28).

Example 11.7 Consider the di¬erence equation un+3 ’un+2 +un+1 ’un = 2n n2 .

The particular solution is of the form un = 2n (b2 n2 + b1 n + b0 ). Substituting this

solution into the equation, we ¬nd 5b2 n2 +(36b2 +5b1 )n+(58b2 +18b1 +5b0 ) = n2 ,

from which, recalling the principle of identity for polynomials, one gets b2 = 1/5,

b1 = ’36/25 and b0 = 358/125. •

Analogous to the homogeneous case, it is possible to express the solution

of (11.28) as

k’1 n

(n) (n’l+k’1)

un = uj ψj + •l ψk’1 , n = 0, 1, . . . (11.42)

j=0 l=k

(i)

where we de¬ne ψk’1 ≡ 0 for all i < 0 and •j ≡ 0 for all j < k.

11.5 Multistep Methods

Let us now introduce some examples of multistep methods (shortly, MS).

De¬nition 11.7 (q-steps methods) A q-step method (q ≥ 1) is one

which, ∀n ≥ q ’ 1, un+1 depends on un+1’q , but not on the values uk

with k < n + 1 ’ q.

488 11. Numerical Solution of Ordinary Di¬erential Equations

A well-known two-step explicit method can be obtained by using the

centered ¬nite di¬erence (10.61) to approximate the ¬rst order derivative

in (11.1). This yields the midpoint method

n≥1

un+1 = un’1 + 2hfn , (11.43)

where u0 = y0 , u1 is to be determined and fk denotes the value f (tk , uk ).

An example of an implicit two-step scheme is the Simpson method, ob-

tained from (11.2) with t0 = tn’1 and t = tn+1 and by using the Cavalieri-

Simpson quadrature rule to approximate the integral of f

h

n≥1

un+1 = un’1 + [fn’1 + 4fn + fn+1 ], (11.44)

3

where u0 = y0 , and u1 is to be determined.

From these examples, it is clear that a multistep method requires q initial

values u0 , . . . , uq’1 for “taking o¬”. Since the Cauchy problem provides

only one datum (u0 ), one way to assign the remaining values consists of

resorting to explicit one-step methods of high order. An example is given by

Heun™s method (11.10), other examples are provided by the Runge-Kutta

methods, which will be introduced in Section 11.8.

In this section we deal with linear multistep methods

p p

un+1 = aj un’j + h bj fn’j + hb’1 fn+1 , n = p, p + 1, . . . (11.45)

j=0 j=0

which are p + 1-step methods, p ≥ 0. For p = 0, we recover one-step

methods.

The coe¬cients aj , bj are real and fully identify the scheme; they are

such that ap = 0 or bp = 0. If b’1 = 0 the scheme is implicit, otherwise it

is explicit.

We can reformulate (11.45) as follows

p+1 p+1

βs f (tn+s , un+s ), n = 0, 1, . . . , Nh ’ (p + 1) (11.46)

±s un+s = h

s=0 s=0

having set ±p+1 = 1, ±s = ’ap’s for s = 0, . . . , p and βs = bp’s for

s = 0, . . . , p+1. Relation (11.46) is a special instance of the linear di¬erence

equation (11.28), where we set k = p + 1 and •n+j = hβj f (tn+j , un+j ), for

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

Also for MS methods we can characterize consistency in terms of the local

truncation error, according to the following de¬nition.

De¬nition 11.8 The local truncation error (LTE) „n+1 (h) introduced by

the multistep method (11.45) at tn+1 (for n ≥ p) is de¬ned through the

11.5 Multistep Methods 489

following relation

®

p p

h„n+1 (h) = yn+1 ’ ° bj yn’j » , n ≥ p, (11.47)

aj yn’j + h

j=0 j=’1

where yn’j = y(tn’j ) and yn’j = y (tn’j ) for j = ’1, . . . , p.

Analogous to one-step methods, the quantity h„n+1 (h) is the residual gen-

erated at tn+1 if we pretend that the exact solution “satis¬es” the numerical

scheme. Letting „ (h) = max|„n (h)|, we have the following de¬nition.

n

De¬nition 11.9 (Consistency) The multistep method (11.45) is consis-

tent if „ (h) ’ 0 as h ’ 0. Moreover, if „ (h) = O(hq ), for some q ≥ 1, then

the method is said to have order q.

A more precise characterization of the LTE can be given by introducing the

following linear operator L associated with the linear MS method (11.45)

p p

L[w(t); h] = w(t + h) ’ aj w(t ’ jh) ’ h bj w (t ’ jh), (11.48)

j=0 j=’1

where w ∈ C 1 (I) is an arbitrary function. Notice that the LTE is exactly

L[y(tn ); h]. If we assume that w is su¬ciently smooth and expand w(t’jh)

and w (t ’ jh) about t ’ ph, we obtain

L[w(t); h] = C0 w(t ’ ph) + C1 hw(1) (t ’ ph) + . . . + Ck hk w(k) (t ’ ph) + . . .

Consequently, if the MS method has order q and y ∈ C q+1 (I), we obtain

„n+1 (h) = Cq+1 hq+1 y (q+1) (tn’p ) + O(hq+2 ).

The term Cq+1 hq+1 y (q+1) (tn’p ) is the so-called principal local truncation

error (PLTE) while Cq+1 is the error constant. The PLTE is widely em-

ployed in devising adaptive strategies for MS methods (see [Lam91], Chap-

ter 3).

Program 92 provides an implementation of the multistep method in the

form (11.45) for the solution of a Cauchy problem on the interval (t0 , T ).

The input parameters are: the column vector a containing the p + 1 co-

e¬cients ai ; the column vector b containing the p + 2 coe¬cients bi ; the

discretization stepsize h; the vector of initial data u0 at the corresponding

time instants t0; the macros fun and dfun containing the functions f and

‚f /‚y. If the MS method is implicit, a tolerance tol and a maximum num-

ber of admissible iterations itmax must be provided. These two parameters

monitor the convergence of Newton™s method that is employed to solve the

490 11. Numerical Solution of Ordinary Di¬erential Equations

nonlinear equation (11.45) associated with the MS method. In output the

code returns the vectors u and t containing the computed solution at the

time instants t.

Program 92 - multistep : Linear multistep methods

function [t,u] = multistep (a,b,tf,t0,u0,h,fun,dfun,tol,itmax)

y = u0; t = t0; f = eval (fun); p = length(a) - 1; u = u0;

nt = ¬x((tf - t0 (1) )/h);

for k = 1:nt

lu=length(u);

G = a™ *u (lu:-1:lu-p)+ h * b(2:p+2)™ * f(lu:-1:lu-p);

lt = length(t0); t0 = [t0; t0(lt)+h]; unew = u (lu);

t = t0 (lt+1); err = tol + 1; it = 0;

while (err > tol) & (it <= itmax)

y = unew; den = 1 - h * b (1) * eval(dfun);

fnew = eval (fun);

if den == 0

it = itmax + 1;

else

it = it + 1;

unew = unew - (unew - G - h * b (1) * fnew)/den;

err = abs (unew - y);

end

end

u = [u; unew]; f = [f; fnew];

end

t = t0;

In the forthcoming sections we examine some families of multistep methods.

11.5.1 Adams Methods

These methods are derived from the integral form (11.2) through an ap-

proximate evaluation of the integral of f between tn and tn+1 . We suppose

that the discretization nodes are equally spaced, i.e., tj = t0 + jh, with

h > 0 and j ≥ 1, and then we integrate, instead of f , its interpolating

polynomial on p + 1 distinct nodes. The resulting schemes are thus consis-

tent by construction and have the following form

p

n ≥ p.

un+1 = un + h bj fn’j , (11.49)

j=’1

The interpolation nodes can be either:

1. tn , tn’1 , . . . , tn’p (in this case b’1 = 0 and the resulting method is

explicit);

or

11.5 Multistep Methods 491

2. tn+1 , tn , . . . , tn’p+1 (in this case b’1 = 0 and the scheme is implicit).

The implicit schemes are called Adams-Moulton methods, while the explicit

ones are called Adams-Bashforth methods.

Adams-Bashforth methods (AB)

Taking p = 0 we recover the forward Euler method, since the interpolating

polynomial of degree zero at node tn is given by Π0 f = fn . For p = 1, the

linear interpolating polynomial at the nodes tn’1 and tn is

fn’1 ’ fn

Π1 f (t) = fn + (t ’ tn ) .

tn’1 ’ tn

Since Π1 f (tn ) = fn and Π1 f (tn+1 ) = 2fn ’ fn’1 , we get

tn+1

h h

[Π1 f (tn ) + Π1 f (tn+1 )] = [3fn ’ fn’1 ] .

Π1 f (t) =

2 2

tn

Therefore, the two-step AB method is

h

[3fn ’ fn’1 ] .

un+1 = un + (11.50)

2

With a similar procedure, if p = 2, we ¬nd the three-step AB method

h

[23fn ’ 16fn’1 + 5fn’2 ] ,

un+1 = un +

12

while for p = 3 we get the four-step AB scheme

h

(55fn ’ 59fn’1 + 37fn’2 ’ 9fn’3 ) .

un+1 = un +

24

In general, q-step Adams-Bashforth methods have order q. The error con-

—

stants Cq+1 of these methods are collected in Table 11.1.

Adams-Moulton methods (AM)

If p = ’1, the Backward Euler scheme is recovered, while if p = 0, we

construct the linear polynomial interpolating f at the nodes tn and tn+1

to recover the Crank-Nicolson scheme (11.9). In the case of the two-step

method, the polynomial of degree 2 interpolating f at the nodes tn’1 , tn ,

tn+1 is generated, yielding the following scheme

h

[5fn+1 + 8fn ’ fn’1 ] .

un+1 = un + (11.51)

12

492 11. Numerical Solution of Ordinary Di¬erential Equations

The methods corresponding to p = 3 and 4 are respectively given by

h

(9fn+1 + 19fn ’ 5fn’1 + fn’2 )

un+1 = un +

24

h

(251fn+1 + 646fn ’ 264fn’1 + 106fn’2 ’ 19fn’3 ) .

un+1 = un +

720

The q-step Adams-Moulton methods have order q + 1 and their error con-

stants Cq+1 are summarized in Table 11.1.

— —

q Cq+1 Cq+1 q Cq+1 Cq+1

’1 ’ 24

1 3 1

1 3

2 2 8

’ 12 ’ 720

5 1 251 19

2 4

12 720

TABLE 11.1. Error constants for Adams-Bashforth methods (having order q) and

Adams-Moulton methods (having order q + 1)

11.5.2 BDF Methods

The so-called backward di¬erentiation formulae (henceforth denoted by

BDF) are implicit MS methods derived from a complementary approach

to the one followed for the Adams methods. In fact, for the Adams meth-

ods we have resorted to numerical integration for the source function f ,

whereas in BDF methods we directly approximate the value of the ¬rst

derivative of y at node tn+1 through the ¬rst derivative of the polynomial

interpolating y at the p + 1 nodes tn+1 , tn , . . . , tn’p+1 .

By doing so, we get schemes of the form

p

un+1 = aj un’j + hb’1 fn+1

j=0

with b’1 = 0. Method (11.8) represents the most elementary example,

corresponding to the coe¬cients a0 = 1 and b’1 = 1.

We summarize in Table 11.2 the coe¬cients of BDF methods that are

zero-stable. In fact, we shall see in Section 11.6.3 that only for p ¤ 5 are

BDF methods zero-stable (see [Cry73]).

11.6 Analysis of Multistep Methods

Analogous to what has been done for one-step methods, in this section

we provide algebraic conditions that ensure consistency and stability of

multistep methods.

11.6 Analysis of Multistep Methods 493