0.5

0

’4 ’3 ’2 ’1 0 1 2

FIGURE 11.10. Regions of absolute stability for s-stage explicit RK methods,

with s = 1, . . . , 4. The plot only shows the portion Im(h») ≥ 0 since the regions

are symmetric about the real axis

We ¬nally notice that the regions of absolute stability for explicit RK meth-

ods can fail to be connected; an example is given in Exercise 14.

11.9 Systems of ODEs

Let us consider the system of ¬rst-order ODEs

y = F(t, y), (11.78)

where F : R — Rn ’ Rn is a given vector function and y ∈ Rn is the

solution vector which depends on n arbitrary constants set by the n initial

518 11. Numerical Solution of Ordinary Di¬erential Equations

conditions

y(t0 ) = y0 . (11.79)

Let us recall the following property (see [PS91], p. 209).

Property 11.5 Let F : R — Rn ’ Rn be a continuous function on D =

[t0 , T ] — Rn , with t0 and T ¬nite. Then, if there exists a positive constant

L such that

F(t, y) ’ F(t, y) ¤ L y ’ y

¯ ¯ (11.80)

holds for any (t, y) and (t, y) ∈ D, then, for any y0 ∈ Rn there exists a

¯

unique y, continuous and di¬erentiable with respect to t for any (t, y) ∈ D,

which is a solution of the Cauchy problem (11.78)-(11.79).

Condition (11.80) expresses the fact that F is Lipschitz continuous with

respect to the second argument.

It is seldom possible to write out in closed form the solution to system

(11.78). A special case is where the system takes the form

y (t) = Ay(t), (11.81)

with A∈ Rn—n . Assume that A has n distinct eigenvalues »j , j = 1, . . . , n;

therefore, the solution y can be written as

n

Cj e»j t vj ,

y(t) = (11.82)

j=1

where C1 , . . . , Cn are some constants and {vj } is a basis formed by the

eigenvectors of A, associated with the eigenvalues »j for j = 1, . . . , n. The

solution is determined by setting n initial conditions.

From the numerical standpoint, the methods introduced in the scalar

case can be extended to systems. A delicate matter is how to generalize the

theory developed about absolute stability.

With this aim, let us consider system (11.81). As previously seen, the

property of absolute stability is concerned with the behavior of the numer-

ical solution as t grows to in¬nity, in the case where the solution of problem

(11.78) satis¬es

y(t) ’ 0 as t ’ ∞. (11.83)

Condition (11.83) is satis¬ed if all the real parts of the eigenvalues of A are

negative since this ensures that

e»j t = eRe»j t (cos(Im»j ) + i sin(Im»i )) ’ 0, as t ’ ∞, (11.84)

11.10 Sti¬ Problems 519

from which (11.83) follows recalling (11.82). Since A has n distinct eigen-

values, there exists a nonsingular matrix Q such that Λ = Q’1 AQ, Λ being

the diagonal matrix whose entries are the eigenvalues of A (see Section 1.8).

Introducing the auxiliary variable z = Q’1 y, the initial system can there-

fore be transformed into

z = Λz. (11.85)

Since Λ is a diagonal matrix, the results holding in the scalar case immedi-

ately apply to the vector case as well, provided that the analysis is repeated

on all the (scalar) equations of system (11.85).

11.10 Sti¬ Problems

Consider a nonhomogeneous linear system of ODEs with constant coe¬-

cients

y (t) = Ay(t) + •(t), with A ∈ Rn—n , •(t) ∈ Rn ,

and assume that A has n distinct eigenvalues »j , j = 1, . . . , n. Then

n

Cj e»j t vj + ψ(t)

y(t) =

j=1

where C1 , . . . , Cn , are n constants, {vj } is a basis formed by the eigenvec-

tors of A and ψ(t) is a particular solution of the ODE at hand. Throughout

the section, we assume that Re»j < 0 for all j.

As t ’ ∞, the solution y tends to the particular solution ψ. We can

therefore interpret ψ as the steady-state solution (that is, after an in¬nite

n

Cj e»j t as being the transient solution (that is, for t ¬nite).

time) and

j=1

Assume that we are interested only in the steady-state. If we employ a

numerical scheme with a bounded region of absolute stability, the stepsize h

is subject to a constraint that depends on the maximum module eigenvalue

of A. On the other hand, the greater this module, the shorter the time

interval where the corresponding component in the solution is meaningful.

We are thus faced with a sort of paradox: the scheme is forced to employ

a small integration stepsize to track a component of the solution that is

virtually ¬‚at for large values of t.

Precisely, if we assume that

σ ¤ Re»j ¤ „ < 0, ∀j = 1, . . . , n (11.86)

and introduce the sti¬ness quotient rs = σ/„ , we say that a linear system

of ODEs with constant coe¬cients is sti¬ if the eigenvalues of matrix A all

have negative real parts and rs 1.

520 11. Numerical Solution of Ordinary Di¬erential Equations

However, referring only to the spectrum of A to characterize the sti¬ness of

a problem might have some drawbacks. For instance, when „ 0, the sti¬-

ness quotient can be very large while the problem appears to be “genuinely”

sti¬ only if |σ| is very large. Moreover, enforcing suitable initial conditions

can a¬ect the sti¬ness of the problem (for example, selecting the data in

such a way that the constants multiplying the “sti¬” components of the

solution vanish).

For this reason, several authors ¬nd the previous de¬nition of a sti¬

problem unsatisfactory, and, on the other hand, they agree on the fact that

it is not possible to exactly state what it is meant by a sti¬ problem. We

limit ourselves to quoting only one alternative de¬nition, which is of some

interest since it focuses on what is observed in practice to be a sti¬ problem.

De¬nition 11.14 (from [Lam91], p. 220) A system of ODEs is sti¬ if,

when approximated by a numerical scheme characterized by a region of

absolute stability with ¬nite size, it forces the method, for any initial con-

dition for which the problem admits a solution, to employ a discretization

stepsize excessively small with respect to the smoothness of the exact so-

lution.

From this de¬nition, it is clear that no conditionally absolute stable method

is suitable for approximating a sti¬ problem. This prompts resorting to

implicit methods, such as MS or RK, which are more expensive than explicit

schemes, but have regions of absolute stability of in¬nite size. However, it

is worth recalling that, for nonlinear problems, implicit methods lead to

nonlinear equations, for which it is thus crucial to select iterative numerical

methods free of limitations on h for convergence.

For instance, in the case of MS methods, we have seen that using ¬xed-

point iterations sets the constraint (11.68) on h in terms of the Lipschitz

constant L of f . In the case of a linear system this constraint is

L ≥ max |»i |,

i=1,... ,n

so that (11.68) would imply a strong limitation on h (which could even

be more stringent than those required for an explicit scheme to be stable).

One way of circumventing this drawback consists of resorting to Newton™s

method or some variants. The presence of Dahlquist barriers imposes a

strong limitation on the use of MS methods, the only exception being BDF

methods, which, as already seen, are θ-stable for p ¤ 5 (for a larger number

of steps they are even not zero-stable). The situation becomes de¬nitely

more favorable if implicit RK methods are considered, as observed at the

end of Section 11.8.4.

The theory developed so far holds rigorously if the system is linear. In

the nonlinear case, let us consider the Cauchy problem (11.78), where the

11.11 Applications 521

function F : R — Rn ’ Rn is assumed to be di¬erentiable. To study its

stability a possible strategy consists of linearizing the system as

y (t) = F(„, y(„ )) + JF („, y(„ )) [y(t) ’ y(„ )] ,

in a neighborhood („, y(„ )), where „ is an arbitrarily chosen value of t

within the time integration interval.

The above technique might be dangerous since the eigenvalues of JF do

not su¬ce in general to describe the behavior of the exact solution of the

original problem. Actually, some counterexamples can be found where:

1. JF has complex conjugate eigenvalues, while the solution of (11.78)

does not exhibit oscillatory behavior;

2. JF has real nonnegative eigenvalues, while the solution of (11.78) does

not grow monotonically with t;

3. JF has eigenvalues with negative real parts, but the solution of (11.78)

does not decay monotonically with t.

As an example of the case at item 3. let us consider the system of

ODEs

®

1 2

’

2t

t3

y =° »y = A(t)y.

t 1

’ ’

2 2t

For t ≥ 1 its solution is

t’3/2 2t’3/2 log t

y(t) = C1 + C2

’ 1 t1/2 t1/2 (1 ’ log t)

2

whose Euclidean norm diverges monotonically for t > (12)1/4 1.86

when C1 = 1, C2 = 0, whilst the eigenvalues of A(t), equal to (’1 ±

2i)/(2t), have negative real parts.

Therefore, the nonlinear case must be tackled using ad hoc techniques, by

suitably reformulating the concept of stability itself (see [Lam91], Chapter

7).

11.11 Applications

We consider two examples of dynamical systems that are well-suited to

checking the performances of several numerical methods introduced in the

previous sections.

522 11. Numerical Solution of Ordinary Di¬erential Equations

11.11.1 Analysis of the Motion of a Frictionless Pendulum

Let us consider the frictionless pendulum in Figure 11.11 (left), whose

motion is governed by the following system of ODEs

y1 = y2 ,

(11.87)

’K sin(y1 ),

y2 =

for t > 0, where y1 (t) and y2 (t) represent the position and angular velocity

of the pendulum at time t, respectively, while K is a positive constant

depending on the geometrical-mechanical parameters of the pendulum. We

consider the initial conditions: y1 (0) = θ0 , y2 (0) = 0.

1/2

πK

A

1/2

’ πK

A™

y1

weight

FIGURE 11.11. Left: frictionless pendulum; right: orbits of system (11.87) in the

phase space

Denoting by y = (y1 , y2 )T the solution to system (11.87), this admits

in¬nitely many equilibrium conditions of the form y = (nπ, 0)T for n ∈ Z,

corresponding to the situations where the pendulum is vertical with zero

velocity. For n even, the equilibrium is stable, while for n odd it is unstable.

These conclusions can be drawn by analyzing the linearized system

0 1 0 1

y = Ae y = y, y = Ao y = y.

’K 0 K 0

√

If n is even, matrix Ae has complex conjugate eigenvalues »1,2 = ±i K

√

and associated eigenvectors y1,2 = (“i/ K, 1)T , while for n odd, Ao

√

real and opposite eigenvalues »1,2 = ± K and eigenvectors y1,2 =

has√

(1/ K, “1)T .

Let us consider two di¬erent sets of initial data: y(0) = (θ0 , 0)T and

y(0) = (π + θ0 , 0)T , where |θ0 | 1. The solutions of the corresponding

linearized system are, respectively,

√ √

y1 (t) = θ0 cos( Kt) y1 (t) = (π + θ0 ) cosh( Kt)

√ √ √ √

,

y2 (t) = ’ Kθ0 sin( Kt) y2 (t) = K(π + θ0 ) sinh( Kt),

11.11 Applications 523

and will be henceforth denoted as “stable” and “unstable”, respectively, for

reasons that will be clear later on. To these solutions we associate in the

plane (y1 , y2 ), called the phase space, the following orbits (i.e., the graphs

obtained plotting the curve (y1 (t), y2 (t)) in the phase space).

2 2

y1 y

√2

+ = 1, (stable case)

θ0 Kθ0

2 2

y1 y2

√

’ = 1, (unstable case).

π + θ0 K(π + θ0 )

√

In the stable case, the orbits are ellipses with period 2π/ K and are cen-

tered at (0, 0)T , while in the unstable case they are √

hyperbolae centered at

(0, 0) and asymptotic to the straight lines y2 = ± Ky1 .

T

The complete picture of the motion of the pendulum in the phase space

is shown in Figure 11.11 (right). Notice that, letting v = |y2 | √ ¬xing

and

the initial position y1 (0) = 0, there exists a limit value vL = 2 K which

corresponds in the ¬gure to the points A and A™. For v(0) < vL , the orbits

are closed, while for v(0) > vL they are open, corresponding to a continuous

revolution of the pendulum, with in¬nite passages (with periodic and non

null velocity) through the two equilibrium positions y1 = 0 and y1 = π.

The limit case v(0) = vL yields a solution such that, thanks to the total

energy conservation principle, y2 = 0 when y1 = π. Actually, these two

values are attained asymptotically only as t ’ ∞.

The ¬rst-order nonlinear di¬erential system (11.87) has been numerically

solved using the forward Euler method (FE), the midpoint method (MP)

and the Adams-Bashforth second-order scheme (AB). In Figure 11.12 we

show the orbits in the phase space that have been computed by the two

methods on the time interval (0, 30) and taking K = 1 and h = 0.1. The

crosses denote initial conditions.

As can be noticed, the orbits generated by FE do not close. This kind

of instability is due to the fact that the region of absolute stability of the

FE method completely excludes the imaginary axis. On the contrary, the

MP method describes accurately the closed system orbits due to the fact

that its region of asymptotic stability (see Section 11.6.4) includes pure

imaginary eigenvalues in the neighborhood of the origin of the complex

plane. It must also be noticed that the MP scheme gives rise to oscillating

solutions as v0 gets larger. The second-order AB method, instead, describes

correctly all kinds of orbits.

11.11.2 Compliance of Arterial Walls

An arterial wall subject to blood ¬‚ow can be modelled by a compliant

circular cylinder of length L and radius R0 with walls made by an incom-

pressible, homogeneous, isotropic, elastic tissue of thickness H. A simple

524 11. Numerical Solution of Ordinary Di¬erential Equations

2

0

’2

’10 ’5 0 5 10

2

0

’2

’10 ’5 0 5 10

2

0

’2

’10 ’5 0 5 10

FIGURE 11.12. Orbits for system (11.87) in the case K = 1 and h = 0.1, com-

puted using the FE method (upper plot), the MP method (central plot) and the

AB method (lower plot), respectively. The initial conditions are θ0 = π/10 and

v0 = 0 (thin solid line), v0 = 1 (dashed line), v0 = 2 (dash-dotted line) and

v0 = ’2 (thick solid line)

model describing the mechanical behavior of the walls interacting with the

blood ¬‚ow is the so called “independent-rings” model according to which

the vessel wall is regarded as an assembly of rings which are not in¬‚uenced

one by the others.

This amounts to neglecting the longitudinal (or axial) inner actions along

the vessel, and to assuming that the walls can deform only in the radial

direction. Thus, the vessel radius R is given by R(t) = R0 + y(t), where y is

the radial deformation of the ring with respect to a reference radius R0 and

t is the time variable. The application of Newton™s law to the independent-

ring system yields the following equation modeling the time mechanical

11.11 Applications 525

behavior of the wall

y (t) + βy (t) + ±y(t) = γ(p(t) ’ p0 ) (11.88)

2

where ± = E/(ρw R0 ), γ = 1/(ρw H) and β is a positive constant. The

physical parameters ρw and E denote the vascular wall density and the

Young modulus of the vascular tissue, respectively. The function p ’ p0

is the forcing term acting on the wall due to the pressure drop between

the inner part of the vessel (where the blood ¬‚ows) and its outer part

(surrounding organs). At rest, if p = p0 , the vessel con¬guration coincides

with the undeformed circular cylinder having radius equal exactly to R0

(y = 0).

Equation (11.88) can be formulated as y (t) = Ay(t) + b(t) where y =

(y, y )T , b = (0, ’γ(p ’ p0 ))T and

0 1

A= . (11.89)

’± ’β

√

The eigenvalues of A are »± = (’β ± β 2 ’ 4±)/2; therefore, if β ≥ 2 ±

both the eigenvalues are real and negative and the system is asymptotically