. 1
( 16)


Numerical Methods for
Elliptic and Parabolic
Partial Differential

Peter Knabner
Lutz Angermann

Texts in Applied Mathematics

J.E. Marsden
L. Sirovich
S.S. Antman

G. Iooss
P. Holmes
D. Barkley
M. Dellnitz
P. Newton
This page intentionally left blank
Peter Knabner Lutz Angermann

Numerical Methods for
Elliptic and Parabolic Partial
Differential Equations

With 67 Figures
Peter Knabner Lutz Angermann
Institute for Applied Mathematics Institute for Mathematics
University of Erlangen University of Clausthal
Martensstrasse 3 Erzstrasse 1
D-91058 Erlangen D-38678 Clausthal-Zellerfeld
Germany Germany
knabner@am.uni-erlangen.de angermann@math.tu-clausthal.de
Series Editors
J.E. Marsden L. Sirovich
Control and Dynamical Systems, 107“81 Division of Applied Mathematics
California Institute of Technology Brown University
Pasadena, CA 91125 Providence, RI 02912
marsden@cds.caltech.edu chico@camelot.mssm.edu
S.S. Antman
Department of Mathematics
Institute for Physical Science
and Technology
University of Maryland
College Park, MD 20742-4015
Mathematics Subject Classification (2000): 65Nxx, 65Mxx, 65F10, 65H10

Library of Congress Cataloging-in-Publication Data
Knabner, Peter.
[Numerik partieller Differentialgleichungen. English]
Numerical methods for elliptic and parabolic partial differential equations /
Peter Knabner, Lutz Angermann.
p. cm. ” (Texts in applied mathematics ; 44)
Include bibliographical references and index.
ISBN 0-387-95449-X (alk. paper)
1. Differential equations, Partial”Numerical solutions. I. Angermann, Lutz. II. Title.
III. Series.
QA377.K575 2003
515′.353”dc21 2002044522

ISBN 0-387-95449-X Printed on acid-free paper.

™ 2003 Springer-Verlag New York, Inc.
All rights reserved. This work may not be translated or copied in whole or in part without the written
permission of the publisher (Springer-Verlag New York, Inc., 175 Fifth Avenue, New York, NY 10010,
USA), except for brief excerpts in connection with reviews or scholarly analysis. Use in connection with
any form of information storage and retrieval, electronic adaptation, computer software, or by similar or
dissimilar methodology now known or hereafter developed is forbidden.
The use in this publication of trade names, trademarks, service marks, and similar terms, even if they are
not identified as such, is not to be taken as an expression of opinion as to whether or not they are subject
to proprietary rights.

Printed in the United States of America.

987654321 SPIN 10867187

Typesetting: Pages created by the authors in 2e using Springer™s svsing6.cls macro.


Springer-Verlag New York Berlin Heidelberg
A member of BertelsmannSpringer Science+Business Media GmbH
Series Preface

Mathematics is playing an ever more important role in the physical and
biological sciences, provoking a blurring of boundaries between scientific
disciplines and a resurgence of interest in the modern as well as the classical
techniques of applied mathematics. This renewal of interest, both in re-
search and teaching, has led to the establishment of the series Texts in
Applied Mathematics (TAM).
The development of new courses is a natural consequence of a high level
of excitement on the research frontier as newer techniques, such as numeri-
cal and symbolic computer systems, dynamical systems, and chaos, mix
with and reinforce the traditional methods of applied mathematics. Thus,
the purpose of this textbook series is to meet the current and future needs
of these advances and to encourage the teaching of new courses.
TAM will publish textbooks suitable for use in advanced undergraduate
and beginning graduate courses, and will complement the Applied Mathe-
matical Sciences (AMS) series, which will focus on advanced textbooks and
research-level monographs.

Pasadena, California J.E. Marsden
Providence, Rhode Island L. Sirovich
College Park, Maryland S.S. Antman
This page intentionally left blank
Preface to the English Edition

Shortly after the appearance of the German edition we were asked by
Springer to create an English version of our book, and we gratefully ac-
cepted. We took this opportunity not only to correct some misprints and
mistakes that have come to our knowledge1 but also to extend the text at
various places. This mainly concerns the role of the ¬nite di¬erence and
the ¬nite volume methods, which have gained more attention by a slight
extension of Chapters 1 and 6 and by a considerable extension of Chapter
7. Time-dependent problems are now treated with all three approaches (¬-
nite di¬erences, ¬nite elements, and ¬nite volumes), doing this in a uniform
way as far as possible. This also made a reordering of Chapters 6“8 nec-
essary. Also, the index has been enlarged. To improve the direct usability
in courses, exercises now follow each section and should provide enough
material for homework.
This new version of the book would not have come into existence without
our already mentioned team of helpers, who also carried out ¬rst versions
of translations of parts of the book. Beyond those already mentioned, the
team was enforced by Cecilia David, Basca Jadamba, Dr. Serge Kr¨utle, a
Dr. Wilhelm Merz, and Peter Mirsch. Alexander Prechtel now took charge
of the di¬cult modi¬cation process. Prof. Paul DuChateau suggested im-
provements. We want to extend our gratitude to all of them. Finally, we

1 Users
of the German edition may consult
viii Preface to the English Edition

thank senior editor Achi Dosanjh, from Springer-Verlag New York, Inc., for
her constant encouragement.

Remarks for the Reader and the Use in Lectures
The size of the text corresponds roughly to four hours of lectures per week
over two terms. If the course lasts only one term, then a selection is nec-
essary, which should be orientated to the audience. We recommend the
following “cuts”:
Chapter 0 may be skipped if the partial di¬erential equations treated
therein are familiar. Section 0.5 should be consulted because of the notation
collected there. The same is true for Chapter 1; possibly Section 1.4 may
be integrated into Chapter 3 if one wants to deal with Section 3.9 or with
Section 7.5.
Chapters 2 and 3 are the core of the book. The inductive presenta-
tion that we preferred for some theoretical aspects may be shortened for
students of mathematics. To the lecturer™s taste and depending on the
knowledge of the audience in numerical mathematics Section 2.5 may be
skipped. This might impede the treatment of the ILU preconditioning in
Section 5.3. Observe that in Sections 2.1“2.3 the treatment of the model
problem is merged with basic abstract statements. Skipping the treatment
of the model problem, in turn, requires an integration of these statements
into Chapter 3. In doing so Section 2.4 may be easily combined with Sec-
tion 3.5. In Chapter 3 the theoretical kernel consists of Sections 3.1, 3.2.1,
Chapter 4 presents an overview of its subject, not a detailed development,
and is an extension of the classical subjects, as are Chapters 6 and 9 and
the related parts of Chapter 7.
In the extensive Chapter 5 one might focus on special subjects or just con-
sider Sections 5.2, 5.3 (and 5.4) in order to present at least one practically
relevant and modern iterative method.
Section 8.1 and the ¬rst part of Section 8.2 contain basic knowledge of
numerical mathematics and, depending on the audience, may be omitted.
The appendices are meant only for consultation and may complete
the basic lectures, such as in analysis, linear algebra, and advanced
mathematics for engineers.
Concerning related textbooks for supplementary use, to the best of our
knowledge there is none covering approximately the same topics. Quite a
few deal with ¬nite element methods, and the closest one in spirit probably
is [21], but also [6] or [7] have a certain overlap, and also o¬er additional
material not covered here. From the books specialised in ¬nite di¬erence
methods, we mention [32] as an example. The (node-oriented) ¬nite volume
method is popular in engineering, in particular in ¬‚uid dynamics, but to
the best of our knowledge there is no presentation similar to ours in a
Preface to the English Edition ix

mathematical textbook. References to textbooks specialised in the topics
of Chapters 4, 5 and 8 are given there.

Remarks on the Notation
Printing in italics emphasizes de¬nitions of notation, even if this is not
carried out as a numbered de¬nition.
Vectors appear in di¬erent forms: Besides the “short” space vectors
x ∈ Rd there are “long” representation vectors u ∈ Rm , which describe
in general the degrees of freedom of a ¬nite element (or volume) approxi-
mation or represent the values on grid points of a ¬nite di¬erence method.
Here we choose bold type, also in order to have a distinctive feature from
the generated functions, which frequently have the same notation, or from
the grid functions.
Deviations can be found in Chapter 0, where vectorial quantities belong-
ing to Rd are boldly typed, and in Chapters 5 and 8, where the unknowns
of linear and nonlinear systems of equations, which are treated in a general
manner there, are denoted by x ∈ Rm .
Components of vectors will be designated by a subindex, creating a
double index for indexed quantities. Sequences of vectors will be supplied
with a superindex (in parentheses); only in an abstract setting do we use

Erlangen, Germany Peter Knabner
Clausthal-Zellerfeld, Germany Lutz Angermann
January 2002
This page intentionally left blank
Preface to the German Edition

This book resulted from lectures given at the University of Erlangen“
Nuremberg and at the University of Magdeburg. On these occasions we
often had to deal with the problem of a heterogeneous audience composed
of students of mathematics and of di¬erent natural or engineering sciences.
Thus the expectations of the students concerning the mathematical accu-
racy and the applicability of the results were widely spread. On the other
hand, neither relevant models of partial di¬erential equations nor some
knowledge of the (modern) theory of partial di¬erential equations could be
assumed among the whole audience. Consequently, in order to overcome the
given situation, we have chosen a selection of models and methods relevant
for applications (which might be extended) and attempted to illuminate the
whole spectrum, extending from the theory to the implementation, with-
out assuming advanced mathematical background. Most of the theoretical
obstacles, di¬cult for nonmathematicians, will be treated in an “induc-
tive” manner. In general, we use an explanatory style without (hopefully)
compromising the mathematical accuracy.
We hope to supply especially students of mathematics with the in-
formation necessary for the comprehension and implementation of ¬nite
element/¬nite volume methods. For students of the various natural or
engineering sciences the text o¬ers, beyond the possibly already existing
knowledge concerning the application of the methods in special ¬elds, an
introduction into the mathematical foundations, which should facilitate the
transformation of speci¬c knowledge to other ¬elds of applications.
We want to express our gratitude for the valuable help that we received
during the writing of this book: Dr. Markus Bause, Sandro Bitterlich,
xii Preface to the German Edition

Dr. Christof Eck, Alexander Prechtel, Joachim Rang, and Dr. Eckhard
Schneid did the proofreading and suggested important improvements. From
the anonymous referees we received useful comments. Very special thanks
go to Mrs. Magdalena Ihle and Dr. Gerhard Summ. Mrs. Ihle transposed
the text quickly and precisely into TEX. Dr. Summ not only worked on the
original script and on the TEX-form, he also organized the complex and
distributed rewriting and extension procedure. The elimination of many
inconsistencies is due to him. Additionally he in¬‚uenced parts of Sec-
tions 3.4 and 3.8 by his outstanding diploma thesis. We also want to thank
Dr. Chistoph Tapp for the preparation of the graphic of the title and for
providing other graphics from his doctoral thesis [70].
Of course, hints concerning (typing) mistakes and general improvements
are always welcome.
We thank Springer-Verlag for their constructive collaboration.
Last, but not least, we want to express our gratitude to our families for
their understanding and forbearance, which were necessary for us especially
during the last months of writing.

Erlangen, Germany Peter Knabner
Magdeburg, Germany Lutz Angermann
February 2000

Series Preface v

Preface to the English Edition vii

Preface to the German Edition xi

0 For Example: Modelling Processes in Porous Media
with Di¬erential Equations 1
0.1 The Basic Partial Di¬erential Equation Models . . . . . 1
0.2 Reactions and Transport in Porous Media . . . . . . . . 5
0.3 Fluid Flow in Porous Media . . . . . . . . . . . . . . . . 7
0.4 Reactive Solute Transport in Porous Media . . . . . . . . 11
0.5 Boundary and Initial Value Problems . . . . . . . . . . . 14

1 For the Beginning: The Finite Di¬erence Method
for the Poisson Equation 19
1.1 The Dirichlet Problem for the Poisson Equation . . ... 19
1.2 The Finite Di¬erence Method . . . . . . . . . . . . ... 21
1.3 Generalizations and Limitations
of the Finite Di¬erence Method . . . . . . . . . . . ... 29
1.4 Maximum Principles and Stability . . . . . . . . . . ... 36

2 The Finite Element Method for the Poisson Equation 46
2.1 Variational Formulation for the Model Problem . . . . . 46
xiv Contents

2.2 The Finite Element Method with Linear Elements . ... 55
2.3 Stability and Convergence of the
Finite Element Method . . . . . . . . . . . . . . . . ... 68
2.4 The Implementation of the Finite Element Method:
Part 1 . . . . . . . . . . . . . . . . . . . . . . . . . ... 74
2.5 Solving Sparse Systems of Linear Equations
by Direct Methods . . . . . . . . . . . . . . . . . . ... 82

3 The Finite Element Method for Linear Elliptic
Boundary Value Problems of Second Order 92
3.1 Variational Equations and Sobolev Spaces . . . . . . .. 92
3.2 Elliptic Boundary Value Problems of Second Order . .. 100
3.3 Element Types and A¬ne
Equivalent Triangulations . . . . . . . . . . . . . . . .. 114
3.4 Convergence Rate Estimates . . . . . . . . . . . . . . .. 131
3.5 The Implementation of the Finite Element Method:
Part 2 . . . . . . . . . . . . . . . . . . . . . . . . . . .. 148
3.6 Convergence Rate Results in Case of
Quadrature and Interpolation . . . . . . . . . . . . . . . 155
3.7 The Condition Number of Finite Element Matrices . . . 163
3.8 General Domains and Isoparametric Elements . . . . . . 167
3.9 The Maximum Principle for Finite Element Methods . . 171

4 Grid Generation and A Posteriori Error Estimation 176
4.1 Grid Generation . . . . . . . . . . . . . . . . . . . . . . . 176
4.2 A Posteriori Error Estimates and Grid Adaptation . . . 185

5 Iterative Methods for Systems of Linear Equations 198
5.1 Linear Stationary Iterative Methods . . . . . . . . . ... 200
5.2 Gradient and Conjugate Gradient Methods . . . . . ... 217
5.3 Preconditioned Conjugate Gradient Method . . . . ... 227
5.4 Krylov Subspace Methods
for Nonsymmetric Systems of Equations . . . . . . ... 233
5.5 The Multigrid Method . . . . . . . . . . . . . . . . ... 238
5.6 Nested Iterations . . . . . . . . . . . . . . . . . . . ... 251

6 The Finite Volume Method 255
6.1 The Basic Idea of the Finite Volume Method . . . . . . . 256
6.2 The Finite Volume Method for Linear Elliptic Di¬eren-
tial Equations of Second Order on Triangular Grids . . . 262

7 Discretization Methods for Parabolic Initial Boundary
Value Problems 283
7.1 Problem Setting and Solution Concept . . . . . . . . . . 283
7.2 Semidiscretization by the Vertical Method of Lines . . . 293
Contents xv

7.3 Fully Discrete Schemes . . . . . .............. 311
7.4 Stability . . . . . . . . . . . . . .............. 315
7.5 The Maximum Principle for the
One-Step-Theta Method . . . . .............. 323
7.6 Order of Convergence Estimates .............. 330

8 Iterative Methods for Nonlinear Equations 342
8.1 Fixed-Point Iterations . . . . . . . . . . . . . . . . . . . . 344
8.2 Newton™s Method and Its Variants . . . . . . . . . . . . 348
8.3 Semilinear Boundary Value Problems for Elliptic
and Parabolic Equations . . . . . . . . . . . . . . . . . . 360

9 Discretization Methods
for Convection-Dominated Problems 368
9.1 Standard Methods and
Convection-Dominated Problems . . . . . . . . . . . . . 368
9.2 The Streamline-Di¬usion Method . . . . . . . . . . . . . 375
9.3 Finite Volume Methods . . . . . . . . . . . . . . . . . . . 383
9.4 The Lagrange“Galerkin Method . . . . . . . . . . . . . . 387

A Appendices 390
A.1 Notation . . . . . . . . . . . . . . . . . . . ........ 390
A.2 Basic Concepts of Analysis . . . . . . . . . ........ 393
A.3 Basic Concepts of Linear Algebra . . . . . ........ 394
A.4 Some De¬nitions and Arguments of Linear
Functional Analysis . . . . . . . . . . . . . ........ 399
A.5 Function Spaces . . . . . . . . . . . . . . . ........ 404

References: Textbooks and Monographs 409

References: Journal Papers 412

Index 415
This page intentionally left blank
For Example:
Modelling Processes in Porous
Media with Di¬erential Equations

This chapter illustrates the scienti¬c context in which di¬erential equation
models may occur, in general, and also in a speci¬c example. Section 0.1
reviews the fundamental equations, for some of them discretization tech-
niques will be developed and investigated in this book. In Sections 0.2 “
0.4 we focus on reaction and transport processes in porous media. These
sections are independent of the remaining parts and may be skipped by
the reader. Section 0.5, however, should be consulted because it ¬xes some
notation to be used later on.

0.1 The Basic Partial Di¬erential Equation Models
Partial di¬erential equations are equations involving some partial deriva-
tives of an unknown function u in several independent variables. Partial
di¬erential equations which arise from the modelling of spatial (and tempo-
ral) processes in nature or technology are of particular interest. Therefore,
we assume that the variables of u are x = (x1 , . . . , xd )T ∈ Rd for d ≥ 1,
representing a spatial point, and possibly t ∈ R, representing time. Thus
the minimal set of variables is (x1 , x2 ) or (x1 , t), otherwise we have ordinary
di¬erential equations. We will assume that x ∈ „¦, where „¦ is a bounded
domain, e.g., a metal workpiece, or a groundwater aquifer, and t ∈ (0, T ] for
some (time horizon) T > 0. Nevertheless also processes acting in the whole
Rd — R, or in unbounded subsets of it, are of interest. One may consult the
Appendix for notations from analysis etc. used here. Often the function u
2 0. Modelling Processes in Porous Media with Di¬erential Equations

represents, or is related to, the volume density of an extensive quantity like
mass, energy, or momentum, which is conserved. In their original form all
quantities have dimensions that we denote in accordance with the Inter-
national System of Units (SI) and write in square brackets [ ]. Let a be
a symbol for the unit of the extensive quantity, then its volume density
is assumed to have the form S = S(u), i.e., the unit of S(u) is a/m3 . For
example, for mass conservation a = kg, and S(u) is a concentration. For
describing the conservation we consider an arbitrary “not too bad” sub-
set „¦ ‚ „¦, the control volume. The time variation of the total extensive
quantity in „¦ is then

‚t S(u(x, t))dx . (0.1)

If this function does not vanish, only two reasons are possible due to con-
” There is an internally distributed source density Q = Q(x, t, u) [a/m3 /s],
being positive if S(u) is produced, and negative if it is destroyed, i.e., one
term to balance (0.1) is „¦ Q(x, t, u(x, t))dx.
” There is a net ¬‚ux of the extensive quantity over the boundary ‚ „¦ of
„¦. Let J = J (x, t) [a/m2 /s] denote the ¬‚ux density, i.e., J i is the amount,
that passes a unit square perpendicular to the ith axis in one second in
the direction of the ith axis (if positive), and in the opposite direction
otherwise. Then another term to balance (0.1) is given by

’ J (x, t) · ν(x)dσ ,

where ν denotes the outer unit normal on ‚„¦. Summarizing the conserva-
tion reads

S(u(x, t))dx = ’ J (x, t) · ν(x)dσ +
‚t Q(x, t, u(x, t))dx . (0.2)
˜ ˜ ˜
„¦ ‚„¦ „¦

The integral theorem of Gauss (see (2.3)) and an exchange of time
derivative and integral leads to

[‚t S(u(x, t)) + ∇ · J (x, t) ’ Q(x, t, u(x, t))]dx = 0 ,

and, as „¦ is arbitrary, also to
‚t S(u(x, t)) + ∇ · J(x, t) = Q(x, t, u(x, t)) for x ∈ „¦, t ∈ (0, T ] . (0.3)
All manipulations here are formal assuming that the functions involved
have the necessary properties. The partial di¬erential equation (0.3) is the
basic pointwise conservation equation, (0.2) its corresponding integral form.
Equation (0.3) is one requirement for the two unknowns u and J , thus it
0.1. The Basic Partial Di¬erential Equation Models 3

has to be closed by a (phenomenological) constitutive law, postulating a
relation between J and u.
Assume „¦ is a container ¬lled with a ¬‚uid in which a substance is dis-
solved. If u is the concentration of this substance, then S(u) = u and a
= kg. The description of J depends on the processes involved. If the ¬‚uid
is at rest, then ¬‚ux is only possible due to molecular di¬usion, i.e., a ¬‚ux
from high to low concentrations due to random motion of the dissolved
particles. Experimental evidence leads to
J (1) = ’K∇u (0.4)
with a parameter K > 0 [m2 /s], the molecular di¬usivity. Equation (0.4)
is called Fick™s law.
In other situations, like heat conduction in a solid, a similar model occurs.
Here, u represents the temperature, and the underlying principle is energy
conservation. The constitutive law is Fourier™s law, which also has the form
(0.4), but as K is a material parameter, it may vary with space or, for
anisotropic materials, be a matrix instead of a scalar.
Thus we obtain the di¬usion equation
‚t u ’ ∇ · (K∇u) = Q . (0.5)
If K is scalar and constant ” let K = 1 by scaling ”, and f := Q is
independent of u, the equation simpli¬es further to
‚t u ’ ∆u = f ,
where ∆u := ∇·(∇u) . We mentioned already that this equation also occurs
in the modelling of heat conduction, therefore this equation or (0.5) is also
called the heat equation.
If the ¬‚uid is in motion with a (given) velocity c then (forced) convection
of the particles takes place, being described by
J (2) = uc , (0.6)
i.e., taking both processes into account, the model takes the form of the
convection-di¬usion equation
‚t u ’ ∇ · (K∇u ’ cu) = Q . (0.7)
The relative strength of the two processes is measured by the P´clete
number (de¬ned in Section 0.4). If convection is dominating one may ignore
di¬usion and only consider the transport equation
‚t u + ∇ · (cu) = Q . (0.8)
The di¬erent nature of the two processes has to be re¬‚ected in the models,
therefore, adapted discretization techniques will be necessary. In this book
we will consider models like (0.7), usually with a signi¬cant contribution
of di¬usion, and the case of dominating convection is studied in Chapter
9. The pure convective case like (0.8) will not be treated.
4 0. Modelling Processes in Porous Media with Di¬erential Equations

In more general versions of (0.7) ‚t u is replaced by ‚t S(u), where S
depends linearly or nonlinearly on u. In the case of heat conduction S is
the internal energy density, which is related to the temperature u via the
factors mass density and speci¬c heat. For some materials the speci¬c heat
depends on the temperature, then S is a nonlinear function of u.
Further aspects come into play by the source term Q if it depends linearly
or nonlinearly on u, in particular due to (chemical) reactions. Examples for
these cases will be developed in the following sections. Since equation (0.3)
and its examples describe conservation in general, it still has to be adapted
to a concrete situation to ensure a unique solution u. This is done by the
speci¬cation of an initial condition
for x ∈ „¦ ,
S(u(x, 0)) = S0 (x)
and by boundary conditions. In the example of the water ¬lled container
no mass ¬‚ux will occur across its walls, therefore, the following boundary
J · ν(x, t) = 0 for x ∈ ‚„¦, t ∈ (0, T ) (0.9)
is appropriate, which ” depending on the de¬nition of J ” prescribes the
normal derivative of u, or a linear combination of it and u. In Section 0.5
additional situations are depicted.
If a process is stationary, i.e. time-independent, then equation (0.3)
reduces to
∇ · J (x) = Q(x, u(x)) for x ∈ „¦ ,
which in the case of di¬usion and convection is speci¬ed to
’∇ · (K∇u ’ cu) = Q .
For constant K ” let K = 1 by scaling ”, c = 0, and f := Q, being
independent of u, this equation reduces to
’∆u = f in „¦ ,
the Poisson equation.
Instead of the boundary condition (0.9), one can prescribe the values of
the function u at the boundary:
for x ∈ ‚„¦ .
u(x) = g(x)
For models , where u is a concentration or temperature, the physical reali-
sation of such a boundary condition may raise questions, but in mechanical
models, where u is to interpreted as a displacement, such a boundary con-
dition seems reasonable. The last boundary value problem will be the ¬rst
model, whose discretization will be discussed in Chapters 1 and 2.
Finally it should be noted that it is advisable to non-dimensionalise the
¬nal model before numerical methods are applied. This means that both
the independent variables xi (and t), and the dependent one u, are replaced
0.2. Reactions and Transport in Porous Media 5

by xi /xi,ref , t/tref , and u/uref , where xi,ref , tref , and uref are ¬xed reference
values of the same dimension as xi , t, and u, respectively. These reference
values are considered to be of typical size for the problems under investiga-
tion. This procedure has two advantages: On the one hand, the typical size
is now 1, such that there is an absolute scale for (an error in) a quantity
to be small or large. On the other hand, if the reference values are chosen
appropriately a reduction in the number of equation parameters like K
and c in (0.7) might be possible, having only fewer algebraic expressions of
the original material parameters in the equation. This facilitates numerical
parameter studies.

0.2 Reactions and Transport in Porous Media
A porous medium is a heterogeneous material consisting of a solid matrix
and a pore space contained therein. We consider the pore space (of the
porous medium) as connected; otherwise, the transport of ¬‚uids in the
pore space would not be possible. Porous media occur in nature and man-
ufactured materials. Soils and aquifers are examples in geosciences; porous
catalysts, chromatographic columns, and ceramic foams play important
roles in chemical engineering. Even the human skin can be considered a
porous medium. In the following we focus on applications in the geosciences.
Thus we use a terminology referring to the natural soil as a porous medium.
On the micro or pore scale of a single grain or pore, i.e., in a range of µm
to mm, the ¬‚uids constitute di¬erent phases in the thermodynamic sense.
Thus we name this system in the case of k ¬‚uids including the solid matrix
as (k + 1)-phase system or we speak of k-phase ¬‚ow.
We distinguish three classes of ¬‚uids with di¬erent a¬nities to the solid
matrix. These are an aqueous phase, marked with the index “w” for water,
a nonaqueous phase liquid (like oil or gasoline as natural resources or con-
taminants), marked with the index “o,” and a gaseous phase, marked with
the index “g” (e.g., soil air). Locally, at least one of these phases has al-
ways to be present; during a transient process phases can locally disappear
or be generated. These ¬‚uid phases are in turn mixtures of several com-
ponents. In applications of the earth sciences, for example, we do not deal
with pure water but encounter di¬erent species in true or colloidal solu-
tion in the solvent water. The wide range of chemical components includes
plant nutrients, mineral nutrients from salt domes, organic decomposition
products, and various organic and inorganic chemicals. These substances
are normally not inert, but are subject to reactions and transformation
processes. Along with di¬usion, forced convection induced by the motion
of the ¬‚uid is the essential driving mechanism for the transport of solutes.
But we also encounter natural convection by the coupling of the dynamics
of the substance to the ¬‚uid ¬‚ow. The description level at the microscale
6 0. Modelling Processes in Porous Media with Di¬erential Equations

that we have used so far is not suitable for processes at the laboratory or
technical scale, which take place in ranges of cm to m, or even for processes
in a catchment area with units of km. For those macroscales new models
have to be developed, which emerge from averaging procedures of the mod-
els on the microscale. There may also exist principal di¬erences among the
various macroscales that let us expect di¬erent models, which arise from
each other by upscaling. But this aspect will not be investigated here fur-
ther. For the transition of micro to macro scales the engineering sciences
provide the heuristic method of volume averaging, and mathematics the
rigorous (but of only limited use) approach of homogenization (see [36] or
[19]). None of the two possibilities can be depicted here completely. Where
necessary we will refer to volume averaging for (heuristic) motivation.
Let „¦ ‚ Rd be the domain of interest. All subsequent considerations are
formal in the sense that the admissibility of the analytic manipulations is
supposed. This can be achieved by the assumption of su¬cient smoothness
for the corresponding functions and domains.
Let V ‚ „¦ be an admissible representative elementary volume in the
sense of volume averaging around a point x ∈ „¦. Typically the shape and
the size of a representative elementary volume are selected in such a manner
that the averaged values of all geometric characteristics of the microstruc-
ture of the pore space are independent of the size of V but depend on
the location of the point x. Then we obtain for a given variable ω± in the
phase ± (after continuation of ω± with 0 outside of ±) the corresponding
macroscopic quantities, assigned to the location x, as the extrinsic phase

ω± := ω±
|V | V

or as the intrinsic phase average

ω± := ω± .
|V± | V±

Here V± denotes the subset of V corresponding to ±. Let t ∈ (0, T ) be
the time at which the process is observed. The notation x ∈ „¦ means the
vector in Cartesian coordinates, whose coordinates are referred to by x,
y, and z ∈ R. Despite this ambiguity the meaning can always be clearly
derived from the context.
Let the index “s” (for solid) stand for the solid phase; then

φ(x) := |V \ Vs | |V | > 0

denotes the porosity, and for every liquid phase ±,

S± (x, t) := |V± | |V \ Vs | ≥ 0
0.3. Fluid Flow in Porous Media 7

is the saturation of the phase ±. Here we suppose that the solid phase is
stable and immobile. Thus
ω± = φS± ω±
for a ¬‚uid phase ± and

S± = 1 . (0.10)

So if the ¬‚uid phases are immiscible on the micro scale, they may be miscible
on the macro scale, and the immiscibility on the macro scale is an additional
assumption for the model.
As in other disciplines the di¬erential equation models are derived here
from conservation laws for the extensive quantities mass, impulse, and en-
ergy, supplemented by constitutive relationships, where we want to focus
on the mass.

0.3 Fluid Flow in Porous Media
Consider a liquid phase ± on the micro scale. In this chapter, for clarity, we
write “short” vectors in Rd also in bold with the exception of the coordinate
vector x. Let ˜± [kg/m3 ] be the (microscopic) density, q ± :=
˜ ˜
· ˜· v · ˜±
[m/s] the mass average mixture velocity based on the particle velocity v · of
a component · and its concentration in solution ˜· [kg/m ]. The transport
theorem of Reynolds (see, for example, [10]) leads to the mass conservation
‚t ˜± + ∇ · (˜± q ± ) = f±
˜ (0.11)
with a distributed mass source density f± . By averaging we obtain from
here the mass conservation law
‚t (φS± ±) ±q±) = f± (0.12)
with ± , the density of phase ±, as the intrinsic phase average of ˜± and
q ± , the volumetric ¬‚uid velocity or Darcy velocity of the phase ±, as the
extrinsic phase average of q ± . Correspondingly, f± is an average mass source
Before we proceed in the general discussion, we want to consider some
speci¬c situations: The area between the groundwater table and the imper-
meable body of an aquifer is characterized by the fact that the whole pore
space is occupied by a ¬‚uid phase, the soil water. The corresponding satu-
ration thus equals 1 everywhere, and with omission of the index equation
(0.12) takes the form
‚t (φ ) + ∇ · ( q) = f . (0.13)
8 0. Modelling Processes in Porous Media with Di¬erential Equations

If the density of water is assumed to be constant, due to neglecting
the mass of solutes and compressibility of water, equation (0.13) simpli¬es
further to the stationary equation
∇·q =f , (0.14)
where f has been replaced by the volume source density f / , keeping the
same notation. This equation will be completed by a relationship that
can be interpreted as the macroscopic analogue of the conservation of mo-
mentum, but should be accounted here only as an experimentally derived
constitutive relationship. This relationship is called Darcy™s law, which
reads as
q = ’K (∇p + gez ) (0.15)
and can be applied in the range of laminar ¬‚ow. Here p [N/m2 ] is the intrinsic
average of the water pressure, g [m/s2 ] the gravitational acceleration, ez the
unit vector in the z-direction oriented against the gravitation,
K = k/µ , (0.16)
a quantity, which is given by the permeability k determined by the solid
phase, and the viscosity µ determined by the ¬‚uid phase. For an anisotropic
solid, the matrix k = k(x) is a symmetric positive de¬nite matrix.
Inserting (0.15) in (0.14) and replacing K by K g, known as hydraulic
conductivity in the literature, and keeping the same notation gives the
following linear equation for
h(x, t) := p(x, t) + z ,
the piezometric head h [m]:
’∇ · (K∇h) = f . (0.17)
The resulting equation is stationary and linear. We call a di¬erential equa-
tion model stationary if it depends only on the location x and not on the
time t, and instationary otherwise. A di¬erential equation and correspond-
ing boundary conditions (cf. Section 0.5) are called linear if the sum or a
scalar multiple of a solution again forms a solution for the sum, respectively
the scalar multiple, of the sources.
If we deal with an isotropic solid matrix, we have K = KI with the d— d
unit matrix I and a scalar function K. Equation (0.17) in this case reads
’∇ · (K∇h) = f . (0.18)
Finally if the solid matrix is homogeneous, i.e., K is constant, we get from
division by K and maintaining the notation f the Poisson equation
’∆h = f , (0.19)
0.3. Fluid Flow in Porous Media 9

which is termed the Laplace equation for f = 0. This model and its more
general formulations occur in various contexts. If, contrary to the above as-
sumption, the solid matrix is compressible under the pressure of the water,
and if we suppose (0.13) to be valid, then we can establish a relationship
φ = φ(x, t) = φ0 (x)φf (p)
with φ0 (x) > 0 and a monotone increasing φf such that with S(p) := φf (p)
we get the equation
φ0 S(p) ‚t p + ∇ · q = f
and the instationary equations corresponding to (0.17)“(0.19), respectively.
For constant S(p) > 0 this yields the following linear equation:
φ0 S ‚t h ’ ∇ · (K∇h) = f , (0.20)
which also represents a common model in many contexts and is known from
corresponding ¬elds of application as the heat conduction equation.
We consider single phase ¬‚ow further, but now we will consider gas as
¬‚uid phase. Because of the compressibility, the density is a function of the
pressure, which is invertible due to its strict monotonicity to
p = P( ) .
Together with (0.13) and (0.15) we get a nonlinear variant of the heat
conduction equation in the unknown :
‚t (φ ) ’ ∇ · K( ∇P ( ) + 2
gez ) = f , (0.21)
which also contains derivatives of ¬rst order in space. If P ( ) = ln(± ) holds
for a constant ± > 0, then ∇P ( ) simpli¬es to ±∇ . Thus for horizontal
¬‚ow we again encounter the heat conduction equation. For the relationship
P ( ) = ± suggested by the universal gas law, ± ∇ = 1 ±∇ 2 remains
nonlinear. The choice of the variable u := 2 would result in u1/2 in the
time derivative as the only nonlinearity. Thus in the formulation in the
coe¬cient of ∇ disappears in the divergence of = 0. Correspondingly,
the coe¬cient S(u) = 1 φu’1/2 of ‚t u in the formulation in u becomes
unbounded for u = 0. In both versions the equations are degenerate, whose
treatment is beyond the scope of this book. A variant of this equation has
gained much attention as the porous medium equation (with convection) in
the ¬eld of analysis (see, for example, [42]).
Returning to the general framework, the following generalization of
Darcy™s law can be justi¬ed experimentally for several liquid phases:
q± = ’ k (∇p± + ± gez ) .
Here the relative permeability kr± of the phase ± depends upon the
saturations of the present phases and takes values in [0, 1].
10 0. Modelling Processes in Porous Media with Di¬erential Equations

At the interface of two liquid phases ±1 and ±2 we observe a di¬erence of
the pressures, the so-called capillary pressure, that turns out experimentally
to be a function of the saturations:

pc ±1 ±2 := p±1 ’ p±2 = F±1 ±2 (Sw , So , Sg ) . (0.22)

A general model for multiphase ¬‚ow, formulated for the moment in terms
of the variables p± , S± , is thus given by the equations

‚t (φS± ±) ± »± k(∇p± + ± gez )) = f± (0.23)

with the mobilities »± := kr± /µ± , and the equations (0.22) and (0.10),
where one of the S± ™s can be eliminated. For two liquid phases w and g,
e.g., water and air, equations (0.22) and (0.10) for ± = w, g read pc =
pg ’ pw = F (Sw ) and Sg = 1 ’ Sw . Apparently, this is a time-dependent,
nonlinear model in the variables pw , pg , Sw , where one of the variables can
be eliminated. Assuming constant densities ± , further formulations based

∇ · q w + q g = fw / + fg / (0.24)
w g

can be given as consequences of (0.10). These equations consist of a sta-
tionary equation for a new quantity, the global pressure, based on (0.24),
and a time-dependent equation for one of the saturations (see Exercise 0.2).
In many situations it is justi¬ed to assume a gaseous phase with constant
pressure in the whole domain and to scale this pressure to pg = 0. Thus
for ψ := pw = ’pc we have

φ‚t S(ψ) ’ ∇ · (»(ψ)k(∇ψ + gez )) = fw / (0.25)

:= w , and S(ψ) := F ’1 (’ψ) as a strictly
with constant pressure
monotone increasing nonlinearity as well as ».
With the convention to set the value of the air pressure to 0, the pressure
in the aqueous phase is in the unsaturated state, where the gaseous phase is
also present, and represented by negative values. The water pressure ψ = 0
marks the transition from the unsaturated to the saturated zone. Thus
in the unsaturated zone, equation (0.25) represents a nonlinear variant
of the heat conduction equation for ψ < 0, the Richards equation. As
most functional relationships have the property S (0) = 0, the equation
degenerates in the absence of a gaseous phase, namely to a stationary
equation in a way that is di¬erent from above.
Equation (0.25) with S(ψ) := 1 and »(ψ) := »(0) can be continued in a
consistent way with (0.14) and (0.15) also for ψ ≥ 0, i.e., for the case of a
sole aqueous phase. The resulting equation is also called Richards equation
or a model of saturated-unsaturated ¬‚ow.
0.4. Reactive Solute Transport in Porous Media 11

0.4 Reactive Solute Transport in Porous Media
In this chapter we will discuss the transport of a single component in a
liquid phase and some selected reactions. We will always refer to water
as liquid phase explicitly. Although we treat inhomogeneous reactions in
terms of surface reactions with the solid phase, we want to ignore exchange
processes between the ¬‚uid phases. On the microscopic scale the mass con-
servation law for a single component · is, in the notation of (0.11) by
omitting the phase index w,
‚t ˜· + ∇ · (˜· q ) + ∇ · J · = Q· ,
J · := ˜· (˜ · ’ q ) [kg/m2 /s]
v (0.26)
represents the di¬usive mass ¬‚ux of the component · and Q· [kg/m3 /s] is
its volumetric production rate. For a description of reactions via the mass
action law it is appropriate to choose the mole as the unit of mass. The
di¬usive mass ¬‚ux requires a phenomenological description. The assump-
tion that solely binary molecular di¬usion, described by Fick™s law, acts
between the component · and the solvent, means that
J · = ’ ˜D· ∇ (˜· / ˜) (0.27)
with a molecular di¬usivity D· > 0 [m2 /s]. The averaging procedure applied
on (0.26), (0.27) leads to
‚t (˜c· ) + ∇ · (qc· ) + ∇ · J (1) + ∇ · J (2) = Q(1) + Q(2)
· ·

for the solute concentration of the component ·, c· [kg/m3 ], as intrinsic
phase average of ˜· . Here, we have J (1) as the average of J · and J (2) ,
the mass ¬‚ux due to mechanical dispersion, a newly emerging term at the
(1) ˜
macroscopic scale. Analogously, Q· is the intrinsic phase average of Q· ,
and Q· is a newly emerging term describing the exchange between the
liquid and solid phases.
The volumetric water content is given by ˜ := φSw with the water
saturation Sw . Experimentally, the following phenomenological descriptions
are suggested:
J (1) = ’˜„ D· ∇c·
with a tortuosity factor „ ∈ (0, 1],
J (2) = ’˜Dmech ∇c· , (0.28)
and a symmetric positive de¬nite matrix of mechanical dispersion D mech ,
which depends on q/˜. Consequently, the resulting di¬erential equation
‚t (˜c· ) + ∇ · (qc· ’ ˜D∇c· ) = Q· (0.29)
12 0. Modelling Processes in Porous Media with Di¬erential Equations

(1) (2)
with D := „ D · + D mech, Q· := Q· + Q· .
Because the mass ¬‚ux consists of qc· , a part due to forced convection, and
of J (1) + J (2) , a part that corresponds to a generalized Fick™s law, an equa-
tion like (0.29) is called a convection-di¬usion equation. Accordingly, for
the part with ¬rst spatial derivatives like ∇ · (qc· ) the term convective part
is used, and for the part with second spatial derivatives like ’∇ · (˜D∇c· )
the term di¬usive part is used. If the ¬rst term determines the character of
the solution, the equation is called convection-dominated. The occurrence
of such a situation is measured by the quantity Pe, the global P´clet num-
ber, that has the form Pe = q L/ ˜D [ - ]. Here L is a characteristic
length of the domain „¦. The extreme case of purely convective transport
results in a conservation equation of ¬rst order. Since the common mod-
els for the dispersion matrix lead to a bound for Pe, the reduction to the
purely convective transport is not reasonable. However, we have to take
convection-dominated problems into consideration.
Likewise, we speak of di¬usive parts in (0.17) and (0.20) and of (nonlin-
ear) di¬usive and convective parts in (0.21) and (0.25). Also, the multiphase
transport equation can be formulated as a nonlinear convection-di¬usion
equation by use of (0.24) (see Exercise 0.2), where convection often dom-
inates. If the production rate Q· is independent of c· , equation (0.29) is
In general, in case of a surface reaction of the component ·, the kinetics of
the reaction have to be described . If this component is not in competition
with the other components, one speaks of adsorption. The kinetic equation
thus takes the general form
‚t s· (x, t) = k· f· (x, c· (x, t), s· (x, t)) (0.30)
with a rate parameter k· for the sorbed concentration s· [kg/kg], which is
given in reference to the mass of the solid matrix. Here, the components
in sorbed form are considered spatially immobile. The conservation of the
total mass of the component undergoing sorption gives
Q(2) = ’ b ‚t s· (0.31)

with the bulk density b = s (1’φ), where s denotes the density of the solid
phase. With (0.30), (0.31) we have a system consisting of an instationary
partial and an ordinary di¬erential equation (with x ∈ „¦ as parameter). A
widespread model by Langmuir reads
f· = ka c· (s· ’ s· ) ’ kd s·
with constants ka , kd that depend upon the temperature (among other
factors), and a saturation concentration s· (cf. for example [24]). If we
assume f· = f· (x, c· ) for simplicity, we get a scalar nonlinear equation in
c· ,
‚t (˜c· ) + ∇ · (qc· ’ ˜D∇c· ) + = Q(1) ,
b k· f· (·, c· ) (0.32)
0.4. Reactive Solute Transport in Porous Media 13

and s· is decoupled and extracted from (0.30). If the time scales of transport
and reaction di¬er greatly, and the limit case k· ’ ∞ is reasonable, then
(0.30) is replaced by
f· (x, c· (x, t), s· (x, t)) = 0 .
If this equation is solvable for s· , i.e.,
s· (x, t) = •· (x, c· (x, t)) ,
the following scalar equation for c· with a nonlinearity in the time
derivative emerges:
+ ∇ · (qc· ’ ˜D∇c· ) = Q(1) .
‚t (˜c· + b •· (·, c· )) ·

If the component · is in competition with other components in the sur-
face reaction, as, e.g., in ion exchange, then f· has to be replaced by a
nonlinearity that depends on the concentrations of all involved components
c1 , . . . , cN , s1 , . . . , sN . Thus we obtain a coupled system in these variables.
Finally, if we encounter homogeneous reactions that take place solely in the
¬‚uid phase, an analogous statement is true for the source term Q· .

0.1 Give a geometric interpretation for the matrix condition of k in (0.16)
and D mech in (0.28).

± ∈ {w, g})
0.2 Consider the two-phase ¬‚ow (with constant ±,

‚t (φS± ) + ∇ · q ± = f± ,
’»± k (∇p± +
q± = ± gez ) ,
Sw + Sg = 1,
pg ’ pw = pc
with coe¬cient functions
± ∈ {w, g}.
pc = pc (Sw ) , »± = »± (Sw ) ,
Starting from equation (0.23), perform a transformation to the new
q = qw + qg , “total ¬‚ow,”

»g ’ »w dpc
1 1
p= (pw + pg ) + dξ , “global pressure,”
2 2 »g + »w dξ

and the water saturation Sw . Derive a representation of the phase ¬‚ows in
the new variables.
14 0. Modelling Processes in Porous Media with Di¬erential Equations

0.3 A frequently employed model for mechanical dispersion is
Dmech = »L |v|2 Pv + »T |v|2 (I ’ Pv )
with parameters »L > »T , where v = q/˜ and Pv = vv T /|v|2 . Here2
»L and »T are the longitudinal and transversal dispersion lengths. Give a
geometrical interpretation.

0.5 Boundary and Initial Value Problems
The di¬erential equations that we derived in Sections 0.3 and 0.4 have the
common form
‚t S(u) + ∇ · (C(u) ’ K(∇u)) = Q(u) (0.33)
with a source term S, a convective part C, a di¬usive part K, i.e., a total
¬‚ux C ’ K and a source term Q, which depend linearly or nonlinearly
on the unknown u. For simpli¬cation, we assume u to be a scalar. The
nonlinearities S, C, K, and Q may also depend on x and t, which shall be
suppressed in the notation in the following. Such an equation is said to be
in divergence form or in conservative form; a more general formulation is
obtained by di¬erentiating ∇ · C(u) = ‚u C(u) · ∇u + (∇ · C)(u) or by

introducing a generalized “source term” Q = Q(u, ∇u). Up to now we have
considered di¬erential equations pointwise in x ∈ „¦ (and t ∈ (0, T )) under
the assumption that all occurring functions are well-de¬ned. Due to the
applicability of the integral theorem of Gauss on „¦ ‚ „¦ (cf. (3.10)), the
integral form of the conservation equation follows straightforwardly from
the above:

(C(u) ’ K(∇u)) · ν dσ = Q(u, ∇u) dx
‚t S(u) dx + (0.34)
˜ ˜
„¦ „¦

with the outer unit normal ν (see Theorem 3.8) for a ¬xed time t or also
in t integrated over (0, T ). Indeed, this equation (on the microscopic scale)
is the primary description of the conservation of an extensive quantity:
Changes in time through storage and sources in „¦ are compensated by the
normal ¬‚ux over ‚ „¦. Moreover, for ‚t S, ∇ · (C ’ K), and Q continuous
on the closure of „¦, (0.33) follows from (0.34). If, on the other hand, F is
a hyperplane in „¦ where the material properties may rapidly change, the
jump condition
[(C(u) ’ K(∇u)) · ν] = 0 (0.35)
for a ¬xed unit normal ν on F follows from (0.34), where [ · ] denotes the
di¬erence of the one-sided limits (see Exercise 0.4).
Since the di¬erential equation describes conservation only in general,
it has to be supplemented by initial and boundary conditions in order to
0.5. Boundary and Initial Value Problems 15

specify a particular situation where a unique solution is expected. Boundary
conditions are speci¬cations on ‚„¦, where ν denotes the outer unit normal
• of the normal component of the ¬‚ux (inwards):
’ (C(u) ’ K(∇u)) · ν = g1 on “1 (0.36)
(¬‚ux boundary condition),
• of a linear combination of the normal ¬‚ux and the unknown itself:
’ (C(u) ’ K(∇u)) · ν + ±u = g2 on “2 (0.37)
(mixed boundary condition),
• of the unknown itself:
u = g3 on “3 (0.38)
(Dirichlet boundary condition).
Here “1 , “2 , “3 form a disjoint decomposition of ‚„¦:
‚„¦ = “1 ∪ “2 ∪ “3 , (0.39)
where “3 is supposed to be a closed subset of ‚„¦. The inhomogeneities
gi and the factor ± in general depend on x ∈ „¦, and for nonstationary
problems (where S(u) = 0 holds) on t ∈ (0, T ). The boundary conditions
are linear if the gi do not depend (nonlinearly) on u (see below). If the gi
are zero, we speak of homogeneous, otherwise of inhomogeneous, boundary
Thus the pointwise formulation of a nonstationary equation (where S
does not vanish) requires the validity of the equation in the space-time
QT := „¦ — (0, T )
and the boundary conditions on the lateral surface of the space-time
ST := ‚„¦ — (0, T ) .
Di¬erent types of boundary conditions are possible with decompositions
of the type (0.39). Additionally, an initial condition on the bottom of the
space-time cylinder is necessary:
for x ∈ „¦ .
S(u(x, 0)) = S0 (x) (0.40)
These are so-called initial-boundary value problems; for stationary prob-
lems we speak of boundary value problems. As shown in (0.34) and (0.35)
¬‚ux boundary conditions have a natural relationship with the di¬erential
equation (0.33). For a linear di¬usive part K(∇u) = K∇u alternatively
we may require
‚νK u := K∇u · ν = g1 on “1 , (0.41)
16 0. Modelling Processes in Porous Media with Di¬erential Equations

and an analogous mixed boundary condition. This boundary condition is
the so-called Neumann boundary condition. Since K is symmetric, ‚νK u =
∇u · Kν holds; i.e., ‚νK u is the derivative in direction of the conormal Kν.
For the special case K = I the normal derivative is given.
In contrast to ordinary di¬erential equations, there is hardly any general
theory of partial di¬erential equations. In fact, we have to distinguish dif-
ferent types of di¬erential equations according to the various described
physical phenomena. These determine, as discussed, di¬erent (initial-)
boundary value speci¬cations to render the problem well-posed. Well-
posedness means that the problem possesses a unique solution (with certain
properties yet to be de¬ned) that depends continuously (in appropriate
norms) on the data of the problem, in particular on the (initial and)
boundary values. There exist also ill-posed boundary value problems for
partial di¬erential equations, which correspond to physical and technical
applications. They require special techniques and shall not be treated here.
The classi¬cation into di¬erent types is simple if the problem is lin-
ear and the di¬erential equation is of second order as in (0.33). By order
we mean the highest order of the derivative with respect to the variables
(x1 , . . . , xd , t) that appears, where the time derivative is considered to be
like a spatial derivative. Almost all di¬erential equations treated in this
book will be of second order, although important models in elasticity the-
ory are of fourth order or certain transport phenomena are modelled by
systems of ¬rst order.
The di¬erential equation (0.33) is generally nonlinear due to the nonlin-
ear relationships S, C, K, and Q. Such an equation is called quasilinear if
all derivatives of the highest order are linear, i.e., we have
K(∇u) = K∇u (0.42)
with a matrix K, which may also depend (nonlinearly) on x, t, and u.
Furthermore, (0.33) is called semilinear if nonlinearities are present only
in u, but not in the derivatives, i.e., if in addition to (0.42) with K being
independent of u, we have
S(u) = Su , C(u) = uc (0.43)
with scalar and vectorial functions S and c, respectively, which may depend
on x and t. Such variable factors standing before u or di¬erential terms are
called coe¬cients in general.
Finally, the di¬erential equation is linear if we have, in addition to the
above requirements,
Q(u) = ’ru + f
with functions r and f of x and t.
In the case f = 0 the linear di¬erential equation is termed homoge-
neous, otherwise inhomogeneous. A linear di¬erential equation obeys the
superposition principle: Suppose u1 and u2 are solutions of (0.33) with the
0.5. Boundary and Initial Value Problems 17

source terms f1 and f2 and otherwise identical coe¬cient functions. Then
u1 + γu2 is a solution of the same di¬erential equation with the source
term f1 + γf2 for arbitrary γ ∈ R. The same holds for linear boundary
conditions. The term solution of an (initial-) boundary value problem is
used here in a classical sense, yet to be speci¬ed, where all the quantities
occurring should satisfy pointwise certain regularity conditions (see De¬ni-
tion 1.1 for the Poisson equation). However, for variational solutions (see
De¬nition 2.2), which are appropriate in the framework of ¬nite element
methods, the above statements are also valid.
Linear di¬erential equations of second order in two variables (x, y) (in-
cluding possibly the time variable) can be classi¬ed in di¬erent types as
To the homogeneous di¬erential equation
‚2 ‚2 ‚2
Lu = a(x, y) 2 u + b(x, y) u + c(x, y) 2 u
‚x ‚x‚y ‚y (0.44)
‚ ‚
+ d(x, y) u + e(x, y) u + f (x, y)u = 0
‚x ‚y
the following quadratic form is assigned:
(ξ, ·) ’ a(x, y)ξ 2 + b(x, y)ξ· + c(x, y)· 2 . (0.45)
According to its eigenvalues, i.e., the eigenvalues of the matrix
a(x, y) 2 b(x, y)
, (0.46)
2 b(x, y) c(x, y)
we classify the types. In analogy with the classi¬cation of conic sections,
which are described by (0.45) (for ¬xed (x, y)), the di¬erential equation
(0.44) is called at the point (x, y)
• elliptic if the eigenvalues of (0.46) are not 0 and have the same sign,
• hyperbolic if one eigenvalue is positive and the other is negative,
• parabolic if exactly one eigenvalue is equal to 0.
For the corresponding generalization of the terms for d + 1 variables and
arbitrary order, the stationary boundary value problems we treat in this
book will be elliptic, of second order, and ” except in Chapter 8 ” also
linear; the nonstationary initial-boundary value problems will be parabolic.
Systems of hyperbolic di¬erential equations of ¬rst order require partic-
ular approaches, which are beyond the scope of this book. Nevertheless,
we dedicate Chapter 9 to convection-dominated problems, i.e., elliptic or
parabolic problems close to the hyperbolic limit case.
The di¬erent discretization strategies are based on various formulations
of the (initial-) boundary value problems: The ¬nite di¬erence method,
which is presented in Section 1, and further outlined for nonstationary prob-
lems in Chapter 7, has the pointwise formulation of (0.33), (0.36)“(0.38)
18 0. Modelling Processes in Porous Media with Di¬erential Equations

(and (0.40)) as a starting point. The ¬nite element method, , which lies in
the focus of our book (Chapters 2, 3, and 7), is based on an integral formu-
lation of (0.33) (which we still have to depict) that incorporates (0.36) and
(0.37). The conditions (0.38) and (0.40) have to be enforced additionally.
Finally, the ¬nite volume method (Chapters 6 and 7) will be derived from
the integral formulation (0.34), where also initial and boundary conditions
come along as in the ¬nite element approach.

0.4 Derive (formally) (0.35) from (0.34).

0.5 Derive the orders of the given di¬erential operators and di¬er-
ential equations, and decide in every case whether the operator is
linear or nonlinear, and whether the linear equation is homogeneous or
(a) Lu := uxx + xuy ,
(b) Lu := ux + uuy ,

(c) Lu := 1 + x2 (cos y)ux + uyxy ’ arctan x u = ln(x2 + y 2 ) ,

(d) Lu := ut + uxxxx + 1 + u = 0 ,
(e) utt ’ uxx + x2 = 0 .

0.6 (a) Determine the type of the given di¬erential operator:
(i) Lu := uxx ’ uxy + 2uy + uyy ’ 3uyx + 4u ,
(ii) Lu = 9uxx + 6uxy + uyy + ux .
(b) Determine the parts of the plane where the di¬erential operator Lu :=
yuxx ’ 2uxy + xuyy is elliptic, hyperbolic, or parabolic.
(c) (i) Determine the type of Lu := 3uy + uxy .
(ii) Compute the general solution of Lu = 0.

0.7 Consider the equation Lu = f with a linear di¬erential operator of
second order, de¬ned for functions in d variables (d ∈ N) in x ∈ „¦ ‚ Rd .
The transformation ¦ : „¦ ’ „¦ ‚ Rd has a continuously di¬erentiable,
nonsingular Jacobi matrix D¦ := ‚¦ .
Show that the partial di¬erential equation does not change its type if it
is written in the new coordinates ξ = ¦(x).
For the Beginning:
The Finite Di¬erence Method for the
Poisson Equation

1.1 The Dirichlet Problem for the Poisson
In this section we want to introduce the ¬nite di¬erence method using
the Poisson equation on a rectangle as an example. By means of this ex-
ample and generalizations of the problem, advantages and limitations of
the approach will be elucidated. Also, in the following section the Poisson
equation will be the main topic, but then on an arbitrary domain. For the
spatial basic set of the di¬erential equation „¦ ‚ Rd we assume as minimal
requirement that „¦ is a domain, where a domain is a nonempty, open, and
connected set. The boundary of this domain will be denoted by ‚„¦, the
closure „¦ ∪ ‚„¦ by „¦ (see Appendix A.2). The Dirichlet problem for the
Poisson equation is then de¬ned as follows: Given functions g : ‚„¦ ’ R
and f : „¦ ’ R, we are looking for a function u : „¦ ’ R such that
’ u = f in „¦ , (1.1)
u = g on ‚„¦ . (1.2)

This di¬erential equation model has already appeared in (0.19) and
(0.38) and beyond this application has an importance in a wide spectrum of
disciplines. The unknown function u can be interpreted as an electromag-
netic potential, a displacement of an elastic membrane, or a temperature.
Similar to the multi-index notation to be introduced in (2.16) (but with
20 1. Finite Di¬erence Method for the Poisson Equation

indices at the top) from now on for partial derivatives we use the following

Notation: For u : „¦ ‚ Rd ’ R we set

‚i u := u for i = 1, . . . , d ,
for i, j = 1, · · · , d ,
‚ij u := u
‚xi ‚xj
∆u := (‚11 + . . . + ‚dd ) u .
The expression ∆u is called the Laplace operator. By means of this, (1.1)
can be written in abbreviated form as
’∆u = f in „¦ . (1.3)
We could also de¬ne the Laplace operator by
∆u = ∇ · (∇u) ,
where ∇u = (‚1 u, . . . , ‚d u) denotes the gradient of a function u, and
∇ · v = ‚1 v1 + · · · + ‚d vd the divergence of a vector ¬eld v. Therefore,
an alternative notation exists, which will not be used in the following:
∆u = ∇2 u. The incorporation of the minus sign in the left-hand side of
(1.3), which looks strange at ¬rst glance, is related to the monotonicity and
de¬niteness properties of ’∆ (see Sections 1.4 and 2.1, respectively).
The notion of a solution for (1.1), (1.2) still has to speci¬ed more pre-
cisely. Considering the equations in a pointwise sense, which will be pursued
in this chapter, the functions in (1.1), (1.2) have to exist, and the equations
have to be satis¬ed pointwise. Since (1.1) is an equation on an open set „¦,
there are no implications for the behaviour of u up to the boundary ‚„¦. To
have a real requirement due to the boundary condition, u has to be at least
continuous up to the boundary, that is, on „¦. These requirements can be
formulated in a compact way by means of corresponding function spaces.
The function spaces are introduced more precisely in Appendix A.5. Some
examples are
u : „¦ ’ R u continuous in „¦ ,
C(„¦) :=
u : „¦ ’ R u ∈ C(„¦) , ‚i u exists in „¦ ,
C 1 („¦) :=
‚i u ∈ C(„¦) for all i = 1, . . . , d .
The spaces C k („¦) for k ∈ N, C(„¦), and C k („¦), as well as C(‚„¦), are
de¬ned analogously. In general, the requirements related to the (contin-
uous) existence of derivatives are called, a little bit vaguely, smoothness
In the following, in view of the ¬nite di¬erence method, f and g will also
be assumed continuous in „¦ and ‚„¦, respectively.
1.2. The Finite Di¬erence Method 21

De¬nition 1.1 Assume f ∈ C(„¦) and g ∈ C(‚„¦). A function u is called
a (classical) solution of (1.1), (1.2) if u ∈ C 2 („¦) © C(„¦), (1.1) holds for all
x ∈ „¦, and (1.2) holds for all x ∈ ‚„¦.

1.2 The Finite Di¬erence Method
The ¬nite di¬erence method is based on the following approach: We are
looking for an approximation to the solution of a boundary value problem
at a ¬nite number of points in „¦ (the grid points). For this reason we
substitute the derivatives in (1.1) by di¬erence quotients, which involve
only function values at grid points in „¦ and require (1.2) only at grid
points. By this we obtain algebraic equations for the approximating values
at grid points. In general, such a procedure is called the discretization of the
boundary value problem. Since the boundary value problem is linear, the
system of equations for the approximate values is also linear. In general, for
other (di¬erential equation) problems and other discretization approaches
we also speak of the discrete problem as an approximation of the continuous
problem. The aim of further investigations will be to estimate the resulting
error and thus to judge the quality of the approximative solution.

Generation of Grid Points
In the following, for the beginning, we will restrict our attention to problems
in two space dimensions (d = 2). For simpli¬cation we consider the case
of a constant step size (or mesh width) h > 0 in both space directions.
The quantity h here is the discretization parameter, which in particular
determines the dimension of the discrete problem.

—¦—¦—¦ —¦ —¦ —¦ —¦ —¦ —¦ • : „¦h
—¦ •• • • • • • —¦
—¦ : ‚„¦h
—¦ •• • • • • • —¦
l=8 2
—¦ •• • • • • • —¦
2 : far from boundary
—¦ •• • • • • • —¦
—¦—¦—¦ —¦ —¦ —¦ —¦ —¦ —¦ 3 : close to boundary

Figure 1.1. Grid points in a square domain.

For the time being, let „¦ be a rectangle, which represents the simplest
case for the ¬nite di¬erence method (see Figure 1.1). By translation of the
coordinate system the situation can be reduced to „¦ = (0, a) — (0, b) with
a, b > 0 . We assume that the lengths a, b, and h are such that
for certain l, m ∈ N.
a = lh, b = mh (1.4)
22 1. Finite Di¬erence Method for the Poisson Equation

We de¬ne
(ih, jh) i = 1, . . . , l ’ 1 , j = 1, . . . , m ’ 1
„¦h :=
(x, y) ∈ „¦ x = ih , y = jh with i, j ∈ Z
as a set of grid points in „¦ in which an approximation of the di¬erential
equation has to be satis¬ed. In the same way,
(ih, jh) i ∈ {0, l} , j ∈ {0, . . . , m} or i ∈ {0, . . . , l} , j ∈ {0, m}
‚„¦h :=
(x, y) ∈ ‚„¦ x = ih , y = jh with i, j ∈ Z
de¬nes the grid points on ‚„¦ in which an approximation of the boundary
condition has to be satis¬ed. The union of grid points will be denoted by
„¦h := „¦h ∪ ‚„¦h .

Setup of the System of Equations
Lemma 1.2 Let „¦ := (x ’ h, x + h) for x ∈ R, h > 0. Then there exists
a quantity R, depending on u and h, the absolute value of which can be
bounded independently of h and such that
(1) for u ∈ C 2 („¦):
u(x + h) ’ u(x) 1
|R| ¤
u (x) = + hR and u ,

h 2
(2) for u ∈ C 2 („¦):
u(x) ’ u(x ’ h) 1
|R| ¤
u (x) = + hR and u ,

h 2
(3) for u ∈ C 3 („¦):
u(x + h) ’ u(x ’ h) 1
|R| ¤
+ h2 R
u (x) = and u ,

2h 6
(4) for u ∈ C 4 („¦):
u(x + h) ’ 2u(x) + u(x ’ h) 1 (4)
+ h2 R and |R| ¤
u (x) = u .

h 12
Here the maximum norm · ∞ (see Appendix A.5) has to be taken over
the interval of the involved points (x, x + h), (x ’ h, x), or (x ’ h, x + h).

Proof: The proof follows immediately by Taylor expansion. As an example
we consider statement 3: From
h2 h3
u(x ± h) = u(x) ± hu (x) + u (x) ± u (x ± ξ± ) for certain ξ± ∈ (0, h)
2 6
the assertion follows by linear combination.
1.2. Derivation and Properties 23

Notation: The quotient in statement 1 is called the forward di¬erence
quotient, and it is denoted by ‚ + u(x). The quotient in statement 2 is
called the backward di¬erence quotient (‚ ’ u(x)), and the one in statement
3 the symmetric di¬erence quotient (‚ 0 u(x)). The quotient appearing in
statement 4 can be written as ‚ ’ ‚ + u(x) by means of the above notation.
In order to use statement 4 in every space direction for the approximation
of ‚11 u and ‚22 u in a grid point (ih, jh), in addition to the conditions of
De¬nition 1.1, the further smoothness properties ‚ (3,0) u, ‚ (4,0) u ∈ C(„¦)
and analogously for the second coordinate are necessary. Here we use, e.g.,
the notation ‚ (3,0) u := ‚ 3 u/‚x3 (see (2.16)).
Using these approximations for the boundary value problem (1.1), (1.2),
at each grid point (ih, jh) ∈ „¦h we get
u ((i + 1)h, jh) ’ 2u(ih, jh) + u ((i ’ 1)h, jh)

u (ih, (j + 1)h) ’ 2u(ih, jh) + u (ih, (j ’ 1)h)
+ = (1.6)
f (ih, jh) + R(ih, jh)h2 .
Here R is as described in statement 4 of Lemma 1.2, a function depending
on the solution u and on the step size h, but the absolute value of which can
be bounded independently of h. In cases where we have less smoothness of
the solution u, we can nevertheless formulate the approximation (1.6) for
’∆u, but the size of the error in the equation is unclear at the moment.
For the grid points (ih, jh) ∈ ‚„¦h no approximation of the boundary
condition is necessary:
u(ih, jh) = g(ih, jh) .
If we neglect the term Rh2 in (1.6), we get a system of linear equations
for the approximating values uij for u(x, y) at points (x, y) = (ih, jh) ∈ „¦h .
They have the form
’ ui,j’1 ’ ui’1,j + 4uij ’ ui+1,j ’ ui,j+1 = fij (1.7)
for i = 1, . . . , l ’ 1 , j = 1, . . . , m ’ 1 ,
if i ∈ {0, l}, j = 0, . . . , m or j ∈ {0, m}, i = 0, . . . , l .
uij = gij (1.8)
Here we used the abbreviations
fij := f (ih, jh), gij := g(ih, jh) . (1.9)
Therefore, for each unknown grid value uij we get an equation. The grid
points (ih, jh) and the approximating values uij located at these have a
natural two-dimensional indexing.
In equation (1.7) for a grid point (i, j) only the neighbours at the four
cardinal points of the compass appear, as it is displayed in Figure 1.2. This
24 1. Finite Di¬erence Method for the Poisson Equation

interconnection is also called the ¬ve-point stencil of the di¬erence method
and the method the ¬ve-point stencil discretization.
T (i,j+1)

(i’1,j) (i,j) (i+1,j)
• • •


Figure 1.2. Five-point stencil.

At the interior grid points (x, y) = (ih, jh) ∈ „¦h , two cases can be
(1) (i, j) has a position such that its all neighbouring grid points lie in
„¦h (far from the boundary).
(2) (i, j) has a position such that at least one neighbouring grid point
(r, s) lies on ‚„¦h (close to the boundary). Then in equation (1.7) the
value urs is known due to (1.8) (urs = grs ), and (1.7) can be modi¬ed
in the following way:
Remove the values urs with (rh, sh) ∈ ‚„¦h in the equations for (i, j)
close to the boundary and add the value grs /h2 to the right-hand
side of (1.7). The set of equations that arises by this elimination of
boundary unknowns by means of Dirichlet boundary conditions we
call (1.7)— ; it is equivalent to (1.7), (1.8).
Instead of considering the values uij , i = 1, . . . , l ’ 1, j = 1, . . . , m ’ 1,

. 1
( 16)