line). As expected, the CN method is far more accurate than the BE method. •

’1

10

’2

10

’3

10

’4

10

’5

10

’6

10

’7

10 1 2

10 10

FIGURE 13.1. Convergence analysis of the θ-method as a function of the number

1/∆t of time steps (represented on the x-axis): θ = 1 (solid line) and θ = 0.5

(dashed line)

13.4 Space-Time Finite Element Methods for the

Heat Equation

An alternative approach for time discretization is based on the use of a

Galerkin method to discretize both space and time variables.

Suppose to solve the heat equation for x ∈ [0, 1] and t ∈ [0, T ]. Let us

denote by Ik = [tk’1 , tk ] the k-th time interval for k = 1, . . . , n with ∆tk =

tk ’ tk’1 ; moreover, we let ∆t = maxk ∆tk ; the rectangle Sk = [0, 1] — Ik is

the so called space-time slab. At each time level tk , we consider a partition

Thk of (0, 1) into mk subintervals Kj = [xk , xk ], j = 0, . . . , mk ’ 1. We

k

j j+1

let hk = xk ’ xk and denote by hk = maxj hk and by h = maxk hk .

j j+1 j j

Let us now associate with Sk a space-time partition Sk = ∪mk Rj where k

j=1

Rj = Kj — Ik and Kj ∈ Thk . The space-time slab Sk is thus decomposed

k k k

k

into rectangles Rj (see Figure 13.2).

For each time slab Sk we introduce the space-time ¬nite element space

Qq (Sk ) = v ∈ C 0 (Sk ), v|Rj ∈ P1 (Kj ) — Pq (Ik ), j = 0, . . . , mk ’ 1

k

k

594 13. Parabolic and Hyperbolic Initial Boundary Value Problems

t

tk

j

Rk Sk

k’1

t

1x

0

FIGURE 13.2. Space-time ¬nite element discretization

where, usually, q = 0 or q = 1. Then, the space-time ¬nite element space

over [0, 1] — [0, T ] is de¬ned as follows

Vh,∆t = v : [0, 1] — [0, T ] ’ R : v|Sk ∈ Yh,k , k = 1, . . . , n ,

where

Yh,k = {v ∈ Qq (Sk ) : v(0, t) = v(1, t) = 0 ∀t ∈ Ik } .

The number of degrees of freedom of Vh,∆t is equal to (q + 1)(mk ’ 1).

The functions in Vh,∆t are linear and continuous in space while they are

piecewise polynomials of degree q in time. These functions are in general

discontinuous across the time levels tk and the partitions Th do not match

k

at the interface between contiguous time levels (see Figure 13.2). For this

reason, we adopt henceforth the following notation

v± = lim v(tk ± „ ), [v k ] = v+ ’ v’ .

k k k

„ ’0

The discretization of problem (13.12) using continuous ¬nite elements in

space of degree 1 and discontinuous ¬nite elements of degree q in time

(abbreviated by cG(1)dG(q) method) is: ¬nd U ∈ Vh,∆t such that

n n’1

‚U

([U k ], V+ )

k

,V + a(U, V ) dt +

‚t

k=1I k=1

k

T

0

∀V ∈ V h,∆t ,

0 0

+(U+ , V+ ) = (f, V ) dt,

0

where

0

V h,∆t = {v ∈ Vh,∆t : v(0, t) = v(1, t) = 0 ∀t ∈ [0, T ]} ,

13.4 Space-Time Finite Element Methods for the Heat Equation 595

1

U’ = u0h , U k = U (x, tk ) and (u, v) = 0 uv dx denotes the scalar product

0

of L2 (0, 1). The continuity of U at each point tk is therefore imposed only

in a weak sense.

To construct the algebraic equations for the unknown U we need expand-

ing it over a basis in time and in space. The single space-time basis func-

tion •k (x, t) can be written as •k (x, t) = •k (x)ψl (t), j = 1, . . . , mk ’ 1,

j

jl jl

k

l = 0, . . . , q, where •j is the usual piecewise linear basis function and ψl is

the l-th basis function of Pq (Ik ).

When q = 0 the solution U is piecewise constant in time. In that case

k

Nh

Uj •k (x), x ∈ [0, 1], t ∈ Ik ,

U k (x, t) = k

j

j=1

where Uj = U k (xj , t) ∀t ∈ Ik . Let

k

Ak = (aij ) = (a(•k , •k )), Mk = (mij ) = ((•k , •k )),

j i j i

«

fk = (fi ) = f (x, t)•k (x)dx dt , Bk,k’1 = (bij ) = ((•k , •k’1 )),

i j i

Sk

denote the sti¬ness matrix, the mass matrix, the data vector and the pro-

k’1

and Vh , respectively, at the time level tk .

k

jection matrix between Vh

Then, letting Uk = (Uj ), at each k-th time level the cG(1)dG(0) method

k

requires solving the following linear system

Mk + ∆tk Ak Uk = Bk,k’1 Uk’1 + fk ,

which is nothing else than the Euler backward discretization scheme with

a modi¬ed right hand side.

When q = 1, the solution is piecewise linear in time. For ease of notation

we let U k (x) = U’ (x, tk ) and U k’1 (x) = U+ (x, tk’1 ). Moreover, we assume

that the spatial partition Thk does not change with the time level and we

let mk = m for every k = 0, . . . , n. Then, we can write

tk ’ t t ’ tk’1

k’1 k

U|Sk = U (x) + U (x) .

∆tk ∆tk

Thus the cG(1)dG(1) method leads to the solution of the following 2 —

2 block-system in the unknowns Uk = (Uik ) and Uk’1 = (Uik’1 ), i =

1, . . . , m ’ 1

± k k

’ 1 M + ∆t A Uk’1 + 1 M + ∆t A Uk = f

k’1

k’1 + Bk,k’1 U’ ,

k k k k

2 3 2 6

1 ∆tk ∆tk

1

Ak Uk’1 + Ak Uk = fk

Mk + Mk +

2 6 2 3

596 13. Parabolic and Hyperbolic Initial Boundary Value Problems

where

f (x, t)•k (x)ψ1 (t)dx dt, fk =

k

f (x, t)•k (x)ψ2 (t)dx dt

k

fk’1 = i i

Sk Sk

and ψ1 (t) = (tk ’ t)/∆tk , ψ2 (t)(t ’ tk’1 )/∆tk are the two basis functions

k k

of P1 (Ik ).

Assuming that Vh,k’1 ‚ Vh,k , it is possible to prove that (see for the

proof [EEHJ96])

u(tn ) ’ U n ¤ C(u0h , f, u, n)(∆t2 + h2 ), (13.25)

L2 (0,1)

where C is a constant that depends on its arguments (assuming that they

are su¬ciently smooth) but not on h and ∆t.

An advantage in using space-time ¬nite elements is the possibility to

perform a space-time grid adaptivity on each time-slab based on a posteriori

error estimates (the interested reader is referred to [EEHJ96] where the

analysis of this method is carried out in detail).

Program 101 provides an implementation of the dG(1)cG(1) method for

the solution of the heat equation on the space-time domain (a, b) — (t0 , T ).

The input parameters are the same as in Program 100.

Program 101 - pardg1cg1 : dG(1)cG(1) method for the heat equation

function [u,x]=pardg1cg1(I,n,u0,f,nu,bc)

nx = n(1); h = (I(2)-I(1))/nx;

x = [I(1):h:I(2)]; t = I(3); um = (eval(u0))™;

nt = n(2); k = (I(4)-I(3))/nt;

e = ones(nx+1,1);

Add = spdiags([(h/12-k*nu/(3*h))*e, (h/3+2*k*nu/(3*h))*e, ...

(h/12-k*nu/(3*h))*e],-1:1,nx+1,nx+1);

Aud = spdiags([(h/12-k*nu/(6*h))*e, (h/3+k*nu/(3*h))*e, ...

(h/12-k*nu/(6*h))*e],-1:1,nx+1,nx+1);

Ald = spdiags([(-h/12-k*nu/(6*h))*e, (-h/3+k*nu/(3*h))*e, ...

(-h/12-k*nu/(6*h))*e],-1:1,nx+1,nx+1);

B = spdiags([h*e/6, 2*h*e/3, h*e/6],-1:1,nx+1,nx+1);

Add(1,1) = 1; Add(1,2) = 0; B(1,1) = 0; B(1,2) = 0;

Aud(1,1) = 0; Aud(1,2) = 0; Ald(1,1) = 0; Ald(1,2) = 0;

Add(nx+1,nx+1)=1; Add(nx+1,nx)=0;

B(nx+1,nx+1)=0; B(nx+1,nx) = 0;

Ald(nx+1,nx+1)=0; Ald(nx+1,nx)=0;

Aud(nx+1,nx+1)=0; Aud(nx+1,nx)=0;

[L,U]=lu([Add Aud; Ald Add]);

x = [I(1)+h:h:I(2)-h]; xx=[I(1),x,I(2)];

for time = I(3)+k:k:I(4)

13.5 Hyperbolic Equations: A Scalar Transport Problem 597

t = time; fq1 = 0.5*k*h*eval(f);

t = time-k; fq0 = 0.5*k*h*eval(f);

rhs0 = [bc(1), fq0, bc(2)]; rhs1 = [bc(1), fq1, bc(2)];

b = [rhs0™; rhs1™] + [B*um; zeros(nx+1,1)];

y = L \ b; u = U \ y; um = u(nx+2:2*nx+2,1);

end

x = [I(1):h:I(2)]; u = um;

Example 13.2 We assess the accuracy of the dG(1)cG(1) method on the same

test problem considered in Example 13.1. In order to neatly identify both spatial

and temporal contributions in the error estimate (13.25) we have performed the

numerical computations using Program 101 by varying either the time step or

the space discretization step only, having chosen in each case the discretization

step in the other variable su¬ciently small that the corresponding error can be

neglected. The convergence behavior in Figure 13.3 shows perfect agreement with

•

the theoretical results (second-order accuracy in both space and time).

’1

10

’2

10

’3

10

’4

10

’5

10 1 2 3

10 10 10

FIGURE 13.3. Convergence analysis for the dG(1)cG(1) method. The solid line

is the time discretization error while the dashed line is the space discretization

error. In the ¬rst case the x-axis denotes the number of time steps while in second

case it represents the number of space subintervals

13.5 Hyperbolic Equations: A Scalar Transport

Problem

Let us consider the following scalar hyperbolic problem

±

‚u + a ‚u = 0 x ∈ R, t > 0,

‚t ‚x (13.26)

u(x, 0) = u0 (x) x ∈ R,

598 13. Parabolic and Hyperbolic Initial Boundary Value Problems

where a is a positive real number. Its solution is given by

u(x, t) = u0 (x ’ at) t ≥ 0,

and represents a travelling wave with velocity a. The curves (x(t), t) in the

plane (x, t), that satisfy the following scalar ordinary di¬erential equation

±

dx(t) = a t > 0

dt (13.27)

x(0) = x0 ,

are called characteristic curves. They are the straight lines x(t) = x0 + at,

t > 0. The solution of (13.26) remains constant along them since

du ‚u ‚u dx

= + =0 on (x(t), t).

dt ‚t ‚x dt

For the more general problem

±

‚u + a ‚u + a u = f x ∈ R, t > 0,

0

‚t ‚x (13.28)

x ∈ R,

u(x, 0) = u0 (x)

where a, a0 and f are given functions of the variables (x, t), the charac-

teristic curves are still de¬ned as in (13.27). In this case, the solutions of

(13.28) satisfy along the characteristics the following di¬erential equation

du

= f ’ a0 u on (x(t), t).

dt

t

Q

t

t=1

¯

t P

P0 x x

0 1

± β

FIGURE 13.4. Left: examples of characteristics which are straight lines issuing

from the points P and Q. Right: characteristic straight lines for the Burgers

equation

Let us now consider problem (13.26) on a bounded interval. For example,

assume that x ∈ [±, β] and a > 0. Since u is constant along the character-

istics, from Figure 13.4 (left) we deduce that the value of the solution at P

13.6 Systems of Linear Hyperbolic Equations 599

attains the value of u0 at P0 , the foot of the characteristic issuing from P .

On the other hand, the characteristic issuing from Q intersects the straight

¯

line x(t) = ± at a certain time t = t > 0. Thus, the point x = ± is an

in¬‚ow point and it is necessary to assign a boundary value for u at x = ±

for every t > 0. Notice that if a < 0 then the in¬‚ow point is x = β.

Referring to problem (13.26) it is worth noting that if u0 is discontinuous

at a point x0 , then such a discontinuity propagates along the characteristics

issuing from x0 . This process can be rigorously formalized by introducing

the concept of weak solutions of hyperbolic problems, see e.g. [GR96]. An-

other reason for introducing weak solutions is that in the case of nonlinear

hyperbolic problems the characteristic lines can intersect: in this case the

solution cannot be continuous and no classical solution does exist.

Example 13.3 (Burgers equation) Let us consider the Burgers equation

‚u ‚u

x∈R

+u =0 (13.29)

‚t ‚x

which is perhaps the simplest nontrivial example of a nonlinear hyperbolic equa-

tion. Taking as initial condition

±

x ¤ 0,

1

1 ’ x 0 ¤ x ¤ 1,

u(x, 0) = u0 (x) =

x ≥ 1,

0

the characteristic line issuing from the point (x0 , 0) is given by

±

x0 ¤ 0,