<< . .

. 76
( : 95)



. . >>

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

<< . .

. 76
( : 95)



. . >>