j=1 j=1

0

Using the same notation as in (13.8) we obtain

.

Mu(t) + Afe u(t) = ffe (t) (13.14)

1

where Afe = (a(•j , •i )), ffe (t) = (F (•i )) and M = (mij ) = ( 0 •j •i dx) for

i, j = 1, . . . , Nh . M is called the mass matrix. Since it is nonsingular, the

system of ODEs (13.14) can be written in normal form as

.

u(t) = ’M’1 Afe u(t) + M’1 ffe (t). (13.15)

To solve (13.15) approximately we can still apply the θ-method and obtain

uk+1 ’ uk

+ Afe θuk+1 + (1 ’ θ)uk = θffe + (1 ’ θ)ffe .

k+1 k

M (13.16)

∆t

As usual, the upper index k means that the quantity at hand is computed at

time tk . As in the ¬nite di¬erence case, for θ = 0, 1 and 1/2 we respectively

obtain the forward Euler, backward Euler and Crank-Nicolson methods,

where the Crank-Nicolson method is the only one which is second-order

accurate with respect to ∆t.

For each k, (13.16) is a linear system whose matrix is

1

K= M + θAfe .

∆t

Since M and Afe are symmetric and positive de¬nite, the matrix K is also

symmetric and positive de¬nite. Thus, its Cholesky decomposition K =

HT H where H is upper triangular (see Section 3.4.2) can be carried out at

t = 0. Consequently, at each time step the following two linear triangular

systems, each of size equal to Nh , must be solved, with a computational

2

cost of Nh /2 ¬‚ops

±

HT y = 1 M ’ (1 ’ θ)A uk + θf k+1 + (1 ’ θ)f k ,

fe fe

fe

∆t

Huk+1 = y.

588 13. Parabolic and Hyperbolic Initial Boundary Value Problems

When θ = 0, a suitable diagonalization of M would allow to decouple

the system equations (13.16). The procedure is carried out by the so-called

mass-lumping in which we approximate M by a nonsingular diagonal matrix

M. In the case of piecewise linear ¬nite elements M can be obtained using

the composite trapezoidal formula over the nodes {xi } to evaluate the

1

integrals 0 •j •i dx, obtaining mij = hδij , i, j = 1, . . . , Nh (see Exercise

˜

2).

Stability Analysis of the θ-Method

13.3.1

Applying the θ-method to the Galerkin problem (13.13) yields

uk+1 ’ uk

+ a θuk+1 + (1 ’ θ)uk , vh

h

h

, vh h

h

∆t (13.17)

= θF k+1 (vh ) + (1 ’ θ)F k (vh ) ∀vh ∈ Vh

1

for k ≥ 0 and with u0 = u0 h , F k (vh ) = 0 f (tk )vh (x)dx. Since we are

h

interested in the stability analysis, we can consider the special case where

F = 0; moreover, for the time being, we focus on the case θ = 1 (implicit

Euler scheme), i.e.

uk+1 ’ uk

∀vh ∈ Vh .

+ a uk+1 , vh = 0

h

h

, vh h

∆t

Letting vh = uk+1 , we get

h

uk+1 ’ uk k+1

+ a(uk+1 , uk+1 ) = 0.

h

h

, uh h h

∆t

From the de¬nition of a(·, ·), it follows that

2

‚uk+1

uk+1 , uk+1 h

a =ν . (13.18)

h h

‚x

L2 (0,1)

Moreover, we remark that (see Exercise 3 for the proof of this result)

2

‚uk+1

¤ uk

uk+1 2 2 (0,1) 2

h

+ 2ν∆t L2 (0,1) . (13.19)

h

L

h

‚x

L2 (0,1)

It follows that, ∀n ≥ 1

2

n’1 n’1 n’1

‚uk+1

¤

uk+1 2 2 (0,1) uk 2

h

+ 2ν∆t L2 (0,1) .

h

L

h

‚x

L2 (0,1)

k=0 k=0 k=0

13.3 Finite Element Approximation of the Heat Equation 589

Since these are telescopic sums, we get

2

n’1

‚uk+1

¤ u0h

un 2 2 (0,1) 2

h

+ 2ν∆t L2 (0,1) , (13.20)

hL

‚x

L2 (0,1)

k=0

which shows that the scheme is unconditionally stable. Proceeding similarly

if f = 0, it can be shown that

2

n’1

‚uk+1

h

un 2

+ 2ν∆t

L2 (0,1)

h

‚x

L2 (0,1)

k=0 (13.21)

n

¤ C(n) ∆t f k 2

2

u0 h + ,

L2 (0,1)

L2 (0,1)

k=1

where C(n) is a constant independent of both h and ∆t.

Remark 13.1 The same kind of stability inequalities (13.20) and (13.21)

can be obtained if a(·, ·) is a more general bilinear form provided that it is

continuous and coercive (see Exercise 4).

To carry out the stability analysis of the θ-method for every θ ∈ [0, 1] we

need de¬ning the eigenvalues and eigenvectors of a bilinear form.

De¬nition 13.1 We say that » is an eigenvalue and w ∈ V is the associ-

ated eigenvector for the bilinear form a(·, ·) : V — V ’ R if

∀v ∈ V,

a(w, v) = »(w, v)

where (·, ·) denotes the usual scalar product in L2 (0, 1).

If the bilinear form a(·, ·) is symmetric and coercive, it has in¬nitely many

real positive eigenvalues that form an unbounded sequence; moreover, its

eigenvectors (called also eigenfunctions) form a basis for the space V .

At a discrete level the corresponding pair »h ∈ R, wh ∈ Vh satis¬es

∀vh ∈ Vh .

a(wh , vh ) = »h (wh , vh ) (13.22)

From the algebraic standpoint, problem (13.22) can be formulated as

Afe w = »h Mw

(where w is the vector of the gridvalues of wh ) and can be regarded

as a generalized eigenvalue problem (see Section 5.9). All the eigenvalues

»1 , . . . , »Nh are positive. The corresponding eigenvectors wh , . . . , wh h form

N

1

h h

590 13. Parabolic and Hyperbolic Initial Boundary Value Problems

a basis for the subspace Vh and can be chosen in such a way as to be or-

j

thonormal, i.e., such that (wh , wh ) = δij , ∀i, j = 1, . . . , Nh . In particular,

i

any function vh ∈ Vh can be represented as

Nh

j

vh (x) = vj wh (x).

j=1

Let us now assume that θ ∈ [0, 1] and focus on the case where the bilinear

form a(·, ·) is symmetric. Although the ¬nal stability result still holds in

the nonsymmetric case, the proof that follows cannot apply since in that

i

case the eigenvectors would no longer form a basis for Vh . Let wh be the

eigenvectors of a(·, ·) whose span forms an orthonormal basis for Vh . Since

at each time step uk ∈ Vh , we can express uk as

h h

Nh

j

uk (x) uk wh (x).

=

h j

j=1

i

Letting F = 0 in (13.17) and taking vh = wh , we ¬nd

Nh

1 j

uk+1 ’ uk i

wh , w h

j

j

∆t j=1

Nh

j

θuk+1 + (1 ’ θ)uk a(wh , wh ) = 0,

i

+ i = 1, . . . , Nh .

j

j

j=1

j

Since wh are eigenfunctions of a(·, ·) we obtain

a(wh , wh ) = »j (wh , wh ) = »j δij = »i ,

j j

i i

h

h h

so that

uk+1 ’ uk

+ θuk+1 + (1 ’ θ)uk »i = 0.

i

i

i h

i

∆t

Solving this equation with respect to uk+1 gives

i

1 ’ (1 ’ θ)»i ∆t

h

k+1

uk

ui = .

i i ∆t

1 + θ»h

In order for the method to be unconditionally stable we must have (see

Chapter 11)

1 ’ (1 ’ θ)»i ∆t

h

< 1,

i ∆t

1 + θ»h

that is

2

2θ ’ 1 > ’ .

»i ∆t

h

13.3 Finite Element Approximation of the Heat Equation 591

If θ ≥ 1/2, this inequality is satis¬ed for any value of ∆t. Conversely, if

θ < 1/2 we must have

2

∆t < .

(1 ’ 2θ)»i

h

Since this relation must hold for all the eigenvalues »i of the bilinear form,

h

it su¬ces requiring that it is satis¬ed for the largest of them, which we

assume to be »Nh . h

We therefore conclude that if θ ≥ 1/2 the θ-method is unconditionally

stable (i.e., it is stable ∀∆t), whereas if 0 ¤ θ < 1/2 the θ-method is stable

only if

2

∆t ¤ .

(1 ’ 2θ)»Nh

h

It can be shown that there exist two positive constants c1 and c2 , indepen-

dent of h, such that

c1 h’2 ¤ »Nh = c2 h’2

h

(see for the proof, [QV94], Section 6.3.2). Accounting for this, we obtain

that if 0 ¤ θ < 1/2 the method is stable only if

∆t ¤ C1 (θ)h2 , (13.23)

for a suitable constant C1 (θ) independent of both h and ∆t.

With an analogous proof, it can be shown that if a pseudo-spectral Galerkin

approximation is used for problem (13.12), the θ’method is uncondition-

ally stable if θ ≥ 1 , while for 0 ¤ θ < 1 stability holds only if

2 2

∆t ¤ C2 (θ)N ’4 , (13.24)

for a suitable constant C2 (θ) independent of both N and ∆t. The di¬erence

between (13.23) and (13.24) is due to the fact that the largest eigenvalue

of the spectral sti¬ness matrix grows like O(N 4 ) with respect to the degree

of the approximating polynomial.

Comparing the solution of the globally discretized problem (13.17) with

that of the semi-discrete problem (13.13), by a suitable use of the stability

result (13.21) and of the truncation time discretization error, the following

convergence result can be proved

u(tk ) ’ uk ¤ C(u0 , f, u)(∆tp(θ) + hr+1 ), ∀k ≥ 1

L2 (0,1)

h

where r denotes the piecewise polynomial degree of the ¬nite element space

Vh , p(θ) = 1 if θ = 1/2 while p(1/2) = 2 and C is a constant that depends

on its arguments (assuming that they are su¬ciently smooth) but not on

592 13. Parabolic and Hyperbolic Initial Boundary Value Problems

h and ∆t. In particular, if f ≡ 0 on can obtain the following improved

estimates

r+1 p(θ)

h ∆t

√

u(t ) ’ ¤C

k

uk L2 (0,1) + u0 L2 (0,1) ,

h

tk

tk

for k ≥ 1, θ = 1 or θ = 1/2. (For the proof of these results, see [QV94], pp.

394-395).

Program 100 provides an implementation of the θ-method for the solution

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

cretization in space is based on piecewise-linear ¬nite elements. The input

parameters are: the column vector I containing the endpoints of the space

interval (a = I(1), b = I(2)) and of the time interval (t0 = I(3), T = I(4));

the column vector n containing the number of steps in space and time; the

macros u0 and f containing the functions u0h and f , the constant viscosity

nu, the Dirichlet boundary conditions bc(1) and bc(2), and the value of

the parameter theta.

Program 100 - thetameth : θ-method for the heat equation

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

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

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

nt = n(2); k = (I(4)-I(3))/nt; e = ones(nx+1,1);

K = spdiags([(h/(6*k)-nu*theta/h)*e, (2*h/(3*k)+2*nu*theta/h)*e, ...

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

B = spdiags([(h/(6*k)+nu*(1-theta)/h)*e, (2*h/(3*k)-nu*2*(1-theta)/h)*e, ...

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

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

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

[L,U]=lu(K);

t = I(3); x=[I(1)+h:h:I(2)-h]; fold = (eval(f))™;

fold = h*fold;

fold = [bc(1); fold; bc(2)];

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

t = time; fnew = (eval(f))™;

fnew = h*fnew;

fnew = [bc(1); fnew; bc(2)];

b = theta*fnew+(1-theta)*fold + B*uold;

y = L \ b; u = U \ y;

uold = u;

end

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

Example 13.1 Let us assess the time-accuracy of the θ-method in the solution

of the heat equation (13.1) on the space-time domain (0, 1) — (0, 1) where f is

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

chosen in such a way that the exact solution is u = sin(2πx) cos(2πt). A ¬xed

spatial grid size h = 1/500 has been used while the time step ∆t is equal to

(10k)’1 , k = 1, . . . , 4. Finally, piecewise ¬nite elements are used for the space

discretization. Figure 13.1 shows the convergence behavior in the L2 (0, 1) norm

(evaluated at time t = 1), as ∆t tends to zero, of the backward Euler method