<< . .

. 86
( : 95)



. . >>

(BE) (θ = 1, solid line) and of the Crank-Nicolson scheme (CN) (θ = 1/2, dashed
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,

<< . .

. 86
( : 95)



. . >>