<< . .

. 44
( : 59)



. . >>

+ b2 (s) — ] e’±(t’s) ds ,
+ [ (s) + b1 (s) —
0,h
0

where |bj (t, v)| ¤ bj (t) — v 0,h for all v ∈ Vh , t ∈ (0, T ), and j = 1, 2. The
¬rst term in the integral can be estimated by means of (7.65), whereas the
310 7. Discretization of Parabolic Problems

estimates of b1 (s) — , b2 (s) result from Lemma 7.14:

’±t
¤
θ(t) θ(0) 0,h e
0,h
t
e’±(t’s) ds .
+ |u (s)|1,∞ + f (s)
+Ch [ u (s) 1,∞ ]
2
0

If u0 ∈ V © H 2 („¦), we can write, by (7.65),
¤ uh0 ’ u0 + (I ’ Rh )u0 ¤ uh0 ’ u0
θ(0) + Ch u0 2 .
0,h 0,h 0,h 0,h

So we get
’±t
+ Ch u0 2 e’±t
¤ uh0 ’ u0
θ(t) 0,h e
0,h
t
e’±(t’s) ds .
+ |u (s)|1,∞ + f (s)
+ [ u (s) 1,∞ ]
2
0
Since
uh (t) ’ u(t) ¤ θ(t) + (Rh ’ I)u(t) 0,h ,
0,h 0,h

the obtained estimate and (7.65) yield the following result.
Theorem 7.15 In addition to the assumptions of Theorem 6.15, let f ∈
C([0, T ], C 1 („¦)), u0 ∈ V © H 2 („¦), and u0h ∈ Vh . Then if u(t) is su¬ciently
smooth, the solution uh (t) of the semidiscrete ¬nite volume method (7.48)
on Donald diagrams satis¬es the following estimate:

e’±t + Ch u0 2 e’±t + u(t)
uh (t) ’ u(t) ¤ uh0 ’ u0
0,h 0,h 2
t
e’±(t’s) ds .
+ |u (s)|1,∞ + f (s)
+ [ u (s) 1,∞ ]
2
0
Remark 7.16 In comparison with the ¬nite element method, the result is
not optimal in h. The reason is that, in general, the ¬nite volume method
does not yield optimal L2 -error estimates even in the elliptic case, but this
type of result is necessary to obtain optimal estimates.


Exercises
7.4 Let A ∈ RM,M be an RM-elliptic matrix and let the symmetric positive
de¬nite matrix B ∈ RM,M have the Cholesky decomposition B = E T E.
Show that the matrix A := E ’T AE ’1 is RM-elliptic.
ˆ

7.5 Prove identity (7.62) by ¬rst proving the corresponding identity for
one space dimension:
l’1
π a
sin2 ν ih = .
h
a 2
i=1
7.3. Fully Discrete Schemes 311



7.6 Let V be a Banach space and a : V — V ’ R a V -elliptic, continuous
bilinear form. Show that the Ritz projection Rh : V ’ Vh in a subspace
Vh ‚ V (cf. De¬nition 7.8) has the following properties:
(i) Rh : V ’ Vh is continuous because of Rh u V ¤ M
u ,
V
±
(ii) Rh yields quasi-optimal approximations; that is,
M
u ’ Rh u ¤ inf u ’ vh .
V V
± vh ∈Vh
Here M and ± denote the constants in the continuity and ellipticity
conditions, respectively.

7.7 Let u ∈ C 1 ([0, T ], V ). Show that Rh u ∈ C 1 ([0, T ], V ) and d
dt Rh u(t) =
d
Rh dt u(t).

7.8 Transfer the derivation of the ¬nite volume method given in Sec-
tion 6.2.2 for the case of an elliptic boundary value problem to the parabolic
initial-boundary value problem (7.1) in divergence form; i.e., convince your-
self that the formalism of obtaining (7.48) indeed can be interpreted as a
¬nite volume semidiscretization of (7.1).



7.3 Fully Discrete Schemes
As we have seen, the application of the vertical method of lines results in
the following situation:
• There is a linear system of ordinary di¬erential equations of high
order (dimension) to be solved.
• There is an error estimate for the solution u of the initial-boundary
value problem (7.1) by means of the exact solution uh of the system
(7.41).
A di¬culty in the choice and in the analysis of an appropriate discretiza-
tion method for systems of ordinary di¬erential equations is that many
standard estimates involve the Lipschitz constant of the corresponding
right-hand side, here q h ’ Ah uh (cf. (7.41)). But this constant is typically
large for small spatial parameters h, and so we would obtain nonrealistic
error estimates (cf. Theorem 3.45).
There are two alternatives. For comparatively simple time discretiza-
tions, certain estimates can be derived in a direct way (i.e., without using
standard estimates for systems of ordinary di¬erential equations). The sec-
ond way is to apply speci¬c time discretizations in conjunction with re¬ned
methods of proof.
312 7. Discretization of Parabolic Problems

Here we will explain the ¬rst way for the so-called one-step-theta method.


One-Step Discretizations in Time, in Particular for the Finite
Di¬erence Method
We start from the problem (7.41), which resulted from spatial discretization
techniques. Provided that T < ∞, the time interval (0, T ) is subdivided
into N ∈ N subintervals of equal length „ := T /N. Furthermore, we set
tn := n„ for n ∈ {0, . . . , N } and un ∈ RM for an approximation of uh (tn ).
h
If the time interval is unbounded, the time step „ > 0 is given, and the
number n ∈ N is allowed to increase without bounded; that is, we set
formally N = ∞.
The values t = tn , where an approximation is to be determined, are
called time levels. The restriction to equidistant time steps is only for the
d
sake of simplicity. We approximate dt uh by the di¬erence quotient
d 1
uh (t) ∼ (uh (tn+1 ) ’ uh (tn )).
dt „
If we interpret this approximation to be at t = tn , we take the forward
di¬erence quotient; at t = tn+1 we take the backward di¬erence quotient;
at t = tn + 1 „ we take the symmetric di¬erence quotient. Again we obtain
2
a generalization and uni¬cation by introducing a parameter ˜ ∈ [0, 1] and
interpreting the approximation to be taken at t = tn + ˜„. As for ˜ = 0
or 1, we are not at a time level, and so we need the further approximation
Ah uh ((n + ˜) „ ) ∼ ˜Ah uh (tn ) + ˜Ah uh (tn+1 ).
Here we use the abbreviation ˜ := 1 ’ ˜. The (one-step-)theta method
de¬nes a sequence of vectors u0 , . . . , uN by, for n = 0, 1, . . . , N ’ 1 ,
h h
1 n+1
uh ’ un + ˜Ah un + ˜Ah un+1 = q h ((n + ˜)„ ) , (7.66)
h h h

u0 = u0 .
h

If we apply this discretization to the more general form (7.45), we get
correspondingly
1
Bh un+1 ’ Bh un + ˜Ah un + ˜Ah un+1 = q h ((n + ˜)„ ) .
ˆ ˆ (7.67)
h h
h h

Analogously to (7.45), the more general form can be transformed to (7.66),
’1
assuming that Bh is regular: either by multiplying (7.67) by Bh or in the
T
case of a decomposition Bh = Eh Eh (for a symmetric positive de¬nite Bh )
’T
by multiplying by Eh and a change of variables to Eh un . We will apply
h
two techniques in the following:
One is based on the eigenvector decomposition of Ah ; thus for (7.67),
this means to consider the generalized eigenvalue problem
ˆ
Ah v = »Bh v . (7.68)
7.3. Fully Discrete Schemes 313

Note that the Galerkin approach for the eigenvalue problems according
to De¬nition 7.6 leads to such a generalized eigenvalue problem with the
ˆ
sti¬ness matrix Ah and the mass matrix Bh .
The other approach is based on the matrix properties (1.32)* or (1.32).
For the most important case,

Bh = diag(bi ) , bi > 0 for i = 1, . . . , M , (7.69)

which corresponds to the mass lumping procedure, the above-mentioned
transformation reduces to a diagonal scaling, which does not in¬‚uence any
of their properties.
Having this in mind, in the following we will consider explicitly only the
formulation (7.66).
In the case ˜ = 0, the explicit Euler method, un can be determined
h
explicitly by

un+1 = „ (q h (tn ) ’ Ah un ) + un = (I ’ „ Ah )un + „ q(tn ) .
h h h
h

Thus the e¬ort for one time step consists of a SAXPY operation, a vector
addition, and a matrix-vector operation. For dimension M the ¬rst of these
is of complexity O(M ), and also the last one if the matrix is sparse in the
sense de¬ned at the beginning of Chapter 5. On the other hand, for ˜ = 0,
the method is implicit, since for each time step a system of linear equations
has to be solved with the system matrix I + ˜„ Ah . Here the cases ˜ = 1,
the implicit Euler method, and ˜ = 1 , the Crank“Nicolson method, will
2
be of interest. Due to our restriction to time-independent coe¬cients, the
matrix is the same for every time step (for constant „ ). If direct methods
(see Section 2.5) are used, then the LU factorization has to be computed
only once, and only forward and backward substitutions with changing
right-hand sides are necessary, where computation for ˜ = 1 also requires
a matrix-vector operation. For band matrices, for example, operations of
the complexity bandwidth — dimension are necessary, which means for the
basic example of the heat equation on a rectangle O(M 3/2 ) operations in-
stead of O(M ) for the explicit method. Iterative methods for the resolution
of (7.66) cannot make use of the constant matrix, but with un there is a
h
good initial iterate if „ is not too large.
Although the explicit Euler method ˜ = 0 seems to be attractive, we will
see later that with respect to accuracy or stability one may prefer ˜ = 1 2
or ˜ = 1.
To investigate further the theta method, we resolve recursively the
relations (7.66) to gain the representation
n
’1
I ’ ˜„ Ah
un = (I + ˜„ Ah ) u0 (7.70)
h
n n’k
’1 ’1
I ’ ˜„ Ah q h tk ’ ˜„ .
+„ (I + ˜„ Ah ) (I + ˜„ Ah )
k=1
314 7. Discretization of Parabolic Problems

Here we use the abbreviation A’n = (A’1 )n for a matrix A. Compar-
ing this with the solution (7.55) of the semidiscrete problem, we see the
approximations

e’Ah tn ∼ Eh,„ ,
n


where
’1
I ’ ˜„ Ah
Eh,„ := (I + ˜„ Ah )

and
tn tn
(tn ’s)/„
e’Ah (tn ’s) q h (s)ds e’Ah „
= q h (s)ds
0 0
n
(t ’s)/„
(I + ˜„ Ah )’1 q h s ’ ˜„ .
∼„ n
Eh,„
k=1
s=k„

(7.71)

The matrix Eh,„ thus is the solution operator of (7.66) for one time step
and homogeneous boundary conditions and right-hand side. It is to be
expected that it has to capture the qualitative behaviour of e’Ah „ that it
is approximating. This will be investigated in the next section.

One-Step Discretizations for the Finite Element Method
The fully discrete scheme can be achieved in two ways: Besides apply-
ing (7.66) to (7.41) in the transformed variable or in the form (7.67), the
discretization approach can applied directly to (7.44):
With ‚U n+1 := (U n+1 ’ U n )/„, f n+s := sf (tn+1 ) + (1 ’ s)f (tn ),
bn+s (v) := sb(tn+1 , v) + (1 ’ s)b(tn , v), b according to (7.43), s ∈ [0, 1],
and with a ¬xed number ˜ ∈ [0, 1], the fully discrete method for (7.44)
then reads as follows:
Find a sequence U 0 , . . . , U N ∈ Vh such that for n ∈ {0, . . . , N ’ 1},

‚U n+1 , vh + a(˜U n+1 + ˜U n , vh ) = bn+˜ (vh )
0
for all vh ∈ Vh , (7.72)
U0 = u0h .

An alternative choice for the right-hand side, closer to the ¬nite di¬erence
method, is the direct evaluation at tn + ˜„, e.g., f (tn + ˜„ ). The version
here is chosen to simplify the order of convergence estimate in Section 7.6.
By representing the U n by means of a basis of Vh as after (7.44), again
we get the form (7.67) (or (7.66) in the transformed variable). Note that
also for ˜ = 0 the problem here is implicit if Bh is not diagonal. Therefore,
mass lumping is often applied, and the scalar product ·, · 0 in (7.72) is
7.4. Stability 315

replaced by an approximation due to numerical quadrature, i.e.,
‚U n+1 , vh + a(˜U n+1 + ˜U n , vh ) = bn+˜ (vh )
0,h
for all vh ∈ Vh , (7.73)
U0 = u0h .
As explained in Section 3.5.2, uh , vh 0,h is the sum over all contributions
from elements K ∈ Th , which takes the form (3.112) for the reference
element. In the case of Lagrange elements and a nodal quadrature formula
we have for the nodal basis functions •i :
•j , •i = •i , •i δij =: bi δij for i, j = 1, . . . , M, (7.74)
0,h 0,h

since for i = j the integrand •i •j vanishes at all quadrature points. In this
case we arrive at the form (7.67) with a matrix Bh satisfying (7.69).

One-Step Discretizations for the Finite Volume Method
As in the previous subsection on the ¬nite element approach, the
semidiscrete formulation (7.48) can be discretized in time directly:
Find a sequence U 0 , . . . , U N ∈ Vh such that for n ∈ {0, . . . , N ’ 1},
‚U n+1 , vh + ah (˜U n+1 + ˜U n , vh ) = f n+˜ , vh 0,h
0,h
for all vh ∈ Vh , (7.75)
U0 = u0h ,
where ‚U n+1 , ˜, f n+˜ are de¬ned as before (7.72).
Remember that here we consider only homogeneous boundary conditions.
If the elements U n , U n+1 are represented by means of a basis of Vh , we
recover the form (7.67).
Since the mass matrix Bh is diagonal, the problem can be regarded as
being explicit for ˜ = 0.


Exercise
7.9 Consider linear simplicial elements de¬ned on a general conforming
triangulation of a polygonally bounded domain „¦ ‚ R2 .
(a) Determine the entries of the mass matrix Bh .
(b) Using the trapezoidal rule, determine the entries of the lumped mass
matrix diag(bi ).


7.4 Stability
In Section 7.3 we have seen that at least if a basis of eigenvectors of the
discretization matrix Ah allows for the solution representation (7.55) for
316 7. Discretization of Parabolic Problems

the semidiscrete method, the qualitative behaviour of e’Ah „ u0 should be
preserved by Eh,„ u0 , being one time step „ for homogeneous boundary
conditions and right-hand side (q h = 0) in the semi- and fully discrete
cases. It is su¬cient to consider the eigenvector w i instead of a general u0 .
Thus, we have to compare
e’Ah „ w i = (e’»i „ )w i (7.76)
with
1 ’ ˜„ »i
’1
I ’ ˜„ Ah
(I + ˜„ Ah ) wi = wi . (7.77)
1 + ˜„ »i
We see that the exponential function is approximated by
1 + (1 ’ ˜)z
R(z) = , (7.78)
1 ’ ˜z
the stability function, at the points z = ’»i „ ∈ C, given by the eigenvalues
»i , and the time step size „ .
For n time steps and q h = 0 we have
n
e’Ah „ wi = e’»i tn wi ∼ R(’»i „ )n w i . (7.79)
Thus, the restriction to eigenvectors wi with eigenvalues »i has diagonal-
ized the system of ordinary di¬erential equations (7.41) for q h = 0 to the
scalar problems
t ∈ (0, T ) ,
ξ + »i ξ = 0, (7.80)
ξ(0) = ξ0
(for ξ0 = 1) with its solution ξ(t) = e’»i tξ0 , for which the one-step-theta
method gives the approximation
ξn+1 = R(’»i „ )ξn = (R(’»i „ ))n+1 ξ0 (7.81)
at t = tn+1 . A basic requirement for a discretion method is the following:
De¬nition 7.17 A one-step method is called nonexpansive if for two nu-
merical approximations un and un , generated under the same conditions
˜h
h
˜
except for two discrete initial values u0 and u0 , respectively, the following
estimate is valid:
|un+1 ’ un+1 | ¤ |un ’ un | , n ∈ {0, . . . , N ’ 1} .
˜h ˜
h

A recursive application of this estimate immediately results in
|un ’ un | ¤ |u0 ’ u0 | , n ∈ {1, . . . , N } .
˜ ˜
Here a general one-step method has the form
n ∈ {0, . . . , N ’ 1} ,
un+1 = un + „ ¦(„, tn , un ) ,
h h
h

with u0 = u0 and a so-called generating function ¦ : R+ — [0, T ) — RM ’
h
RM that characterizes the particular method. The generating function of
7.4. Stability 317

the one-step-theta method applied to the system (7.41) is

<< . .

. 44
( : 59)



. . >>