<< . .

. 75
( : 95)

. . >>


’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


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
Cj e»j t vj ,
y(t) = (11.82)

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¬-
y (t) = Ay(t) + •(t), with A ∈ Rn—n , •(t) ∈ Rn ,
and assume that A has n distinct eigenvalues »j , j = 1, . . . , n. Then
Cj e»j t vj + ψ(t)
y(t) =

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
Cj e»j t as being the transient solution (that is, for t ¬nite).
time) and
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-
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
® 
1 2

 2t 
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)

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

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 ,
’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.


’ πK

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 =
(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
+ = 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 .

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
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




’10 ’5 0 5 10




’10 ’5 0 5 10




’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)
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

<< . .

. 75
( : 95)

. . >>