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]).