<< . .

. 91
( : 95)



. . >>


[u]j = u+ (xj ) ’ u’ (xj ), u± (xj ) = lim u(xj + sa), j = 1, . . . , n ’ 1.
± s’0


Then, for any t ∈ [0, T ] the stability estimate for problem (13.68) reads
« 
t n’1
 uh („ ) a(xj )[uh („ )]2  d„
2
2
uh (t) + +
L2 (±,β)
L2 (±,β) j
j=0
® 
0
(13.69)
t

¤ C ° u0,h + a•2 („ ) d„ » .
2
2
+ f („ ) L2 (±,β)
L2 (±,β)
0


As for the convergence analysis, the following error estimate can be proved
for continuous ¬nite elements of degree r, r ≥ 1 (see [QV94], Section 14.3.1)
« 1/2
t

+  a|u(±, „ ) ’ uh (±, „ )|2 d„ 
max u(t) ’ uh (t) L2 (±,β)
t∈[0,T ]
0
= O( u0 ’ u0,h + hr ),
L2 (±,β)


If, instead, discontinuous ¬nite elements of degree r are used, r ≥ 0, the
convergence estimate becomes (see [QV94], Section 14.3.3 and the refer-
ences therein)
« T

+
max u(t) ’ uh (t) u(t) ’ uh (t) 2
dt
L2 (±,β) L2 (±,β)
t∈[0,T ]
0
1/2
T n’1

a(xj ) [u(t) ’ uh (t)]2 dt = O( u0 ’ u0,h + hr+1/2 ).
+ L2 (±,β)
j
0 j=0



13.10.2 Time Discretization
The time discretization of the ¬nite element schemes introduced in the
previous section can be carried out by resorting either to ¬nite di¬erences
or ¬nite elements. If an implicit ¬nite di¬erence scheme is adopted, both
method (13.65) and (13.68) are unconditionally stable.
As an example, let us use the backward Euler method for the time dis-
cretization of problem (13.65). We obtain for each n ≥ 0: ¬nd un+1 ∈ Vh
h
13.10 Finite Element Approximation of Hyperbolic Equations 621

such that
β β
‚un+1
1
(un+1 ’ un )vh dx + h
a vh dx
h
h
∆t ‚x
± ± (13.70)
β β

∀vh ∈ Vh
an+1 un+1 vh dx = f n+1 vh dx in
+ 0 h
± ±

with un+1 (±) = •n+1 and u0 = u0h . If f ≡ 0 and • ≡ 0, taking vh = un+1
h
h h
in (13.70) we can obtain
1
’ un ¤0
un+1 + a(β)(un+1 (β))2 + µ0 un+1
2 2 2
L2 (±,β) L2 (±,β) L2 (±,β)
h
h h h
2∆t
∀n ≥ 0. Summing for n from 0 to m ’ 1 yields for m ≥ 1
« 
m m
+ 2∆t  a(β)(uj+1 (β))2  ¤ u0
uj 2 2 (±,β)
um 2 2 (±,β) 2
+ L2 (±,β) .
hL h
hL h
j=1 j=1


In particular, we conclude that

¤ u0 ∀m ≥ 0.
um L2 (±,β) L2 (±,β)
h h



On the other hand, explicit schemes for the hyperbolic equations are
subject to a stability condition: for example, in the case of the forward
Euler method the stability condition is ∆t = O(∆x). In practice, this
restriction is not as severe as happens in the case of parabolic equations
and for this reason explicit schemes are often used in the approximation of
hyperbolic equations.

Programs 102 and 103 provide an implementation of the discontinous
Galerkin-¬nite element method of degree 0 (dG(0)) and 1 (dG(1)) in space
coupled with the backward Euler method in time for the solution of (13.26)
on the space-time domain (±, β) — (t0 , T ).

Program 102 - ipeidg0 : dG(0) implicit Euler
function [u,x]=ipeidg0(I,n,a,u0,bc)
nx = n(1); h = (I(2)-I(1))/nx;
x = [I(1)+h/2:h:I(2)]; t = I(3); u = (eval(u0))™;
nt = n(2); k = (I(4)-I(3))/nt;
lambda = k/h; e = ones(nx,1);
A=spdiags([-a*lambda*e, (1+a*lambda)*e],-1:0,nx,nx);
[L,U]=lu(A);
for t = I(3)+k:k:I(4)
622 13. Parabolic and Hyperbolic Initial Boundary Value Problems

f = u;
if a > 0
f(1) = a*bc(1)+f(1);
elseif a <= 0
f(nx) = a*bc(2)+f(nx);
end
y = L \ f; u = U \ y;
end




Program 103 - ipeidg1 : dG(1) implicit Euler



function [u,x]=ipeidg1(I,n,a,u0,bc)
nx = n(1); h = (I(2)-I(1))/nx;
x = [I(1):h:I(2)]; t = I(3); um = (eval(u0))™;
u = []; xx=[];
for i = 1:nx+1
u = [u, um(i), um(i)]; xx = [xx, x(i), x(i)];
end
u = u™; nt = n(2); k = (I(4)-I(3))/nt;
lambda = k/h; e = ones(2*nx+2,1);
B = spdiags([1/6*e,1/3*e,1/6*e],-1:1,2*nx+2,2*nx+2);
dd = 1/3+0.5*a*lambda; du = 1/6+0.5*a*lambda;
dl = 1/6-0.5*a*lambda; A=sparse([]);
A(1,1) = dd; A(1,2) = du; A(2,1) = dl; A(2,2) = dd;
for i=3:2:2*nx+2
A(i,i-1) =-a*lambda;
A(i,i) = dd; A(i,i+1) = du;
A(i+1,i) = dl; A(i+1,i+1) = A(i,i);
end
[L,U]=lu(A);
for t = I(3)+k:k:I(4)
f = B*u;
if a > 0
f(1) = a*bc(1)+f(1);
elseif a <= 0
f(nx) = a*bc(2)+f(nx);
end
y = L \ f;
u = U \ y;
end
x = xx;
13.11 Applications 623

13.11 Applications
13.11.1 Heat Conduction in a Bar
Consider a homogeneous bar of unit length with thermal conductivity ν,
which is connected at the endpoints to an external thermal source at a ¬xed
temperature, say u = 0. Let u0 (x) be the temperature distribution along
the bar at time t = 0 and f = f (x, t) be a given heat production term.
Then, the initial-boundary value problem (13.1)-(13.4) provides a model of
the time evolution of the temperature u = u(x, t) throughout the bar.
In the following, we study the case where f ≡ 0 and the temperature of
the bar is suddenly raised at the points around 1/2. A rough mathematical
model for this situation is provided, for instance, by taking u0 = K in
a certain subinterval [a, b] ⊆ [0, 1] and equal to 0 outside, where K is a
given positive constant. The initial condition is therefore a discontinuous
function.
We have used the θ-method with θ = 0.5 (Crank-Nicolson method, CN)
and θ = 1 (Backward Euler method, BE). Program 100 has been run
with h = 1/20, ∆t = 1/40 and the obtained solutions at time t = 2 are
shown in Figure 13.10. The results show that the CN method su¬ers a clear
instability due to the low smoothness of the initial datum (about this point,
see also [QV94], Chapter 11). On the contrary, the BE method provides a
stable solution which decays correctly to zero as t grows since the source
term f is null.


1
0.8
0.9

0.8
0.6
0.7

0.6
0.4
0.5

0.4
0.2
0.3

0.2
0
0.1

’0.2 0
0
0
0.5
0.5

1
1
0 0
0.2
1.5 1.5 0.4
0.5 0.6
0.8
2 1 2 1




FIGURE 13.10. Solutions for a parabolic problem with discontinuous initial da-
tum: CN method (left) and BE method (right)



13.11.2 A Hyperbolic Model for Blood Flow Interaction with
Arterial Walls
Let us consider again the problem of the ¬‚uid-structure interaction in a
cylindrical artery considered in Section 11.11.2, where the simple indepen-
dent rings model (11.88) was adopted.
624 13. Parabolic and Hyperbolic Initial Boundary Value Problems

If the axial action due to the tension between the di¬erent rings is
no longer neglected, denoting by z the longitudinal coordinate, equation
(11.88) modi¬es into

‚2· ‚ 2 · HE
ρw H 2 ’ σz 2 + 2 · = P ’ P0 , t > 0, 0<z<L (13.71)
‚t ‚z R0

where σz is the radial component of the axial stress and L is the length of
the cylindrical arterial district which is considered. In particular, neglecting
the third term on the left hand side and letting γ 2 = σz /(ρw H), f =
(P ’ P0 )/(ρw H), we recover the wave equation (13.33).

We have performed two sets of numerical experiments using the Leap-
Frog (LF) and Newmark (NW) methods. In the ¬rst example the space-
time domain of integration is the cylinder (0, 1) — (0, 1) and the source
term is f = (1 + π 2 γ 2 )e’t sin(πx) in such a way that the exact solution is
u(x, t) = e’t sin(πx). Table 13.3 shows the estimated orders of convergence
of the two methods, denoted by pLF and pN W respectively.
To compute these quantities we have ¬rst solved the wave equation on
(k)
four grids with sizes ∆x = ∆t = 1/(2k · 10), k = 0, . . . , 3. Let uh denote
the numerical solution corresponding to the space-time grid at the k-th
(0)
re¬ning level. Moreover, for j = 1, . . . , 10, let tj = j/10 be the time
discretization nodes of the grid at the coarsest level k = 0. Then, for each
level k, the maximum nodal errors ek on the k-th spatial grid have been
j
(0) (k)
evaluated at each time tj in such a way that the convergence order pj
can be estimated as

log(e0 /ek )
j j
(k)
pj = , k = 1, 2, 3.
k)
log(2

The results clearly show second-order convergence for both the methods,
as theoretically expected.

In the second example we have taken the following expressions for the
coe¬cient and source term: γ 2 = σz /(ρw H), with σz = 1 [Kgs’2 ], f =
(x∆p · sin(ω0 t))/(ρw H). The parameters ρw , H and the length L of the
vessel are as in Section 11.11.2. The space-time computational domain is
(0, L) — (0, T ), with T = 1 [s].
The Newmark method has been ¬rst used with ∆x = L/10 and ∆t =
T /100; the corresponding value of γ» is 3.6515, where » = ∆t/∆x. Since the
Newmark method is unconditionally stable, no spurious oscillations arise as
is con¬rmed by Figure 13.11, left. Notice the correct periodical behaviour
of the solution with a period corresponding to one heart beat; notice also
that with the present values of ∆t and ∆x the Leap-Frog method cannot
be employed since the CFL condition is not satis¬ed. To overcome this
13.12 Exercises 625

(0) (1) (2) (3) (0) (1) (2) (3)
tj pLF pLF pLF tj pN W pN W pN W
0.1 2.0344 2.0215 2.0151 0.1 1.9549 1.9718 1.9803
0.2 2.0223 2.0139 2.0097 0.2 1.9701 1.9813 1.9869
0.3 2.0170 2.0106 2.0074 0.3 1.9754 1.9846 1.9892
0.4 2.0139 2.0087 2.0061 0.4 1.9791 1.9869 1.9909
0.5 2.0117 2.0073 2.0051 0.5 1.9827 1.9892 1.9924
0.6 2.0101 2.0063 2.0044 0.6 1.9865 1.9916 1.9941
0.7 2.0086 2.0054 2.0038 0.7 1.9910 1.9944 1.9961
0.8 2.0073 2.0046 2.0032 0.8 1.9965 1.9979 1.9985
0.9 2.0059 2.0037 2.0026 0.9 2.0034 2.0022 2.0015
1.0 2.0044 2.0028 2.0019 1.0 2.0125 2.0079 2.0055
TABLE 13.3. Estimated orders of convergence for the Leap-Frog (LF) and New-
mark (NW) methods


problem, we have therefore chosen a much smaller time-step ∆t = T /400,
in such a way that γ» 0.9129 and the Leap-Frog scheme can be applied.
The obtained result is shown in Figure 13.11, right; a similar solution has
been computed by using the Newmark method with the same values of the
discretization parameters.

x 10’4 x 10’4
3
3

2
2

1
1

0
0

’1
’1

’2
’2
0.06
’3
’3 0.04
0.04
1
0.8 0.03 1 0.8
0.6 0.02
0.02 0.6
0.4 0.4
0.01
0.2 0.2 0
0 0
0



FIGURE 13.11. Computed solutions using the NM method on a space-time grid
with ∆t = T /100 and ∆x = L/10 (left) and the LF scheme on a space-time grid
with the same value of ∆x but with ∆t = T /400 (right)




13.12 Exercises
1. Apply the θ-method (13.9) to the approximate solution of the scalar Cauchy
problem (11.1) and using the analysis of Section 11.3 prove that the local
truncation error is of the order of ∆t + h2 if θ = 1 while it is of the order
2
of ∆t2 + h2 if θ = 1 .
2

2. Prove that in the case of piecewise linear ¬nite elements, the mass-lumping
process described in Section 13.3 is equivalent to computing the integrals
626 13. Parabolic and Hyperbolic Initial Boundary Value Problems

1
mij = 0 •j •i dx by the trapezoidal quadrature formula (9.11). This, in
particular, shows that the diagonal matrix M is nonsingular.
[Hint: ¬rst, verify that exact integration yields
±
 1 i = j,
h 2
mij =
6
1 i = j.

Then, apply the trapezoidal rule to compute mij recalling that •i (xj ) =
δij .]

<< . .

. 91
( : 95)



. . >>