<< . .

. 51
( : 59)

. . >>

G(ξ) := (Gj (ξ))j with Gj (ξ) := ψ ξi •i •j dx .
„¦ i=1

Note that this notation di¬ers from that in Section 2.2 and the subsequent
chapters: There we denoted S by Ah and b ’ G(ξ) by q h . For reasons of
brevity we omit the index h.
For the moment we want to suppose that the mapping G can be evaluated
exactly. The system of equations (8.26) with
and g(ξ) := G(ξ) ’ b
A := S
8.3. Semilinear Elliptic and Parabolic Problems 361

is of the type introduced in (8.18) in the variable ξ. Thus we may apply,
besides the Newton iteration, the ¬xed-point iteration, introduced in (8.19),
and the variants of Newton™s method, namely the modi¬ed and inexact
versions with their already discussed advantages and drawbacks. We have
to examine the question of how the properties of the matrix will change by
¯ ¯
the transition from A to A + DG(ξ), where ξ stands for the current iterate.
We have
DG(ξ) = ψ (¯)•i •j dx ,
u (8.27)

where u = P ξ = i=1 ξi •i ∈ Vh denotes the function belonging to the
¯ ¯
representing vector ξ. This means that DG(ξ) is symmetric and positive
semide¬nite, respectively de¬nite, if the following condition for ± = 0,
respectively ± > 0, holds:
There exists some ± ≥ 0 such that ψ (u) ≥ ± for all u ∈ R . (8.28)
More precisely, we have for · ∈ RM , if (8.28) is valid,
2 2
¯ ψ (¯) |P ·| dx ≥ ± P ·
·T DG(ξ)· = u .
For such a monotone nonlinearity the properties of de¬niteness of the sti¬-
ness matrix S may be “enforced”. If, on the other hand, we want to make
use of the properties of an M-matrix that can be ensured by the conditions
(1.32) or (1.32)— , then it is not clear whether these properties are conserved
¯ ¯
after addition of DG(ξ). This is due to the fact that DG(ξ) is a sparse ma-
trix of the same structure as S, but it also entails a spatial coupling that
is not contained in the continuous formulation (8.25).

Numerical Quadrature
Owing to the above reason, the use of a node-oriented quadrature rule for
the approximation of G(ξ) is suggested, i.e., a quadrature formula of the
ωi f (ai ) for f ∈ C(„¦)
Q(f ) := (8.29)

with weights ωi ∈ R. Such a quadrature formula results from

for f ∈ C(„¦) ,
Q(f ) := I(f ) dx (8.30)
I : C(„¦) ’ Vh ,
¯ I(f ) := f (ai )•i ,

is the interpolation operator of the degrees of freedom. For this considera-
tion we thus assume that only Lagrangian elements enter the de¬nition of
362 8. Iterative Methods for Nonlinear Equations

Vh . In the case of (8.30) the weights in (8.29) are hence given by

ωi = •i dx .

This corresponds to the local description (3.116). More speci¬cally, we get,
for example, for the linear approach on simplices as a generalization of the
composite trapezoidal rule,
|K| ,
ωi = (8.31)
with ai ∈K

with d denoting the spatial dimension and Th the underlying triangulation.
Approximation of the mapping G by a quadrature rule of the type (8.29)
˜ ˜ ˜
G(ξ) = Gj (ξ) with Gj (ξ) = ωj ψ(ξj ) ,

because of •j (ai ) = δij . We see that the approximation G has the property
˜ ˜
that Gj depends only on ξj . We call such a G a diagonal ¬eld. Qualitatively,
this corresponds better to the continuous formulation (8.25) and leads to
the fact that DG(ξ) is diagonal:
˜¯ ¯
DG(ξ)ij = ωj ψ (ξj )δij . (8.32)
If we impose that all quadrature weights ωi are positive, which is the case
in (8.31) and also in other examples in Section 3.5.2, all of the above con-
˜¯ ˜¯
siderations about the properties of DG(ξ) and S + DG(ξ) remain valid;

additionally, if S is an M-matrix, because the conditions (1.32) or (1.32)
are ful¬lled, then S + DG(ξ) remains an M-matrix, too. This is justi¬ed
by the following fact (compare [34] and [5]; cf. (1.33) for the notation):
If A is an M-matrix and B ≥ A with bij ¤ 0 for i = j,
then B is an M-matrix as well.

Conditions of Convergence
Comparing the requirements for the ¬xed-point iteration and Newton™s
method stated in the (convergence) Theorems 8.4 and 8.11, we observe that
the conditions in Theorem 8.4 can be ful¬lled only in special cases, where
S ’1 DG(ξ) is small according to a suitable matrix norm (see Lemma 8.3).
But it is also di¬cult to draw general conclusions about requirement (iii)
in Theorem 8.11, which together with h < 1 quanti¬es the closeness of the
initial iterate to the solution. The postulation (i), on the other hand, is
met for (8.27) and (8.32) if ψ is Lipschitz continuous (see Exercise 8.7).
Concerning the postulation (ii) we have the following: Let ψ be monotone
nondecreasing (i.e., (8.28) holds with ± ≥ 0) and let S be symmetric and
positive de¬nite, which is true for a problem without convection terms
8.3. Semilinear Elliptic and Parabolic Problems 363

(compare (3.27)). Then we have in the spectral norm
S ’1 = 1/»min(S) .

Here »min (S) > 0 denotes the smallest eigenvalue of S. Hence
(S + DG(ξ))’1 = 1/»min S + DG(ξ) ¤ 1/»min(S) = S ’1
¯ ,
2 2

and consequently, requirement (ii) is valid with β = S ’1 2 .
Concerning the choice of the initial iterate, there is no generally successful
strategy. We may choose the solution of the linear subproblem, i.e.,
Sξ(0) = b . (8.34)
Should it fail to converge even with damping, then we may apply, as a
generalization of (8.34), the continuation method to the family of problems
f (», ξ) := S + »G(ξ) ’ b = 0
with continuation parameter » ∈ [0, 1]. If all these problems have solutions
ξ = ξ» so that Df (ξ; ») exists and is nonsingular in a neighbourhood of
ξ» , and if there exists a continuous solution trajectory without bifurcation,
then [0, 1] can be discretized by 0 = »0 < »1 < · · · < »N = 1, and solutions
ξ»i of f (ξ; »i ) = 0 can be obtained by performing a Newton iteration with
the (approximative) solution for » = »i’1 as starting iterate. Since the ξ »i
for i < N are just auxiliary means, they should be obtained rather coarsely,
i.e., with one or two Newton steps. The stated conditions are ful¬lled under
the supposition (8.28). If this condition of monotonicity does not hold, we
may encounter a bifurcation of the continuous solution (see, for example,
[29, pp. 28 ¬.]).

Instationary Problems
The elliptic boundary value problem (8.25) corresponds to the parabolic
initial value problem
‚t u(x, t) + Lu(x, t) + ψ(u(x, t)) = 0 for (x, t) ∈ QT (8.35)
with linear boundary conditions according to (3.18)“(3.20) and the initial
for x ∈ „¦ .
u(x, 0) = u0 (x) (8.36)
We have already met an example for (8.35), (8.36) in (0.32). Analo-
gously to (8.26) and (7.45), the Galerkin discretization in Vh (i.e., the
semidiscretization) leads to the nonlinear system of ordinary di¬erential
ξ(t) + Sξ(t) + G(ξ(t)) = β(t) for t ∈ (0, T ] ,
B ξ(0) = ξ0
for the representing vector ξ(t) of the approximation uh (·, t) = i=1 ξi (t)•i ,
where u0h = i=1 ξ0i •i is an approximation of the initial value u0 (see
364 8. Iterative Methods for Nonlinear Equations

Section 7.2). The matrix B is the mass matrix
B= •j , •i ,
0 ij

and β(t) contains the contributions of the inhomogeneous boundary
conditions analogously to b in (8.26).
To obtain the fully discrete scheme we use the one-step-theta method as
in Section 7.3. Here we allow the time step size „n to vary in each step, in
particular determined by a time step control before the execution of the
nth time step. So, if the approximation U n is known for t = tn , then the
approximation U n+1 for t = tn+1 := tn + „n is given in generalization of
(7.72) as the solution of
1 ’ U n , vh + a ˜U n+1 + (1 ’ ˜)U n , vh
„n U
, vh = ˜β(tn+1 ) + (1 ’ ˜)β(tn ).

Here ˜ ∈ [0, 1] is the ¬xed parameter of implicity. For the choice of ψ n+˜
we have two possibilities:
ψ n+˜ = ˜ψ(U n+1 ) + (1 ’ ˜)ψ(U n ) (8.38)
ψ n+˜ = ψ ˜U n+1 + (1 ’ ˜)U n . (8.39)
In the explicit case, i.e., ˜ = 0, (8.37) represents a linear system of equa-
tions for U n+1 (with the system matrix B) and does not have to be treated
further here. In the implicit case ˜ ∈ (0, 1] we obtain again a nonlinear
system of the type (8.18), i.e.,
Aξ + g(ξ) = 0 ,
in the variable ξ = ξn+1 , where ξn+1 is the representation vector of U n+1 :
M n+1
U n+1 = i=1 ξi •i . Now we have for the variant (8.38),
A := B + ˜„n S , (8.40)
:= ˜„n G(ξ) ’ b ,
g(ξ) (8.41)
(B ’ (1 ’ ˜)„n S) ξ n ’ (1 ’ ˜)„n G(ξ n )
b :=
+ ˜β(tn+1 ) + (1 ’ ˜)β(tn ) .
For the variant (8.39) g changes to
g(ξ) := „n G (˜ξ + (1 ’ ˜)ξn ) ’ b ,
and in the de¬nition of b the second summation term drops out. The vector
ξn is the representation vector of the already known approximation U n .
8.3. Semilinear Elliptic and Parabolic Problems 365

Numerical Quadrature
As in the stationary case we can approximate g by a quadrature rule of the
form (8.29), which leads to
g(ξ) = ˜„n G(ξ) ’ b
in (8.38) and to
g (ξ) = „n G (˜ξ + (1 ’ ˜)ξn ) ’ b
in (8.39). The functional matrices of g and g are thus equal for (8.38) and
(8.39), except to the point where ψ is being evaluated. Consequently, it
su¬ces in the following to refer to (8.38). Based on the same motivation,
a quadrature rule of the form (8.29) can be applied to the mass matrix
B. Such a mass lumping results in a diagonal approximation of the mass
B = diag(ωi ) .
In contrast to the stationary case we get the factor ˜„n in front of the
nonlinearity, where the time step size „n may be chosen arbitrarily small.
Of course, we have to take into account that the number of time steps
necessary to achieve a ¬xed time T is respectively raised. All of the above
considerations about the matrix properties of A + Dg(ξ) are conserved,
where A is no longer the sti¬ness matrix, but represents the linear combina-
tion (8.40) with the mass matrix. This reduces the requirements concerning
the V -ellipticity of a (see (3.27)) and thus the positive de¬niteness of A.
Admittedly, A is not necessarily an M-matrix if S is one, because the
conditions (1.32) or (1.32)— are not valid. Here the approximation B is ad-
vantageous, because using nonnegative weights will conserve this property
due to (8.33).

Conditions of Convergence
Clear di¬erences arise in answering the question of how to ensure the con-
vergence of the iteration schemes. Even for the ¬xed-point iteration it is
true that the method converges globally if only the time step size „n is
chosen small enough. We want to demonstrate this in the following by an
example of a quadrature with nonnegative weights in the mass matrix and
the nonlinearity. Therefore, the Lipschitz constant of A’1 g is estimated
according to Lemma 8.3. Let the norm be a matrix norm induced by a
p-norm | · |p and let A be nonsingular. We get
A’1 I + ˜„n B ’1 S B ’1 ˜„n sup |ψ (s)| B
¤ ˜ ˜ ˜
sup D˜(ξ)
ξ∈R M

I + ˜„n B ’1 S
¤ ˜„n sup |ψ (s)| κ(B)
˜ ˜
I + ˜„n B ’1 S
=: C„n .
366 8. Iterative Methods for Nonlinear Equations

Thus we assume the boundedness of ψ on R (which may even be weakened).
For a given ‘ ∈ (0, 1) choose „n su¬ciently small such that
˜„n B ’1 S ¤ ‘

holds. With Lemma (A3.11) it follows that
I + ˜„n BS ,
and thus we obtain
as a Lipschitz constant for A’1 g. We see that by choosing „n su¬ciently
small, the contraction property of A’1 g can be guaranteed. From this fact
a (heuristic) step size control can be deduced that reduces the step size
when a lack of convergence is detected and repeats the step, and in case of
satisfactory convergence increases the time step size.
Nevertheless, in general, Newton™s method is to be preferred: Here we
can expect that the quality of the initial iterate ξ(0) = ξ n for time step
(n + 1) improves the smaller we choose „n . The step size control mentioned
above may thus be chosen here, too (in conjunction with the enlargement
of the range of convergence via damping). Nonetheless, a problem only to
be solved in numerical practice consists in coordinating the control param-
eters of the time step control, the damping strategy, and eventually the
termination of the inner iteration in such a way that overall, an e¬cient
algorithm is obtained.

8.7 Study the Lipschitz property of DG de¬ned by (8.27) and of DG
de¬ned by (8.32), provided ψ is Lipschitz.

8.8 Decide whether A’1 g is contractive in case of (8.40)“(8.42).

8.9 The boundary value problem
’u + eu = 0 in (0, 1), u(0) = u(1) = 0,
is to be discretized by a ¬nite element method using continuous, piecewise
linear functions on equidistant grids. Quadrature is to be done with the
trapezoidal rule.
(a) Compute the matrix Ah ∈ Rm,m and the nonlinear vectorvalued
function Fh : Rm ’ Rm , in a matrix-vector notation
Ah Uh + Fh (Uh ) = 0
8.3. Semilinear Elliptic and Parabolic Problems 367

of the discretization. Here Uh ∈ Rm denotes the vector of un-
known nodal values of the approximative solution and, for uniqueness
of the representation, the elements of Ah are independent of the
discretization parameter h.
(b) Study the convergence of the iterative procedure
(k+1) (k) (k)
= (2 + h2 )I ’ Ah Uh ’ Fh Uh
(2 + h2 )Uh
(±) ,
(k+1) (k+1) (k)
= (2I ’ Ah ) Uh .
(β) 2Uh + Fh Uh
Discretization Methods
for Convection-Dominated Problems

9.1 Standard Methods and Convection-Dominated
As we have seen in the introductory Chapter 0, the modelling of transport
and reaction processes in porous media results in di¬erential equations of
the form
‚t u ’ ∇ · (K∇u ’ cu) = f ,
which is a special case of the form (0.33). Similar equations occur in the
modelling of the heat transport in ¬‚owing water, the carrier transport
in semiconductors, and the propagation of epidemics. These application-
speci¬c equations often share the property that their so-called global P´clet
c ∞ diam(„¦)
Pe := (9.1)
K ∞

is signi¬cantly larger than one. For example, representative values range
from 25 (transport of a dissolved substance in ground water) up to about
107 (modelling of semiconductors). In such cases, the equations are called
Therefore, in what follows, the Dirichlet boundary value problem intro-
duced in Section 3.2 will be looked at from the point of view of large global
P´clet numbers, whereas in Section 9.4, the initial boundary value problem
from Chapter 7 will be considered from this aspect.
9.1. Standard Methods 369

Let „¦ ‚ Rd denote a bounded domain with a Lipschitz continuous
boundary. Given a function f : „¦ ’ R, a function u : „¦ ’ R is to
be determined such that
Lu =f in „¦,
u =0 on “,
where again
Lu := ’∇ · (K∇u) + c · ∇u + ru ,
with su¬ciently smooth coe¬cients
K : „¦ ’ Rd,d , c : „¦ ’ Rd , r : „¦ ’ R.
Unfortunately, standard discretization methods (¬nite di¬erence, ¬nite
element, and ¬nite volume methods) fail when applied to convection-
dominated equations. At ¬rst glance, this seems to be a contradiction to the
theory of these methods presented in the preceding chapters, because there
we did not have any restriction on the global P´clet number. This apparent
contradiction may be explained as follows: On the one hand, the theoretical
results are still true for the convection-dominated case, but on the other
hand, some assumptions of the statements therein (such as “for su¬ciently
small h”) lack sharpness. This, in turn, may lead to practically unrealistic
conditions (cf. the later discussion of the estimate (9.13)). For example, it

<< . .

. 51
( : 59)

. . >>