ńņš. 21 |

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

hĪ±

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 |