„¦ 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)

ij

„¦

M¯

¯

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 .

0

„¦

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

type

M

ωi f (ai ) for f ∈ C(„¦)

¯

Q(f ) := (8.29)

i=1

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

for f ∈ C(„¦) ,

¯

Q(f ) := I(f ) dx (8.30)

„¦

where

M

I : C(„¦) ’ Vh ,

¯ I(f ) := f (ai )•i ,

i=1

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,

1

|K| ,

ωi = (8.31)

d+1

K∈Th

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)

gives

˜ ˜ ˜

G(ξ) = Gj (ξ) with Gj (ξ) = ωj ψ(ξ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,

(8.33)

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) .

2

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

condition

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

equations

d

ξ(t) + Sξ(t) + G(ξ(t)) = β(t) for t ∈ (0, T ] ,

B ξ(0) = ξ0

dt

M

for the representing vector ξ(t) of the approximation uh (·, t) = i=1 ξi (t)•i ,

M

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+1

„n U

(8.37)

0

, vh = ˜β(tn+1 ) + (1 ’ ˜)β(tn ).

n+˜

+ψ

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)

or

ψ 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)

with

(B ’ (1 ’ ˜)„n S) ξ n ’ (1 ’ ˜)„n G(ξ n )

b :=

(8.42)

+ ˜β(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

matrix

˜

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

’1

A’1 I + ˜„n B ’1 S B ’1 ˜„n sup |ψ (s)| B

¤ ˜ ˜ ˜

sup D˜(ξ)

g

s∈R

ξ∈R M

’1

I + ˜„n B ’1 S

¤ ˜„n sup |ψ (s)| κ(B)

˜ ˜

s∈R

’1

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

1

’1

¤

˜

I + ˜„n BS ,

1’‘

and thus we obtain

C„n

γ=

1’‘

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.

Exercises

˜

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

9

Discretization Methods

for Convection-Dominated Problems

9.1 Standard Methods and Convection-Dominated

Problems

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

e

number

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

convection-dominated.

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

e

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 „¦,

(9.2)

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

e

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