<< . .

. 70
( : 95)



. . >>

un+k + ±k’j (n)un+k’j = 0, n = 0, 1, . . . . (11.38)
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

<< . .

. 70
( : 95)



. . >>