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