<<

. 21
( 26)



>>


we get

K2 = fn + hc2 (fn,t + K1 fn,y ) + O(h2 ).

We have denoted by fn,z (for z = t or z = y) the partial derivative of f
with respect to z evaluated at (tn , yn ). Then

un+1 = yn + hfn (b1 + b2 ) + h2 c2 b2 (fn,t + fn fn,y ) + O(h3 ).

If we perform the same expansion on the exact solution, we ¬nd

h2 h2
yn + O(h3 ) = yn + hfn + (fn,t + fn fn,y ) + O(h3 ).
yn+1 = yn + hyn +
2 2
Forcing the coe¬cients in the two expansions above to agree, up to higher-
order terms, we obtain that the coe¬cients of the RK method must satisfy
1
b1 + b2 = 1, c2 b2 = 2 .
Thus, there are in¬nitely many 2-stage explicit RK methods with second-
order accuracy. Two examples are the Heun method (11.10) and the modi-
¬ed Euler method (11.91). Of course, with similar (and cumbersome) com-
putations in the case of higher-stage methods, and accounting for a higher
number of terms in Taylor™s expansion, one can generate higher-order RK
methods. For instance, retaining all the terms up to the ¬fth one, we get
scheme (11.73).


11.8.2 Stepsize Adaptivity for RK Methods
Since RK schemes are one-step methods, they are well-suited to adapting
the stepsize h, provided that an e¬cient estimator of the local error is
available. Usually, a tool of this kind is an a posteriori error estimator,
since the a priori local error estimates are too complicated to be used in
practice. The error estimator can be constructed in two ways:

- using the same RK method, but with two di¬erent stepsizes (typically 2h
and h);
- using two RK methods of di¬erent order, but with the same number s of
stages.

In the ¬rst case, if an RK method of order p is being used, one pretends
that, starting from an exact datum un = yn (which would not be available
if n ≥ 1), the local error at tn+1 is less than a ¬xed tolerance. The following
relation holds

yn+1 ’ un+1 = ¦(yn )hp+1 + O(hp+2 ), (11.74)

where ¦ is an unknown function evaluated at yn . (Notice that, in this
special case, yn+1 ’ un+1 = h„n+1 (h)).
11.8 Runge-Kutta (RK) Methods 513

Carrying out the same computation with a stepsize of 2h, starting from
tn’1 , and denoting by un+1 the computed solution, yields

yn+1 ’ un+1 = ¦(yn’1 )(2h)p+1 + O(hp+2 ) = ¦(yn )(2h)p+1 + O(hp+2 )
(11.75)

having expanded also yn’1 with respect to tn . Subtracting (11.74) from
(11.75), we get

(2p+1 ’ 1)hp+1 ¦(yn ) = un+1 ’ un+1 + O(hp+2 ),

from which
un+1 ’ un+1
yn+1 ’ un+1 = E.
(2p+1 ’ 1)

If |E| is less than the ¬xed tolerance µ, the scheme moves to the next time
step, otherwise the estimate is repeated with a halved stepsize. In general,
the stepsize is doubled whenever |E| is less than µ/2p+1 .
This approach yields a considerable increase in the computational e¬ort,
due to the s ’ 1 extra functional evaluations needed to generate the value
un+1 . Moreover, if one needs to half the stepsize, the value un must also
be computed again.
An alternative that does not require extra functional evaluations consists
of using simultaneously two di¬erent RK methods with s stages, of order
p and p + 1, respectively, which share the same set of values Ki . These
methods are synthetically represented by the modi¬ed Butcher array

c A
2
bT
(11.76)
bT 2
2
ET

where the method of order p is identi¬ed by the coe¬cients c, A and b,
while that of order p + 1 is identi¬ed by c, A and b, and where E = b ’ b.
Taking the di¬erence between the approximate solutions at tn+1 pro-
duced by the two methods provides an estimate of the local truncation
error for the scheme of lower order. On the other hand, since the coe¬-
s
cients Ki coincide, this di¬erence is given by h i=1 Ei Ki and thus it does
not require extra functional evaluations.
Notice that, if the solution un+1 computed by the scheme of order p is
used to initialize the scheme at time step n + 2, the method will have order
p, as a whole. If, conversely, the solution computed by the scheme of order
p + 1 is employed, the resulting scheme would still have order p + 1 (exactly
as happens with predictor-corrector methods).
The Runge-Kutta Fehlberg method of fourth-order is one of the most
popular schemes of the form (11.76) and consists of a fourth-order RK
514 11. Numerical Solution of Ordinary Di¬erential Equations

scheme coupled with a ¬fth-order RK method (for this reason, it is known
as the RK45 method). The modi¬ed Butcher array for this method is shown
below
0 0 0 0 0 0 0
1 1
0 0 0 0 0
4 4
3 3 9
0 0 0 0
8 32 32

’ 7200
12 1932 7296
0 0 0
13 2197 2197 2197

’8 ’ 4104
439 3680 845
1 0 0
216 513

’ 27 ’ 3544 ’ 11
1 8 1859
2 0
2 2565 4104 40


’1
25 1408 2197
0 0
216 2565 4104 5

’ 50
16 6656 28561 9 2
0
135 12825 56430 55


’ 4275 ’ 75240
1 128 2197 1 2
0
360 50 55

This method tends to underestimate the error. As such, its use is not com-
pletely reliable when the stepsize h is large.

Remark 11.5 MATLAB provides a package tool funfun, which, besides
the two classical Runge-Kutta Fehlberg methods, RK23 (second-order and
third-order pair) and RK45 (fourth-order and ¬fth-order pair), also imple-
ments other methods suitable for solving sti¬ problems. These methods are
derived from BDF methods (see [SR97]) and are included in the MATLAB
program ode15s.


11.8.3 Implicit RK Methods
Implicit RK methods can be derived from the integral formulation of the
Cauchy problem (11.2). In fact, if a quadrature formula with s nodes in
(tn , tn+1 ) is employed to approximate the integral of f (which we assume,
for simplicity, to depend only on t), we get
tn+1
s
f („ ) d„ h bj f (tn + cj h)
j=1
tn

having denoted by bj the weights and by tn + cj h the quadrature nodes. It
can be proved (see [But64]) that for any RK formula (11.70)-(11.71), there
exists a correspondence between the coe¬cients bj , cj of the formula and
the weights and nodes of a Gauss quadrature rule.
In particular, the coe¬cients c1 , . . . , cs are the roots of the Legendre
polynomial Ls in the variable x = 2c ’ 1, so that x ∈ [’1, 1]. Once the
11.8 Runge-Kutta (RK) Methods 515

s coe¬cients cj have been found, we can construct RK methods of order
2s, by determining the coe¬cients aij and bj as being the solutions of the
linear systems
s
ck’1 aij = (1/k)ck , k = 1, 2, . . . , s,
i
j
j=1
s
ck’1 bj = 1/k, k = 1, 2, . . . , s.
j
j=1

The following families can be derived:

1. Gauss-Legendre RK methods, if Gauss-Legendre quadrature nodes are
used. These methods, for a ¬xed number of stages s, attain the maximum
possible order 2s. Remarkable examples are the one-stage method (implicit
midpoint method) of order 2
1 1
2 2
1 1
un+1 = un + hf tn + 2 h, 2 (un + un+1 ) ,
1

and the 2-stage method of order 4, described by the following Butcher array
√ √
3’ 3 3’2 3
1
6 4 12
√ √
3+ 3 3+2 3 1
6 12 4

1 1
2 2


2. Gauss-Radau methods, which are characterized by the fact that the
quadrature nodes include one of the two endpoints of the interval (tn , tn+1 ).
The maximum order that can be achieved by these methods is 2s ’ 1, when
s stages are used. Elementary examples correspond to the following Butcher
arrays

’ 12
1 5 1
3 12
3 1
1
0 1 11 4 4
, ,
1 1
3 1
4 4

and have order 1, 1 and 3, respectively. The Butcher array in the middle
represents the backward Euler method.

3. Gauss-Lobatto methods, where both the endpoints tn and tn+1 are quadra-
ture nodes. The maximum order that can be achieved using s stages is
2s ’ 2. We recall the methods of the family corresponding to the following
516 11. Numerical Solution of Ordinary Di¬erential Equations

Butcher arrays
’1
1 1
0 6 3 6
1
0 0 0 0 0
’ 12
2 1 1 5 1
1 1 1
1 1 0 2 6 12
2 2 2
, , 1 2 1
1 6 3 6
1 1 1 1
2 2 2 2 1 2 1
6 3 6

which have order 2, 2 and 3, respectively. The ¬rst array represents the
Crank-Nicolson method.

As for semi-implicit RK methods, we limit ourselves to mentioning the
case of DIRK methods (diagonally implicit RK), which, for s = 3, are
represented by the following Butcher array
1+µ 1+µ
0 0
2 2
’µ 1+µ
1
0
2 2 2
1 + µ ’1 ’ 2µ
1’µ 1+µ
2 2


1’
1 1 1
6µ2 3µ2 6µ2

The parameter µ represents one of the three roots of 3µ3 ’ 3µ ’ 1 = 0 (i.e.,
√ √ √
(2/ 3) cos(10o ), ’(2/ 3) cos(50o ), ’(2/ 3) cos(70o )). The maximum or-
der that has been determined in the literature for these methods is 4.


11.8.4 Regions of Absolute Stability for RK Methods
Applying an s-stage RK method to the model problem (11.24) yields
s s
Ki = un + h» aij Kj , un+1 = un + h» bi Ki , (11.77)
i=1 i=1

that is, a ¬rst-order di¬erence equation. If K and 1 are the vectors of com-
ponents (K1 , . . . , Ks )T and (1, . . . , 1)T , respectively, then (11.77) becomes
un+1 = un + h»bT K,
K = un 1 + h»AK,
from which, K = (I ’ h»A)’1 1un and thus
un+1 = 1 + h»bT (I ’ h»A)’1 1 un = R(h»)un
where R(h») is the so-called stability function.
The RK method is absolutely stable, i.e., the sequence {un } satis¬es
(11.25), i¬ |R(h»)| < 1. Its region of absolute stability is given by
A = {z = h» ∈ C such that |R(h»)| < 1} .
11.9 Systems of ODEs 517

If the method is explicit, A is strictly lower triangular and the function R
can be written in the following form (see [DV84])

det(I ’ h»A + h»1bT )
R(h») = .
det(I ’ h»A)
Thus since det(I’h»A) = 1, R(h») is a polynomial function in the variable
h», |R(h»)| can never be less than 1 for all values of h». Consequently, A
can never be unbounded for an explicit RK method.
In the special case of an explicit RK of order s = 1, . . . , 4, one gets (see
[Lam91])
s
1
(h»)k .
R(h») =
k!
k=0

The corresponding regions of absolute stability are drawn in Figure 11.10.
Notice that, unlike MS methods, the regions of absolute stability of RK
methods increase in size as the order grows.

4

3.5

3 s=4

s=3
2.5

2
s=2
1.5

1
s=1
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
with y(t) decaying exponentially to zero as t ’ ∞. Conversely, if 0 <
stable √
β < 2 ± the eigenvalues are complex conjugate and damped oscillations
arise in the solution which again decays exponentially to zero as t ’ ∞.
Numerical approximations have been carried out using both the backward
Euler (BE) and Crank-Nicolson (CN) methods. We have set y(t) = 0
and used the following (physiological) values of the physical parameters:
L = 5 · 10’2 [m], R0 = 5 · 10’3 [m], ρw = 103 [Kgm’3 ], H = 3 · 10’4 [m] and
E = 9 · 105 [N m’2 ], from which γ 3.3[Kg ’1 m’2 ] and ± = 36 · 106 [s’2 ].
A sinusoidal function p ’ p0 = x∆p(a + b cos(ω0 t)) has been used to model
the pressure variation along the vessel direction x and time, where ∆p =
0.25 · 133.32 [N m’2 ], a = 10 · 133.32 [N m’2 ], b = 133.32 [N m’2 ] and the
pulsation ω0 = 2π/0.8 [rad s’1 ] corresponds to a heart beat.
The results reported below refer to the ring coordinate x = L/2. The

two (very di¬erent) cases (1) β = ± [s’1 ] and (2) β = ± [s’1 ] have been
analyzed; it is easily seen that in case (2) the sti¬ness quotient (see Section
11.10) is almost equal to ±, thus the problem is highly sti¬. We notice also
that in both cases the real parts of the eigenvalues of A are very large,
so that an appropriately small time step should be taken to accurately
describe the fast transient of the problem.
In case (1) the di¬erential system has been studied on the time interval
[0, 2.5 · 10’3 ] with a time step h = 10’4 . We notice that the two eigenvalues
of A have modules equal to 6000, thus our choice of h is compatible with
the use of an explicit method as well.
Figure 11.13 (left) shows the numerical solutions as functions of time.
The solid (thin) line is the exact solution while the thick dashed and solid
526 11. Numerical Solution of Ordinary Di¬erential Equations

lines are the solutions given by the CN and BE methods, respectively. A far
better accuracy of the CN method over the BE is clearly demonstrated; this
is con¬rmed by the plot in Figure 11.13 (right) which shows the trajectories
of the computed solutions in the phase space. In this case the di¬erential
system has been integrated on the time interval [0, 0.25] with a time step
h = 2.5 · 10’4 . The dashed line is the trajectory of the CN method while
the solid line is the corresponding one obtained using the BE scheme. A
strong dissipation is clearly introduced by the BE method with respect to
the CN scheme; the plot also shows that both methods converge to a limit
cycle which corresponds to the cosine component of the forcing term.
x 10’5
14 0.6

12
0.5

10
0.4
8
0.3
6
0.2
4
0.1
2

0
0

’0.1
’2 0 0.5 1 1.5
0.5 1 1.5 2 2.5
0 ’3 ’4
x 10 x 10



FIGURE 11.13. Transient simulation (left) and phase space trajectories (right)


In the second case (2) the di¬erential system has been integrated on the
time interval [0, 10] with a time step h = 0.1. The sti¬ness of the problem is
demonstrated by the plot of the deformation velocities z shown in Figure
11.14 (left). The solid line is the solution computed by the BE method
while the dashed line is the corresponding one given by the CN scheme; for
the sake of graphical clarity, only one third of the nodal values have been
plotted for the CN method. Strong oscillations arise since the eigenvalues of
matrix A are »1 = ’1, »2 = ’36·106 so that the CN method approximates
the ¬rst component y of the solution y as
k
1 + (h»1 )/2
k ≥ 0,
CN
(0.9048)k ,
yk =
1 ’ (h»1 )/2
which is clearly stable, while the approximate second component z(= y ) is
k
1 + (h»2 )/2
k≥0
CN
(’0.9999)k ,
zk =
1 ’ (h»2 )/2
which is obviously oscillating. On the contrary, the BE method yields
k
1
k ≥ 0,
BE
(0.9090)k ,
yk =
1 ’ h»1
11.12 Exercises 527

and
k
1
k≥0
CN
(0.2777)k ,
zk =
1 ’ h»2
which are both stable for every h > 0. According to these conclusions the
¬rst component y of the vector solution y is correctly approximated by
both the methods as can be seen in Figure 11.14 (right) where the solid
line refers to the BE scheme while the dashed line is the solution computed
by the CN method.
’4
x 10
’4
1.2
x 10
2


1.5 1


1
0.8

0.5
0.6

0

0.4
’0.5


0.2
’1


0
’1.5
0 1 2 3 4 5 6 7 8 9 10
0 1 2 3 4 5 6 7 8 9 10



FIGURE 11.14. Long-time behavior of the solution: velocities (left) and displace-
ments (right)




11.12 Exercises
1. Prove that Heun™s method has order 2 with respect to h.
[Suggerimento : notice that h„n+1 = yn+1 ’ yn ’ h¦(tn , yn ; h) = E1 + E2 ,
tn+1
f (s, y(s))ds ’ h [f (tn , yn ) + f (tn+1 , yn+1 )]
where E1 = and E2 =
2
tn
{[f (tn+1 , yn+1 ) ’ f (tn+1 , yn + hf (tn , yn ))]}, where E1 is the error due to
h
2
numerical integration with the trapezoidal method and E2 can be bounded
by the error due to using the forward Euler method.]
2. Prove that the Crank-Nicoloson method has order 2 with respect to h.
[Solution : using (9.12) we get, for a suitable ξn in (tn , tn+1 )

h3
h
[f (tn , yn ) + f (tn+1 , yn+1 )] ’
yn+1 = yn + f (ξn , y(ξn ))
2 12

or, equivalently,

yn+1 ’ yn h2
1
= [f (tn , yn ) + f (tn+1 , yn+1 )] ’ f (ξn , y(ξn )). (11.90)
h 2 12

Therefore, relation (11.9) coincides with (11.90) up to an in¬nitesimal of
order 2 with respect to h, provided that f ∈ C 2 (I).]
528 11. Numerical Solution of Ordinary Di¬erential Equations

3. Solve the di¬erence equation un+4 ’ 6un+3 + 14un+2 ’ 16un+1 + 8un = n
subject to the initial conditions u0 = 1, u1 = 2, u2 = 3 and u3 = 4.
[Solution : un = 2n (n/4 ’ 1) + 2(n’2)/2 sin(π/4) + n + 2.]
4. Prove that if the characteristic polynomial Π de¬ned in (11.30) has simple
roots, then any solution of the associated di¬erence equation can be written
in the form (11.32).
[Hint : notice that a generic solution un+k is completely determined by
the initial values u0 , . . . , uk’1 . Moreover, if the roots ri of Π are distinct,
j j
there exist unique k coe¬cients ±i such that ±1 r1 + . . . + ±k rk = uj with
j = 0, . . . , k ’ 1 . . . ]
5. Prove that if the characteristic polynomial Π has simple roots, the matrix
R de¬ned in (11.37) is not singular.
[Hint: it coincides with the transpose of the Vandermonde matrix where
xj is replaced by rj (see Exercise 2, Chapter 8).]
i
i

6. The Legendre polynomials Li satisfy the di¬erence equation

(n + 1)Ln+1 (x) ’ (2n + 1)xLn (x) + nLn’1 (x) = 0

with L0 (x) = 1 and L1 (x) = x (see Section 10.1.2). De¬ning the generating
function F (z, x) = ∞ Pn (x)z n , prove that F (z, x) = (1 ’ 2zx + z 2 )’1/2 .
n=0
7. Prove that the gamma function


e’t tz’1 dt, z ∈ C,
“(z) = Rez > 0
0


is the solution of the di¬erence equation “(z + 1) = z“(z)
[Hint : integrate by parts.]
8. Study, as functions of ± ∈ R, stability and order of the family of linear
multistep methods

un+1 = ±un + (1 ’ ±)un’1 + 2hfn + [fn’1 ’ 3fn ] .
2
9. Consider the following family of linear multistep methods depending on
the real parameter ±
± ±
un+1 = un + h[(1 ’ )f (xn , un ) + f (xn+1 , un+1 )].
2 2
Study their consistency as a function of ±; then, take ± = 1 and use the
corresponding method to solve the Cauchy problem

y (x) = ’10y(x), x > 0,
y(0) = 1.

Determine the values of h in correspondance of which the method is abso-
lutely stable.
[Solution : the only consistent method of the family is the Crank-Nicolson
method (± = 1).]
11.12 Exercises 529

10. Consider the family of linear multistep methods
h
(2(1 ’ ±)fn+1 + 3±fn ’ ±fn’1 )
un+1 = ±un +
2
where ± is a real parameter.
(a) Analyze consistency and order of the methods as functions of ±, de-
termining the value ±— for which the resulting method has maximal
order.
(b) Study the zero-stability of the method with ± = ±— , write its charac-
teristic polynomial Π(r; h») and, using MATLAB, draw its region of
absolute stability in the complex plane.
11. Adams methods can be easily generalized, integrating between tn’r and
tn+1 with r ≥ 1. Show that, by doing so, we get methods of the form
p
un+1 = un’r + h bj fn’j
j=’1


and prove that for r = 1 the midpoint method introduced in (11.43) is
recovered (the methods of this family are called Nystron methods.)
12. Check that Heun™s method (11.10) is an explicit two-stage RK method and
write the Butcher arrays of the method. Then, do the same for the modi¬ed
Euler method, given by
h h
n ≥ 0. (11.91)
un+1 = un + hf (tn + , un + fn ),
2 2

[Solution : the methods have the following Butcher arrays

0 0 0
0 0 0
1 1
0 .]
1 1 0 22 2
2
3
3
2 12 1
0 1
2 2


13. Check that the Butcher array for method (11.73) is given by

0 0 0 0 0
2 12
1
0 0 0
2 2
2 12
1
0 0 0
2 2
1 0 0 1 0
2 12
1 1 1
6 3 3 6


14. Write a MATLAB program to draw the regions of absolute stability for a
RK method, for which the function R(h») is available. Check the code in
the special case of

R(h») = 1 + h» + (h»)2 /2 + (h»)3 /6 + (h»)4 /24 + (h»)5 /120 + (h»)6 /600

and verify that such a region is not connected.
530 11. Numerical Solution of Ordinary Di¬erential Equations

15. Determine the function R(h») associated with the Merson method, whose
Butcher array is

0 0 0 0 0 0
1 1
0 0 0 0
3 3
1 1 1
0 0 0
3 6 6
1 1 3
0 0 0
2 8 8
’2
1 3
1 0 2 0
2
1 2 1
0 0
6 3 6

4 i
+ (h»)5 /144.]
[Solution : one gets R(h») = 1 + i=1 (h») /i!
12
Two-Point Boundary Value Problems




This chapter is devoted to the analysis of approximation methods for two-
point boundary value problems for di¬erential equations of elliptic type.
Finite di¬erences, ¬nite elements and spectral methods will be considered.
A short account is also given on the extension to elliptic boundary value
problems in two-dimensional regions.


12.1 A Model Problem
To start with, let us consider the two-point boundary value problem
’u (x) = f (x), 0 < x < 1, (12.1)

u(0) = u(1) = 0. (12.2)
From the fundamental theorem of calculus, if u ∈ C 2 ([0, 1]) and satis¬es
the di¬erential equation (12.1) then
x

u(x) = c1 + c2 x ’ F (s) ds
0
s
where c1 and c2 are arbitrary constants and F (s) = f (t) dt. Using
0
integration by parts one has
x x x

F (s) ds = [sF (s)]x ’ (x ’ s)f (s) ds.
sF (s) ds =
0
0 0 0
532 12. Two-Point Boundary Value Problems

The constants c1 and c2 can be determined by enforcing the boundary
conditions. The condition u(0) = 0 implies that c1 = 0, and then u(1) = 0
1
yields c2 = 0 (1 ’ s)f (s) ds. Consequently, the solution of (12.1)-(12.2)
can be written in the following form
1 x

(1 ’ s)f (s) ds ’ (x ’ s)f (s) ds
u(x) = x
0 0

or, more compactly,
1

u(x) = G(x, s)f (s) ds, (12.3)
0

where, for any ¬xed x, we have de¬ned

s(1 ’ x) if 0 ¤ s ¤ x,
G(x, s) = (12.4)
x(1 ’ s) if x ¤ s ¤ 1.

The function G is called Green™s function for the boundary value problem
(12.1)-(12.2). It is a piecewise linear function of x for ¬xed s, and vice versa.
It is continuous, symmetric (i.e., G(x, s) = G(s, x) for all x, s ∈ [0, 1]), non
1
negative, null if x or s are equal to 0 or 1, and 0 G(x, s) ds = 1 x(1 ’ x).
2
The function is plotted in Figure 12.1.
We can therefore conclude that for every f ∈ C 0 ([0, 1]) there is a unique
solution u ∈ C 2 ([0, 1]) of the boundary value problem (12.1)-(12.2) which
admits the representation (12.3). Further smoothness of u can be derived
by (12.1); indeed, if f ∈ C m ([0, 1]) for some m ≥ 0 then u ∈ C m+2 ([0, 1]).
0.25




0.2




0.15




0.1




0.05




0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1



FIGURE 12.1. Green™s function for three di¬erent values of x: x = 1/4 (solid
line), x = 1/2 (dashed line), x = 3/4 (dash-dotted line)
12.2 Finite Di¬erence Approximation 533

An interesting property of the solution u is that if f ∈ C 0 ([0, 1]) is a
nonnegative function, then u is also nonnegative. This is referred to as the
monotonicity property, and follows directly from (12.3), since G(x, s) ≥ 0
for all x, s ∈ [0, 1]. The next property is called the maximum principle and
states that if f ∈ C 0 (0, 1),
1
¤
u f (12.5)
∞ ∞
8
= max |u(x)| is the maximum norm. Indeed, since G is non-
where u ∞
0¤x¤1
negative,
1 1
1
|u(x)| ¤ G(x, s)|f (s)| ds ¤ f x(1 ’ x) f
G(x, s) ds =
∞ ∞
2
0 0

from which the inequality (12.5) follows.


12.2 Finite Di¬erence Approximation
n
We introduce on [0, 1] the grid points {xj }j=0 given by xj = jh where
n ≥ 2 is an integer and h = 1/n is the grid spacing. The approximation to
n
the solution u is a ¬nite sequence {uj }j=0 de¬ned only at the grid points
(with the understanding that uj approximates u(xj )) by requiring that
uj+1 ’ 2uj + uj’1
’ for j = 1, . . . , n ’ 1
= f (xj ), (12.6)
h2
and u0 = un = 0. This corresponds to having replaced u (xj ) by its second
order centred ¬nite di¬erence (10.65) (see Section 10.10.1).
If we set u = (u1 , . . . , un’1 )T and f = (f1 , . . . , fn’1 )T , with fi = f (xi ),
it is a simple matter to see that (12.6) can be written in the more compact
form

Afd u = f , (12.7)

where Afd is the symmetric (n ’ 1) — (n ’ 1) ¬nite di¬erence matrix de¬ned
as

Afd = h’2 tridiagn’1 (’1, 2, ’1). (12.8)

This matrix is diagonally dominant by rows; moreover, it is positive de¬nite
since for any vector x ∈ Rn’1
n’1
’2
(xi ’ xi’1 )2 .
T
x2 x2
x Afd x = h + +
1 n’1
i=2
534 12. Two-Point Boundary Value Problems

This implies that (12.7) admits a unique solution. Another interesting prop-
erty is that Afd is an M-matrix (see De¬nition 1.25 and Exercise 2), which
guarantees that the ¬nite di¬erence solution enjoys the same monotonic-
ity property as the exact solution u(x), namely u is nonnegative if f is
nonnegative. This property is called discrete maximum principle.

In order to rewrite (12.6) in operator form, let Vh be a collection of discrete
functions de¬ned at the grid points xj for j = 0, . . . , n. If vh ∈ Vh , then
vh (xj ) is de¬ned for all j and we sometimes use the shorthand notation vj
0
instead of vh (xj ). Next, we let Vh be the subset of Vh containing discrete
functions that are zero at the endpoints x0 and xn . For a function wh we
de¬ne the operator Lh by
wj+1 ’ 2wj + wj’1
(Lh w)(xj ) = ’ j = 1, . . . , n ’ 1
, (12.9)
h2
and reformulate the ¬nite di¬erence problem (12.6) equivalently as: ¬nd
uh ∈ Vh such that
0


for j = 1, . . . , n ’ 1.
(Lh uh )(xj ) = f (xj ) (12.10)
Notice that, in this formulation, the boundary conditions are taken care of
by the requirement that uh ∈ Vh .0

Finite di¬erences can be used to provide approximations of higher-order
di¬erential operators than the one considered in this section. An example
is given in Section 4.7.2 where the ¬nite di¬erence centred discretization of
the fourth-order derivative ’u(iv) (x) is carried out by applying twice the
discrete operator Lh (see also Exercise 11). Again, extra care is needed to
properly handle the boundary conditions.


12.2.1 Stability Analysis by the Energy Method
For two discrete functions wh , vh ∈ Vh we de¬ne the discrete inner product
n
(wh , vh )h = h ck wk vk ,
k=0

with c0 = cn = 1/2 and ck = 1 for k = 1, . . . , n ’ 1. This is nothing but the
composite trapezoidal rule (9.13) which is here used to evaluate the inner
1
product (w, v) = 0 w(x)v(x)dx. Clearly,
1/2
vh = (vh , vh )h
h

is a norm on Vh .

Lemma 12.1 The operator Lh is symmetric, i.e.
∀ wh , vh ∈ Vh ,
0
(Lh wh , vh )h = (wh , Lh vh )h
12.2 Finite Di¬erence Approximation 535

and is positive de¬nite, i.e.

(Lh vh , vh )h ≥ 0 ∀vh ∈ Vh ,
0


with equality only if vh ≡ 0.
Proof. From the identity

wj+1 vj+1 ’ wj vj = (wj+1 ’ wj )vj + (vj+1 ’ vj )wj+1 ,

upon summation over j from 0 to n ’ 1 we obtain the following relation for all
wh , vh ∈ Vh
n’1 n’1
(wj+1 ’ wj )vj = wn vn ’ w0 v0 ’ (vj+1 ’ vj )wj+1
j=0 j=0


which is referred to as summation by parts. Using summation by parts twice, and
setting for ease of notation w’1 = v’1 = 0, for all wh , vh ∈ Vh we obtain
0


n’1
’1
= ’h [(wj+1 ’ wj ) ’ (wj ’ wj’1 )] vj
(Lh wh , vh )h
j=0
n’1
= h’1 (wj+1 ’ wj )(vj+1 ’ vj ).
j=0


From this relation we deduce that (Lh wh , vh )h = (wh , Lh vh )h ; moreover, taking
wh = vh we obtain
n’1
’1
(vj+1 ’ vj )2 .
(Lh vh , vh )h = h (12.11)
j=0


This quantity is always positive, unless vj+1 = vj for j = 0, . . . , n ’ 1, in which
3
case vj = 0 for j = 0, . . . , n since v0 = 0.

For any grid function vh ∈ Vh we de¬ne the following norm
0

± 1/2
 n’1 2
vj+1 ’ vj
|||vh |||h = h . (12.12)
 
h
j=0


Thus, (12.11) is equivalent to

(Lh vh , vh )h = |||vh |||2 for all vh ∈ Vh .
0
(12.13)
h


Lemma 12.2 The following inequality holds for any function vh ∈ Vh
0


1
¤ √ |||vh |||h .
vh (12.14)
h
2
536 12. Two-Point Boundary Value Problems

Proof. Since v0 = 0, we have
j’1
vk+1 ’ vk
for all j = 1, . . . , n ’ 1.
vj = h
h
k=0

Then,
2
j’1
vk+1 ’ vk
2 2
vj =h .
h
k=0

Using the Minkowski inequality
2
m m
¤m p2
pk (12.15)
k
k=1 k=1

which holds for every integer m ≥ 1 and every sequence {p1 , . . . , pm } of real
numbers (see Exercise 4), we obtain
n’1 j’1
n’1
vk+1 ’ vk 2
¤h
2 2
vj j .
h
j=1 j=1 k=0


Then for every vh ∈ Vh we get
0


n’1 n’1 n’1
(n ’ 1)n
vk+1 ’ vk 2
¤h |||vh |||2 .
vh 2 2 2
= h2
=h vj jh
h h
h 2
j=1 j=1 k=0

3
Inequality (12.14) follows since h = 1/n.

(1)
Remark 12.1 For every vh ∈ Vh , the grid function vh whose grid values
0

are (vj+1 ’ vj )/h, j = 0, . . . , n ’ 1, can be regarded as a discrete derivative
of vh (see Section 10.10.1). Inequality (12.14) can thus be rewritten as
1 (1)
¤ √ vh ∀vh ∈ Vh .
0
vh h h
2
It can be regarded as the discrete counterpart in [0, 1] of the following
Poincar´ inequality: for every interval [a, b] there exists a constant CP > 0
e
such that

¤ CP v (1)
v (12.16)
L2 (a,b) L2 (a,b)

for all v ∈ C 1 ([a, b]) such that v(a) = v(b) = 0 and where · is the
L2 (a,b)
norm in L2 (a, b) (see (8.25)).


Inequality (12.14) has an interesting consequence. If we multiply every
equation of (12.10) by uj and then sum for j from 0 on n ’ 1, we obtain

(Lh uh , uh )h = (f, uh )h .
12.2 Finite Di¬erence Approximation 537

Applying to (12.13) the Cauchy-Schwarz inequality (1.14) (valid in the
¬nite dimensional case), we obtain

|||uh |||2 ¤ fh uh
h h
h

where fh ∈ Vh is the grid function such that fh (xj ) = f (xj ) for all j =
1, . . . , n.
Owing to (12.14) we conclude that
1
¤
uh fh (12.17)
h h
2
from which we deduce that the ¬nite di¬erence problem (12.6) has a unique
solution (equivalently, the only solution corresponding to fh = 0 is uh = 0).
Moreover, (12.17) is a stability result, as it states that the ¬nite di¬erence
solution is bounded by the given datum fh .
To prove convergence, we ¬rst introduce the notion of consistency. Ac-
cording to our general de¬nition (2.13), if f ∈ C 0 ([0, 1]) and u ∈ C 2 ([0, 1])
is the corresponding solution of (12.1)-(12.2), the local truncation error is
the grid function „h de¬ned by

„h (xj ) = (Lh u)(xj ) ’ f (xj ), j = 1, . . . , n ’ 1. (12.18)

<<

. 21
( 26)



>>