<< . .

. 85
( : 95)



. . >>

j
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

<< . .

. 85
( : 95)



. . >>