Principal Encyclopedia of Computational Mechanics (3 Volume Set)

Encyclopedia of Computational Mechanics (3 Volume Set)

, , ,
Over the past four decades computational methods in applied mechanics have developed into valuable tools that are widely used across both industry and academia. The applications are numerous: aerospace structures, civil engineering structures, geotechnics, flow problems, automotive industry, geo-environmental modelling, biomechanics, electromagnetism, metal forming, to name but a few.This three volume set provides the most comprehensive and up-to-date collection of knowledge about this increasingly important area of engineering science.  The Encyclopedia provides a well-rounded and practical knowledge base that will serve as a foundation for the reader's research and practice in developing designs and in understanding, assessing and managing numerical analysis systems.Containing over 70 in-depth and thoroughly cross referenced articles on key topics from internationally renowned researchers, the Encyclopedia of Computational Mechanics will cover three key areas.Volume One: Fundamentals will cover the basic concepts behind discretization, interpolation, error estimation, solvers, computer algebra and geometric modelling.Volume Two: Solids and Volume Three: Fluids will build on this foundation with extensive, in-depth coverage of industrial applications. The main readership for this book will be researchers, research students (PhD. D. and postgraduate) and professional engineers in industrial and governmental laboratories.  Academic interest will stem from civil, mechanical, geomechanical, biomedical, aerospace and chemical engineering departments, through to the fields of applied mathematics, computer science and physics.
ISBN 10:
ISBN 13:
PDF, 57.87 MB
Descarga (pdf, 57.87 MB)

You may be interested in


Most frequently terms

You can write a book review and share your experiences. Other readers will always be interested in your opinion of the books you've read. Whether you've loved the book or not, if you give your honest and detailed thoughts then people will find new books that are right for them.
Chapter 1
Fundamentals: Introduction and Survey
Erwin Stein
University of Hannover, Hannover, Germany

1 Motivation and Scope
2 Stages of Development and Features of
Computational Mechanics
3 Survey of the Chapters of Volume 1
4 What We Do Expect


In the ‘Encyclopedia of Computational Mechanics’ (ECM),
Volume 1 ‘Fundamentals’ includes 26 chapters. It contains the basic methodological, analytical, algorithmic, and
implementation topics of computational mechanics.
The main goals of the ECM are to provide first-class
up-to-date representations of all major computer-oriented
numerical methods and related special features for mechanical problems in space and time, their a priori and a posteriori error analysis as well as various convergent and efficient
self-controlling adaptive discretization strategies and to further provide the wide range of robust and efficient direct and
iterative solvers as well as challenging applications in all
relevant technological areas. Geometrical representations of
technical objects, mesh generations, and mesh adaptivity as
well as the visualization of input- and output data are also
important topics.
The now already ‘classical’ discretization methods using
finite differences, finite elements, finite volumes, and
boundary elements were generalized with new conceptual
ideas into various directions in the last decade, such as
Encyclopedia of Computational Mechanics, Edited by Erwin
Stein, René de Borst and Thomas J.R. Hughes. Volume 1: Fundamentals.  2004 John Wiley & Sons, Ltd. ISBN: 0-470-84699-2.

meshfree methods, spectral, and wavelet techniques as well
as discrete finite element algorithms, which are presented
here. Error analysis and adaptivity are essential features
in general.

One can say that we are now in about the third period
of development of computer-based numerical methods,
especially based on weighted residuals, such as the finite
element method (FEM) and its genera; lizations as well as
the boundary integral equation method (BIEM or simply
BEM) and various couplings of both. The finite difference
method (FDM) further plays an important role, especially
for time integrations.
The first period was from about 1960 to 1975, with
a lot of separate engineering approximations for specific
mathematical models, especially FEMs for linear elastic
static systems, like beams and plates with plane stress
states and bending, as well as eigenvalue analysis of stability and vibration problems with applications to structural mechanics. Parallel to this, FEMs for aero- and
hydrostatic problems and also for hydrodynamic processes
were developed.
The second period from about 1976 to 1990 was characterized by rigorous mathematical analysis of Ritz–Galerkintype discretization methods with trial and test functions in
Sobolev spaces within finite subdomains, in order to analyze elliptic boundary-value problems, together with the a
priori and a posteriori error analysis of the FEM and the
BEM in their various forms for large classes of problems,


Fundamentals: Introduction and Survey

such as boundary-value problems of symmetric, positive
definite elliptic operators of second order, and also for
parabolic and hyperbolic operators with operator splits in
space and time, solving systems of ordinary differential
equations in time by a finite difference method. Parallel
to this, sophisticated engineering developments took place
toward complicated linear and nonlinear problems of the
classical partial differential equations (PDEs) in mathematical physics with large dimensions of the algebraic equations, motivated and driven by the fast growth of computer
power and memory as well as the availability of efficient
software systems and, of course, by technological needs and
motivations. Numerous chapters in Volumes 2 and 3 show
the wide range of challenging applications of FEM, BEM,
and various other problem-oriented discretization methods.
These new integral methods of weighted residuals are characterized by two important properties:






Methodical width: This is the intended simple logical and algorithmic structure (e.g. with symmetry properties and well-posedness in regular cases)
and the possibility of extensions and generalizations
within a large class of similar problems, including
higher dimensions, and thus forming the frame for a
box-type program structure. This is mainly achieved
by operations within the finite element subdomains
only, that is, without needing neighborhood information on element level, and thus allowing unified
assembling and solution procedures of the global
systems of algebraic equations. The methods yield
relatively small condition numbers of the algebraic
equation systems and thus provide robust solutions in
regular cases.
Methodical depth: This means the rather simple
extension of methods, algorithms, and computer programs to more complicated – especially geometrically
and physically nonlinear – problems and to physically
coupled problems. This also holds for the implementation of sensitivity analysis within the solution of
optimization problems.

These two properties (a) and (b) are the reasons for the
tremendous development in and the flexible availability of
related program systems for applications in science and
In the third period, from 1991 until now, new tasks,
challenges, and research directions can be observed in computational mechanics (more general in applied physics) and
in computational mathematics that can be summarized as

Meshfree and particle methods, finite elements with
discontinuities for damage and fracture.

Error-controlled adaptive modeling and approximation
of physical events near to nature, also scale-bridging
modeling on different space and timescales, including
homogenizations between them.
Adaptive micromechanical modeling and computation
in material science and engineering, including damage,
phase changes, and various failure processes.
New types of generalized FEM and BEM with hierarchical, spectral, and wavelet-based interpolations.
Modeling and simulation of multiphysics phenomena
in science and engineering.
Complex models and simulations in biomechanics and
human medicine.
New generalized methods for geometrical modeling,
mesh generation, and mesh adaptivity.
New direct and iterative solvers with multilevel and
domain decomposition methods.
Advanced visualization of objects, processes, and
numerical results, 3D-animation of virtual reality.

With the advanced current hardware and software tools
on hand, about 10 million unknowns of a complex problem can be computed today in reasonable time, using
problem-oriented iterative algebraic solvers and preconditioners with advanced data management for high-end
machines with parallel or serial scalar and vector processors. Personal computers also enable us to solve hundreds
of thousands of unknowns together with error estimation
and adaptive mesh refinements. With these tools, it has
become possible to realize the verification and even the
restricted validation of engineering systems and processes,
taking into account disturbed input data and deterministic or statistic imperfections of structures and materials.
This leads us to new paradigms in computational mechanics, namely, guaranteeing reliability, safety, and efficiency
of results very near to the physical reality of the investigated objects. And because of this progress, computational
mechanics helps to simulate virtual products and processes
without the necessity of many physical experiments and
thus reduces costs and development time of new products

Volume 1 can be classified into the four groups: discretization methods (Chapter 2 to Chapter 15 of this Volume
(14 chapters)); geometrical modeling, mesh generation, and
visualization (Chapter 16 to Chapter 18 of this Volume (3 chapters)); Solvers (Chapter 19 to Chapter 23 of
this Volume (5 chapters)); and time-dependent problems
(Chapter 24 to Chapter 26 of this Volume (3 chapters)).

Fundamentals: Introduction and Survey 3
The first group, discretization methods, begins with
Finite difference methods by Owe Axelsson in which
elliptic, parabolic, and hyperbolic problems of second and
fourth order as well as convection–diffusion problems are
treated in a systematic way, including error analysis and
adaptivity, emphasizing computational issues.
Next, FEMs are presented in six chapters, beginning
with Interpolation in finite element spaces by Thomas
Apel, with a survey of different types of test and trial
functions, investigating the interpolation error as a basis for
a priori and a posteriori error estimates of finite element
methods. For a priori estimates, nodal interpolants are
used as well as the maximum available regularity of the
solution to get optimal error bounds. A posteriori error
estimates of the residual type need local interpolation
error representations for functions from the Sobolev space
W1,2 (). Different interpolation operators and related error
estimates are presented for the h-version of the usually used
2D- and 3D-finite elements.
The following chapter, Finite element methods, by
Susanne Brenner and Carsten Carstensen, treats the displacement method (primal finite element method) for
boundary-value problems of second-order elliptic PDE’s
as well as a priori and a posteriori error estimates of
the weak solutions and related h-adaptivity, including nonconforming elements and algorithmic aspects. This basic
chapter is followed by The p-version of the finite element
method for elliptic problems, by Barna Szabó, Alexander Düster, and Ernst Rank, in which hierarchical shape
functions of order p are used as test and trial interpolations of the finite elements instead of nodal basis functions.
Exponential convergence rates in conjunction with sufficient h-refinement in subdomains with large gradients of
the solution are advantageous against the h-version. Boundary layers of dimensionally reduced models (by appropriate
kinematic hypotheses), which need the solutions of the
expanded mathematical model, can be represented in a
consistent way by using adequate p-orders, see also Chapter 8, this Volume, by Monique Dauge et al. The arising
problems are: (i) the fast integration of fully populated
element stiffness matrices, (ii) relatively large algebraic
systems with strongly populated global stiffness matrices,
(iii) the problem of geometric boundary representations
without producing artificial singularities, (iv) hp-adaptivity
for 3D-systems as well as (v) the efficient implementation
of anisotropic p-extensions that are efficient for geometrically and physically anisotropic problems like thin plates
and shells, for example, with anisotropic layers of composites. All these problems have been tackled successfully
such that the p-type finite element method – in connection
with some available computer programs – has reached the
necessary maturity for engineering practice.

In Chapter 6 to Chapter 8 of this Volume, problemoriented effective test and trial spaces are introduced for
BVPs of PDEs.
Chapter 6, this Volume, by Claudio Canuto and Alfio
Quarteroni, is devoted to the high-order trigonometric and
orthogonal Jacobi polynomial expansions to be applied to
generalized Galerkin methods for periodic and nonperiodic
problems, with numerical integration via Gaussian integration points in order to achieve high rates of convergence
in total.
Chapter 7, this Volume, by Albert Cohen, Wolfgang
Dahmen, and Ronald De Vore, represents matrix compression methods for the BIEM based on wavelet coordinates
with application to time-dependent and stationary problems.
Wavelets also yield sparsity for the conserved variables of
problems with hyperbolic conservation laws. In addition, a
new adaptive algorithm is derived for sparse functions and
operators of linear and nonlinear problems.
Chapter 8, this Volume, by Monique Dauge, Erwan
Faou, and Zohar Yosibash, treats known and new methods for consistent reductions of the 3-D theory of elasticity
to 2-D theories of thin-walled plates and shells by expansions with respect to small parameters, without applying
the traditional kinematic and static hypotheses. A polynomial representation of the displacements is presumed,
depending on the thickness direction, generating singularly
perturbed boundary layers in the zero thickness limit. This
favors (hierarchical) p-extensions in the thickness direction,
yielding hierarchical plate and shell models. Finite element
computations show convergence properties and the efficiency of this important problem of boundary layer analysis
of plates and shells.
Chapter 9 to Chapter 11 of this Volume treat Generalized finite element methods. Chapter 9, this Volume,
by Ferdinando Auricchio, Franca Brezzi, and Carlo Lovadina, gives a systematic survey and some new results on
the stability of saddle-point problems in finite dimensions
for some classical mechanical problems, like thermal diffusion, the Stokes equations, and the Lamé equations.
Mixed methods yield the mathematical basis for problems
with locking and other numerical instability phenomena,
like nearly incompressible elastic materials, the Reissner
and Mindlin plate equations, and the Helmholtz equation.
From an engineering point of view, reduced integration
schemes and stabilization techniques get a sound foundation by the problem-dependent inf–sup condition. Chapter 10, this Volume, by Antonio Huerta, Ted Belytschko,
Sonia Fernández-Méndez, and Timon Rabczuk, provides
an advanced and systematic representation of different versions and alternatives of the so-called meshfree and particle
methods, known as moving least squares, partition of unity
FEM, corrected gradient methods, particle-in-cell methods


Fundamentals: Introduction and Survey

and so on. The method was originally invented for moving
singularities, and discontinuities like crack propagation in
solids, in order to avoid frequent complicated and costly
remeshings. These methods are based on Ritz–Galerkinand Petrov–Galerkin-type weighted residua or collocation
concepts and generalizations, such as Lagrange multiplier
and penalty methods. Radial basis functions are a good tool
(without having a compact support) as well as hierarchical
enrichments of particles.
The error-controlled approximation of the essential
boundary conditions of a boundary-value problem and, of
course, the related a priori and a posteriori error analysis as well as the relatively large condition number of the
algebraic systems combined with big computational effort
are crucial points.
In generally speaking, meshfree methods are now superior or at least equivalent to classical FEMs for some of the
addressed specific types of problems.
The last of the three chapters dealing with generalized FEMs, is Chapter 11, this Volume, by Nenced J.N.
Bicanic. Instead of constructing a convergent and stable
numerical method for the approximated solution of, for
example, a boundary-value problem for a continuous differential operator, a direct computational simulation of an a
priori discrete system with embedded discontinuous deformations, cracks, fragmentations, and so on is treated here.
This also includes assemblies of particles of different shapes
with their various contact problems, compactions, and other
scenarios of real processes. This is a rather new area of
computational mechanics, so far mostly treated on an engineering level, that is, without mathematical analysis. The
question arises how ‘convergence’ and ‘numerical stability’ can be defined and analyzed herein. But there is no
doubt that this type of direct computational simulation
of technological problems will play an important role in
the future.
The two chapters that follow are devoted to Boundary
element methods and their coupling with finite element
methods. Chapter 12, this Volume, by George C. Hsiao
and Wolfgang L. Wendland, represents variationally based
Galerkin-BEMs for elliptic boundary-value problems of
second order in a mathematically rigorous way, classified
by the Sobolev index. Various boundary integral equations
can be derived, introducing fundamental solutions, Green’s
representation formula, Cauchy data, and four boundary
integral operators. Out of this reservoir, several numerical
methods and algorithms for boundary elements are presented and discussed. The main features such as stability,
consistency, and convergence as well as adequate solvers,
condition numbers, and efficiency aspects are well treated.
Of course, error analysis and adaptivity play an important role. BEM has advantages over FEM in the case of

complicated boundaries, for example, mechanical problems
with edge notches and regular inner domains, and with
respect to dimensional reduction by one. Efficient recursive integration formulas and solvers for the fully populated
system matrices are available.
Chapter 13, this Volume, by Ernst Stephan, treats the
obvious variant of combining the different strengths of
both methods by symmetric couplings, which, of course,
need considerable algorithmic efforts and adequate solvers
with problem-dependent preconditioners. Special features
are Signorini-type contact problems using both primal and
dual-mixed finite element approximations. Recent features
are adaptive hp-methods. There seems to be a lack of
available software for 2D- – and even more for 3D- –
Chapter 14, this Volume, by J. Donea et al., is to be
seen separate from the previous presentations of different variational discretization methods, as it treats various
coupled processes – for example, fluid-solid interaction –
by suitable coordinates and metrics for each of the constituents – for example, Lagrangian coordinates for solids
and Eulerian coordinates for fluids – using the well-known
tangential push-forward and pull-back mappings between
the two descriptions via the deformation gradient. The
profit of computational efficiency and robustness can be
significant, for example, for the rolling contact of a tire
on a street. The ALE-concept, its analysis, the algorithms,
and important applications for linear and especially nonlinear static/dynamic problems in solid and fluid mechanics
are systematically presented and illustrated by adequate
examples. Also, smoothing and adaptive techniques for the
finite element meshes are discussed. It is remarkable how
quickly the ALE concept was implemented in commercial programs.
Chapter 15, this Volume, by Timothy Barth and Mario
Ohlberger, also stands separate from the scheme of finite
domain and boundary element approximations. Finite volume elements were invented in fluid mechanics and are also
applied now in other branches like biology, and in solid
mechanics, too. The advantage of finite volume approximations in comparison with the usual finite element discretizations in finite subdomains is the intrinsic fulfillment
of local conservation properties, like mass conservation or
entropy growth. Finite volume elements usually also yield
robust algebraic systems for unstructured meshes; they are
especially favorable for nonlinear hyperbolic conservation
systems in which the gradients of the solution functions
can blow up in time. Integral conservation laws and discrete volume methods are applied using various meshing
techniques of cell- and vertex-centered control volumes. A
priori and a posteriori error estimates are presented and

Fundamentals: Introduction and Survey 5
applied for adaptivity. Also, solvers in space and time
are discussed.
Chapter 17 to Chapter 18 of this Volume are devoted
to the Computer representation and visualization of
topology, geometry, meshes, and computed data.
Chapter 16, this Volume, by Franz-Erich Wolter, Niklas
Peinecke, and Martin Reuter, treats a subject growing in
importance in computational mechanics as technical objects
and their physical substructures become more and more
complicated. The presented methods and realizations can
be classified as computational methods for topology and
geometry with volume- and boundary-based (direct and
indirect) representations of objects and also with a rather
new type of modeling, using medial axes and surfaces for
describing objects that are mainly one or two dimensional
in their appearance. This medial modeling allows a natural
transition to finite element meshes. Of course, additional
attributes can be included in the geometry, like photometric
or toughness properties.
In Chapter 17, this Volume, by P.L. George et al., the
techniques of planar, surface, and volume meshing as well
as adaptive global and local remeshing for discretization
methods are outlined, aiming at automatic self-controlled
algorithms for various types of mesh generations and
their visualizations. Hard problems arise with large 3Dmeshes, requiring spatial decomposition, as well as related
automatic remeshing and local mesh adaptivity. Moving
boundaries and the meshing errors of specific elements
with respect to the given analytic or free-form geometry
are crucial problems.
The main methods for structured and unstructured mesh
generations are presented, where unstructured meshes are
constructed in a purely algebraic way or are based on appropriate PDE-solutions. Hierarchical spatial decompositions
are used for arbitrary shaped domains with unstructured
meshes, like the quadtree and octree method, advancing
front strategies and Delaunay type methods.
Chapter 18, this Volume, by William J. Schroeder and
Mark S. Shephard, also treats a crucial subject in modern
science and technology. Selected interactive visualization
of real and virtual objects and processes conveys intrinsic conceiving of the essentials that is hardly possible
with data files only. Visualization concerns geometry with
attribute data, meshes, and results that may need scalar,
vector, and tensor graphics with special features like producing streams of particles or physical quantities through
a total system or through a control volume. The presented
visualization algorithms give information about what is possible today.
Chapter 19 to Chapter 23 of this Volume treat the
crucial problem of stable robust and efficient solvers for
the various discrete algebraic systems introduced in the

first 14 chapters. Regarding the mostly high dimensions
of algebraic equation systems, which are needed today
in science and engineering and which can be solved
now in reasonable execution time with the teraflop generation of computers, a variety of sophisticated and
problem-oriented types of direct and iterative solvers
with adapted preconditioners were developed and are presented here.
Chapter 19, this Volume, by Henk A. van der Vorst,
presents direct elimination and iterative solution methods where the latter usually needs efficient preconditioning operators for efficient solutions. All important iterative solvers – based on Krylov projections – are treated,
like Conjugate Gradients, MINRES, OMR, Bi-CGSTAB,
and GMRES.
Special eigenvalue problems of a Hermitian matrix are
analyzed mainly with iterative QR-methods, emphasizing
tridiagonal and (upper) Hessenberg matrices. The Krylov
subspace approach is an efficient strategy that is applied
with four different versions. Several preconditioners are
presented and discussed.
In Chapter 20, this Volume, by Wolfgang Hackbusch,
fast iterative solvers for discretized linear and nonlinear
elliptic problems are treated, yielding (optimal) linear
complexity of the computational effort in regular cases,
which, of course, is of dominant importance to systems
with millions of unknowns for which multiprocessor and
massively parallel computers are also efficient.
Owing to the smoothing property of the Jacobi or
Gauss–Seidel iteration for elliptic PDEs on fine grids,
only data of coarse grid points are prolongated after
some smoothing steps, and this is repeated through several
grids, for example, in a W-cycle with four to five grid
levels, such that the solution takes place only with a
small equation system. The backward computation with
restriction operators and further smoothing steps finishes
an iteration cycle. It is efficient to produce a hierarchy of
discretizations that is easy for regular and nested grids. Such
a hierarchy may also be a side-product of adaptive mesh
Major issues of the chapter are the complete algorithmic
boxes as well as the analysis and error analysis of multigrid
methods for FEM and BEM.
Chapter 21, this Volume, by Wolfgang Hackbusch,
presents efficient solvers for fully populated matrices as
they arise in the boundary integral equation method.
The goal is the reduction of O(n2 ) arithmetic operations for standard matrix–vector multiplications to nearly
O(n) operations. The essential parts of the algorithm are
the far-field expansion and the panel cluster tree. A generalized variant is the construction of hierarchical matrices (Hmatrices) with different matrix–vector and matrix–matrix


Fundamentals: Introduction and Survey

operations only. The computation of inverse stiffness
matrices in FEM (e.g. for multiload cases) can be efficiently
computed by cluster techniques.
Chapter 22, this Volume, by V.G. Korneev and Ulrich
Langer, treats effective iterative solvers for large algebraic
systems by the alternating Schwarz method and advanced
substructuring techniques, emphasizing efficient problemdependent preconditioners. Nonoverlapping Schwarz methods are favored against overlapping methods, which need
more effort for computer implementation and the control of
the solution process. Of course, multiprocessor computers
are most suitable for parallel solution within the decomposed domains.
In Chapter 23, this Volume, by Werner C. Rheinboldt,
strategies for efficient solvers of highly nonlinear algebraic
systems, especially with physical instabilities like bifurcation and turning points, are investigated. These solvers are
based on Newton’s iterative method, its variants, and inexact Newton methods.
The problem of solution instabilities, which depend on
one scalar parameter, is analyzed by homotopy methods
and continuation methods. Also, the bifurcation behavior
of parameterized systems is investigated.
The last three chapters of Volume 1, Chapter 24, Chapter 25, and Chapter 26 of this Volume, are devoted to the
fundamentals of numerical methods for Time-dependent
Chapter 24, this Volume, by Kenneth Eriksson, Claes
Johnson, and Anders Logg, is based on duality techniques,
by solving associated linearized dual problems for the a
posteriori error analysis and adaptivity, using the residuals
of Galerkin approximations with shape functions, continuous in space and discontinuous in time, yielding the crucial
stability factors and discretization error estimates in adequate norms.
Parabolic initial boundary value problems are usually
stiff, and the main problem is the control of error accumulation in time. This is treated successfully with implicit
and explicit time-stepping methods for some classical
parabolic equations, like the instationary heat equation and
the reaction-diffusion problem.
In Chapter 2, Volume 2, by Martin Costabel, parabolic
and hyperbolic initial boundary value problems with transient solutions in time are considered, like heat con-

duction, diffusion, acoustic scattering, and elastic waves.
The following three approaches are critically compared:
space–time integral equations, Laplace-transform methods,
and time-stepping methods; many advanced mathematical
tools are necessary for the analysis, especially the error
analysis, which is treated here in a systematic way and
illustrated by examples.
Chapter 26, this Volume, by Leszek Demkowicz,
treats finite element approximations of the time-harmonic
Maxwell equations. With the stabilized variational formulations and Nedelec’s three fundamental elements, hpdiscretizations and hp-adaptivity are presented. Tetrahedral
elements of the first and second type, hexahedral elements
of the first type, prismatic elements as well as parametric
elements are treated.
The three Nedelec elements deal with the exact sequence
of gradient, curl, and divergence operators, yielding the
null-space, in addition to projection-based interpolations of
finite elements such that the element shape functions are
defined as a dual basis to the d.o.f.-functionals, aiming at
locality, global continuity, and optimality, and lastly the de
Rham commuting diagram property of the analytical and the
finite dimensional solution spaces, applying the gradient,
curl, and divergence operators.

We are convinced that the above sketched 26 chapters of
Volume 1 are a sound basis for today’s and tomorrow’s
computational mechanics, integrating mathematics, computer science, and physics, especially mechanics, as well as
challenging industrial applications that need high computer
power. Computer implementations will only be competitive
in engineering praxis if they are robust, stable, and efficient,
also concerning the requested logical clearness, as well as
the width and depth of the algorithms, as explained above.
An important benefit of this encyclopedia is seen in the
combination of Volume 1 with Volumes 2 and 3, which
are devoted to computational solid and fluid mechanics
such that users interested in theoretical issues and those
interested in practical issues can both get the information
they want, together with any secondary or background

Chapter 2
Finite Difference Methods
Owe Axelsson
University of Nijmegen, Nijmegen, The Netherlands

1 Introduction
2 Two-point Boundary Value Problems
3 Finite Difference Methods for Elliptic
4 Finite Difference Methods for Parabolic
5 Finite Difference Methods for Hyperbolic
6 Convection–Diffusion Problems
7 A Summary of Difference Schemes
Further Reading


part ni,j =1 aij (x)(∂ 2 u/∂xi ∂xj ). With no limitation, we can
assume that the coefficient matrix A(x) = [aij (x)]ni,j =1 is
As is well known, the differential equation is elliptic in
x if all eigenvalues of A have the same sign, hyperbolic
if one eigenvalue has an opposite sign to the others, and
parabolic if one eigenvalue is zero and the remaining have
the same sign.
We do not consider other cases here. Note that a differential equation can be of different types in different parts of the domain () of definition. For these
classes of differential equations, we need different types
of boundary conditions in order for the problem to be well

Definition 1. A boundary value problem is well posed if

Although, in general, more restricted in their applicability, finite difference methods provide a simple and readily
formulated approximation framework for various types of
partial differential equations. They can, hence, often compete with other methods such as finite element methods,
which are based on certain variational formulations of the
differential equations, and where the preparation work to
construct the corresponding systems of algebraic equations
is more involved.
When constructing the approximation scheme, it is
important to know the type of the given differential
equation. Partial differential equations of second order
are classified according to the type of their principal
Encyclopedia of Computational Mechanics, Edited by Erwin
Stein, René de Borst and Thomas J.R. Hughes. Volume 1: Fundamentals.  2004 John Wiley & Sons, Ltd. ISBN: 0-470-84699-2.

(i) Existence: There exists a solution that satisfies the
equation and each given boundary condition (assuming they are sufficiently smooth).
(ii) Uniqueness: It has at most one solution.
(iii) Stability: Its solution is a continuous function of the
given data, that is, small changes in the given data
entail small changes in the solution.
Although it is of interest in practice to consider also
certain ill-posed problems, these will not be dealt with here.
Solutions of partial differential equations of different
types have different properties. For instance, the solution to
elliptic and parabolic problems at any given point depends
on the solution at all other points in . However, for
hyperbolic problems, the domain of dependence is a subset
of . Although some problems may belong to the class of
elliptic problems, they may exhibit a dominating hyperbolic
nature in most parts of . To illustrate the aforesaid,
consider the following examples.


Finite Difference Methods
the corresponding homogeneous equation. By multiplying
the equation by v and integrating, we get

The problem
Lu ≡ −εu + v · ∇u = 0 in ,

u = g on ∂


where |v|  ε, ε > 0 has a dominating hyperbolic nature in
the interior domain, except near the outflow boundary part,
where v · n ≥ 0, and where the diffusion part dominates and
the solution has a thin layer. Here n is the outward-pointing
normal vector to .
For the Poisson problem
−u = ρ(x),



where we ignore the influence of boundary data, the solution is

|x − ξ|−1 ρ(ξ) dξ
u(x) =
called the gravitational or electric potential, where ρ is the
density of mass charge. Hence, the solution at any point in
 depends on all data. On the other hand, the first-order
hyperbolic problem
aux + buy = 0,

x > 0, y > 0


where a, b are positive constants, and with boundary conditions u(0, y) = g(y), y ≥ 0 and u(x, 0) = h(x), x ≥ 0 has
a solution
x≥ y
 h x − by
u(x, y) =
g y − x
x< y
Therefore, the solution is constant along the lines bx −
ay = const, called characteristic lines, so the solution at
any point depends just on a single-boundary data.
Depending on its type, there exist various methods to
prove uniqueness and stability of a boundary value problem. For elliptic and parabolic problems, one can use a
maximum principle, which shows pointwise stability. As
an example, consider −u = f in , u = g on ∂, where
f ≤ 0. Then, assuming first the contrary, it is readily proven
that the maximum of u is taken on the boundary ∂. This
also shows that any perturbation of boundary data by some
amount ε cannot result in any larger sized perturbation of
the solution in the interior of the domain. An alternative
way of proving the maximum principle is via a Green’s
function (fundamental solution) representation of the solution such as in (2).
For hyperbolic problems, one uses the energy integral
method. Let u1 and u2 be two solutions of ut = aux ,
a > 0, 0 < x < 1, t > 0, where u(0, t) = 0, which correspond to different initial data. Then, v = u1 − u2 satisfies



vt v dx +

or, by letting E(t) =
we find



avx v dx = 0

v 2 dx (called the energy integral),

1 d
E(t) + v 2 (1, t) = 0,
2 dt


that is, E  (t) ≤ 0. Hence, E(t) ≤ E(0), t > 0, where
E(0) = 0 (u1 (x, 0) − u2 (x, 0))2 dx, which shows stability
and uniqueness.
Clearly, the energy method is also applicable for more
general problems, where the analytical solution is not
known. For the discretized problems, one can use a discrete
maximum principle, enabling error estimates in maximum
norm. For parabolic and hyperbolic problems, we can also
use stability estimates based on eigenvalues of the corresponding difference matrix or, more generally, based on
energy type estimates. Finally, for constant coefficient problems, Fourier methods can be used.
Let the domain of definition of the continuous problem
be discretized by a rectangular mesh with a mesh size h. In
order to compute a numerical solution of a partial differential equation like (1), we replace the partial derivatives in
the differential equation by finite differences. In this way,
we obtain a discrete operator equation
Lh uh = fh
We will be concerned with the well-posedness of the
discrete equation, using a discrete maximum principle, for
instance, as well as with estimates of the discretization
error eh = u − uh , where u is the solution of the original
continuous partial differential equation. This estimation will
be based on the truncation error τh = Lh u − fh .
The following finite difference approximations will be
used throughout this text:
u (x)  Dx+ u(x) :=

u (x)  Dx− u(x) :=

u (x)  Dx0 u(x) :=

[u(x + h) − u(x)],
the forward difference


[u(x) − u(x − h)],
the backward difference


[u(x + h) − u(x − h)],
the first-order central difference (6)

Finite Difference Methods 9
Note that Dx0 u = (1/2)(Dx+ + Dx− )u. More generally,
one can use the so-called “θ-method”

u (x) 


+ (1 −

θ)Dx− ]u(x)


where θ is a method parameter. Thus,
u (x) = [θDx+ + (1 − θ)Dx− ]u(x) +

(1 − 2θ)u (x)

[(1 − θ)u(3) (x − η2 h) + θu(3) (x + η1 h)]

where 0 < ηi < 1, i = 1, 2, so

if |1 − 2θ| = O(h)

O(h2 )

[θu(x + h) + (1 − 2θ)u(x) − (1 − θ)u(x − h)]
+ (1 − 2θ)u (x) − u(3) (x + η3 h), −1 < η3 < 1

u (x) =

Note that for θ = 1/2, we get (6), for θ = 1, we get (4),
and for θ = 0, we get (5). Hence, the θ-method generalizes
the previous methods.
An approximation for the second derivative is obtained
u (x)  Dx+ Dx− u(x) ≡ Dxx

the central difference of second order (8)
u(x) = Dx+ Dx− u(x) = h−2 [u(x + h) −
We have then Dxx
2u(x) + u(x − h)], thus

− u(4)
(x + η4 h),
12 x


+ O(h ),

u (x) =

u (x) =


+ u(x − h, y) + u(x, y + h) + u(x, y − h) − 4u(x, y)]
(5) is called the 5-point difference operator.
Various difference methods adjusted to the type of the
problem, with discretization error estimates based on truncation errors, are presented.
When deriving truncation error estimates for difference
approximations, one uses Taylor expansion (assuming sufficient regularity of u)
(k − 1)!

× hk−1 u(k−1) (xi ) + R(xi , h, k)

If 0 ≤ θ ≤ 1, we get

(5) u := [Dx+ Dx− + Dy+ Dy− ]u(x, y) = h−2 [u(x + h, y)

u(xi + h) = u(xi ) + hu (xi ) + · · · +

u (x) = [θDx+ + (1 − θ)Dx− ]u(x)

For hx = hy = h, we have


where the remainder term R(xi , h, k) can be written as
R(xi , h, k) = (1/k!)hk u(k) (ξi ), ξ ∈ (xi , xi + h) or in the
x +h
alternative form R(xi , h, k) = xi i [1/(k − 1)!](xi + h −
s)k−1 u(k) (s) ds.

The most common among problems of applied mathematics type that appear in physics, engineering, and so on are
boundary value problems for partial differential equations.
As an introduction to difference methods for such problems,
we consider here the corresponding problem in one dimension, the two-point linear differential equation problem:
Find u ∈ C 2 [a, b] such that
Lu ≡ −(k(x)u ) + p(x)u = f (x),

a < x < b (11)

with boundary conditions
−1 < η4 < 1

r0 (u) ≡ γ0 u(a) − δ0 k(a)u (a) = α



r1 (u) ≡ γ1 u(b) + δ1 k(b)u (b) = β

Similar expressions hold for Dy+ , Dy− , Dy0 , and so on. In
particular, if u(4)
x , uy ∈ C(),
Dx+ Dx− u(x, y) + Dy+ Dy− u(x, y)
= h−2
x [u(x + hx , y) + u(x − hx , y) − 2u(x, y)]
+ h−2
y [u(x, y + hy ) + u(x, y − hy ) − 2u(x, y)]
= uxx (x, y) + uyy (x, y) + O(h2x ) + O(h2y ), hx , hy → 0

Here, u = du/dx, k(x) ≥ k0 > 0, a ≤ x ≤ b and k ∈
C 1 (a, b), and p, f ∈ C(a, b) are given real-valued functions and γ0 , δ0 , γ1 , δ1 , α, β are given real numbers. The
operator L is self-adjoint, that is, a Lu v dx = a Lv u dx
The solution u will then be a twice continuously differentiable function.
Such problems arise, for instance, if we let u be the
displacement of a (thin) elastic string subjected to forces
with distribution defined by f . In the simplest model, k
is constant, p(x) ≡ 0, and the string is fixed at (a, α),


Finite Difference Methods




Find u ∈ C 2 [a, b] such that
Lu ≡ −u + p(x)u = f (x),

f (x )

u(a) = α, u(b) = β


u (x )

Figure 1. An elastic string subjected to forces.

(b, β), see Figure 1. This equation also arises if we let u be
the temperature in a heat-conducting rod that is perfectly
insulated except at the endpoints, where the temperatures
α, β are given. The material coefficients k and p are often
assumed to be constants.
Existence and uniqueness
For boundary value problems, there are intrinsic questions
regarding existence and uniqueness of solutions. There
holds the following theorem.
Theorem 1 (Fredholm’s alternative) Equation (11) has
exactly one solution if and only if the homogeneous problem
(with f ≡ 0 and homogeneous boundary conditions α =
β = 0) only has the trivial solution.
We now give conditions under which the boundary value
problem (11) has exactly one solution.
Theorem 2. If k(x) > 0, p(x) ≥ 0, γi2 + δ2i > 0, δi ≥ 0,
γi ≥ 0, i = 0, 1 and if any of the following conditions holds,
(i) γ0 = 0, (ii) γ1 = 0, (iii) p(x) > 0 for some x ∈ (a, b),
then the boundary value problem (11) has exactly one

That the homogeneous problem −(ku ) + p(x)u = 0,
a < x < b, r0 (u) = r1 (u) = 0, where the boundary operators r0 , r1 are defined in (12), has only the trivial
solution, can be shown by partial integration of
Lu = 0.


a < x < b (13)

Difference approximations

The boundary value problem (11) cannot, in general,
be solved analytically. We shall now present a simple but efficient numerical method. In order to simplify
the presentation, we shall mostly consider the following

where p(x) ≥ 0, a < x < b, and max(α, β) ≥ 0.
Let x0 = a, xi = xi−1 + h, i = 1, 2, . . . , N , xN+1 = b,
where h = (b − a)/(N + 1) be a uniform partitioning π =
πh of the interval [a, b]. In order to find an approximate
solution at the points xi , we approximate u by finite
−u (xi )  h−2 [−u(xi−1 ) + 2u(xi ) − u(xi+1 )]
Letting uh be the resulting approximation of u at xi ,
i = 1, 2, . . . , N , we get
Lh uh (xi ) ≡ h−2 [−uh (xi−1 ) + 2uh (xi ) − uh (xi+1 )]
+ p(xi )uh (xi ) = f (xi ),
i = 1, 2, . . . , N, uh (a) = α, uh (b) = β


This is a linear system of N equations in the N unknowns
uh (xi ), i = 1, 2, . . . , N . That the corresponding matrix is
nonsingular follows from a discrete maximum principle.
Lemma 1 (Discrete maximum principle) Let Lh be
defined by (14) and vh be any function with max{vh (x0 ),
vh (xN+1 )} ≥ 0 defined on xi for which Lh vh (xi ) ≤ 0, i =
1, 2, . . . , N . Then,
max vh (xi ) = max vh (xi )



that is, the largest value of vh is taken at one of the
Next, we show the corresponding error estimates for τh
and eh . We now define the truncation error of the difference
approximation as the error (or defect) we get when we
substitute the difference approximation uh in (14) for the
exact solution u of (13).
Definition 2. The function τh,i = (Lh u)i − fi , i = 1, . . . ,
N is referred to as the truncation error τh,i of the difference
approximation (14) at xi .
Definition 3. The function eh (xi ) = u(xi ) − uh (xi ), i =
0, 1, . . . , N + 1 is referred to as the discretization error of
the difference approximation (14) at point xi .

Finite Difference Methods 11
Note that the truncation error and the discretization error
are related by
Lh eh = τh


It turns out that τh is easily estimated via Taylor expansions.
We see from (15) that in order to estimate eh , we need a
bound of the inverse of the discrete operator Lh (if it exists).
Such a bound is provided by the next barrier function
lemma, which can be proven by the discrete maximum
In order to estimate eh , it is convenient to use the
following notations:
|vh |πh = max |vh (xi )|,


Lemma 2 (Barrier function lemma) Assume that there
exists a function wh (called barrier function), defined on
xi , i = 0, 1, . . . , N + 1, which satisfies Lh wh > 0, wh ≥ 0
and wh (x0 ) = 0. Then, any function vh defined on xi , i =
0, 1, . . . , N + 1, for which vh (x0 ) = 0 satisfies
maxπh wh
minπh Lh wh

|Lh vh |πh


1 2 (4)
h u (ξi ),

For a uniform mesh and sufficiently smooth solutions,
the order of accuracy of the approximate solution can be
easily improved using the classical trick of Richardson
extrapolation and nested meshes involving approximations
at same meshpoints for two different values of h.
Theorem 4. Let u ∈ C 6 (a, b) and let uh be the solution of
(14). Then,

h → 0,

i = 1, 2, . . . , N

where ϕ is a function that is independent of h.

−ϕ + pϕ = −

In order to prove that the discretization error is also
O(h2 ), we use Lemma 2. As a barrier function, we use
b + a − 2xi



Note that wh (a) = wh (b) = 0.
Theorem 3. The discretization error of the difference
approximation (14) satisfies
|u − uh |πh ∪∂πh ≤

2.2 Richardson extrapolation

Proof. Letting ϕ be the solution of the auxiliary problem
xi−1 < ξi < xi+1

Proof. Since fi = f (xi ) = p(xi )u(xi ) − u (xi ) and (Lh u)i
= h−2 [−u(xi−1 ) + 2u(xi ) − u(xi+1 )] + p(xi )u(xi ), (a) follows. (b) follows now by Taylor expansion around xi .
Hence, if u ∈ C (4) (a, b), the truncation error is O(h2 ),

h → 0.

wh (xi ) = 1 −

u(xi ) − uh (xi ) = h2 ϕ(xi ) + O(h4 ),

Lemma 3. (a) τh,i = h−2 [−u(xi−1 ) + 2u(xi ) − u(xi−1 )]
+ u (xi ).
(b) If u ∈ C (4) (a, b), then
τh,i = −

(b − a)2
|Lh (u − uh )|πh
(b − a)2
(b − a)2
|Lh u − f |πh =
|τh |πh
(b − a)2 h2
max |u(4) (ξ)|

|u − uh |πh ∪∂πh ≤ |u − uh |∂πh +

|vh |∂πh = max |vh (xi )|

|vh |πh ∪∂πh ≤ |vh |∂πh +

Proof. Note that wh , defined by (18), satisfies 1 ≥ wh
≥ 0 on πh and, as is easily seen, Lh wh = 8/(b − a)2 +
pi wh (xi ) ≥ 8/(b − a)2 > 0. By Lemmas 2 and 3(b), we
now get

[(b − a)h]2 max |u(4) (ξ)|


1 (4)
u ,

a < x < b, ϕ(a) = ϕ(b) = 0

and assuming that ϕ ∈ C 4 (a, b), which holds if u ∈ C 6
(a, b), it follows that Lh (e − h2 ϕ)i = O(h4 ), h → 0, where
e = u − uh . Applying the Barrier Lemma completes the

This argument can be repeated so that if u is smooth
enough, one can prove an h2 expansion of the error u(xi ) −
uh (xi ). Hence, repeated Richardson extrapolation may also
be applied. Such error expansions do not, however, always
exist. For instance, if one of the boundary points is not a
meshpoint in the uniform mesh, the distance δh, 0 < δ < 1
of this boundary point to the nearest meshpoint is not a
smooth function of h.
The extrapolation procedure has advantages over conventional higher-order methods. Thus, the basis difference
formula can be very simple, which makes repetition with a
new mesh size easy and the method automatically provides
error estimates.


Finite Difference Methods

For extensions of such asymptotic error expansions and
Richardson extrapolation for linear finite elements; see
Blum, Lin and Rannacher (1986).

results on difference methods for elliptic problems can be
found in Section 6.

3.1 Difference approximations

Computing high-order approximations of the
derivatives of the solution

We now show that not only the solution but also its
derivatives can be computed to high accuracy, using the
already computed approximate values of the solution.
Theorem 5. Let u ∈ C 6 (a, b) be the solution of −u +
p(x)u = f , a < x < b, u(a) = α, u(b) = β and let uh be
the solution of the discrete problem, discretized using central
differences on a uniform mesh. Then,
u (xi ) =

Consider first the Poisson problem u = f with given
Dirichlet boundary conditions on a unit square domain,
aligned with the coordinate system, and a uniform rectangular mesh h with mesh sizes h1 = ρ1 h and h2 = ρ2 h
in the x- and y-direction respectively. Here, ρ1 and ρ2
are given positive numbers, chosen such that 1/ρi h is an
integer, and h is a small positive parameter, which we
intend to make arbitrarily small to get a sufficiently accurate
numerical solution. In each interior meshpoint in h (see
Figure 2a), we use the five-point difference approximation,

uh (xi+1 ) − uh (xi−1 )
+ O(h2 )

x, y + h2

Proof. Theorem 4 shows that uh (xi ) = u(xi ) − h2 ϕ(xi ) +
O(h4 ), where ϕ does not depend on h and ϕ ∈ C 2 (a, b).
Hence, using Taylor series expansions,

x + h1 , y

x, y − h2

In a similar way, assuming a correspondingly higher
order of regularity of the solution, even higher-order derivatives can be computed with error O(h2 ). This result is not
obvious, because to compute an approximation of u , we
make use of approximations of u, divided by h or a higher
power of h. That we do not loose one or more power(s) of
h in the order of approximation is due to the existence of
an h-expansion of the errors.


x, y


x − h1, y

u(xi+1 ) − u(xi−1 )
uh (xi+1 ) − uh (xi−1 )
+ O(h3 ) = u (xi )
− h2
+ u(3) (η1 ) − h2 ϕ (xi ) + O(h3 )
where η1 ∈ (xi−1 , xi+1 ).







We present in this section various difference methods for
the numerical solution of partial differential equations of
elliptic type. Discretization errors are derived for operators
of positive type. The derivations are done for problems
in two space dimensions but most results can be readily
extended to problems in three space dimensions. Further







Figure 2. Local difference meshes (a) rectangular domain, (b)
curved boundary.

Finite Difference Methods 13
+ −
+ −
h u := (Dh1 Dh1 + Dh2 Dh2 )u(x, y) and on the boundary
points, we use the given boundary values.
Let hi = 1/(Ni + 1), i = 1, 2 and, hence, let ρ1 /ρ2 =
(N1 + 1)/(N2 + 1), where Ni + 1, i = 1, 2 are the number
of meshpoints in the two coordinate directions. Let uh
denote the corresponding approximate solution. We then get
a system of N1 N2 linear algebraic equations of the form

uh (x − h1 , y) − 2uh (x, y) + uh (x + h1 , y)
Lh uh := h−2

+ h−2
2 [uh (x, y − h2 ) − 2uh (x, y) + uh (x, y + h2 )]
= f (x, y),

x, y ∈ h


If we order the meshpoints in lexicographic (horizontal)
order, the corresponding matrix takes the form
Ah = block tridiag(Ii , A(i) , Ii )
where A(i) = h−2
1 tridiag (1, −2, 1) − 2h2 I and Ii =
and the identity matrix I
2 I , i = 1, 2, . . . , N2 . Here, A
have order N1 × N1 and there are N2 block-rows in Ah .
Systems with matrix Ah can be solved efficiently using various methods such as fast direct methods or multigrid and
algebraic multilevel iteration methods; see, for example,
Axelsson and Barker (1984) and Hackbusch (1985).
Assuming that u ∈ C 4 (), the local truncation error
Lh u − f of the difference approximation is readily found
to be

1 2 (4)
h1 ux (x + ξ1 , y) + h22 u(4)
Lh u − f =
y (x, y + ξ2 )
= O(h ),

−hi < ξi < hi , i = 1, 2

Curved boundaries
For a more general domain , for instance with a curved
boundary such as illustrated in Figure 2b, we must modify
the approximation scheme at interior points in h next to
the boundary. There are two such efficient schemes. The
first uses a generalized five-point difference approximation
with, in general, nonequidistant meshlengths (see Shortley
and Weller, 1938):

[u (E) − uh (P )] −
[u (P ) − uh (W )]
Lh uh :=
hE h
hW h

[u (N ) − uh (P )]
hW + hE
hN h

= f (P ), P ∈ h
− [uh (P ) − uh (S)]
hS + hN
where hE , hW , hN , hS denote the distances from P to the
surrounding points in the East, West, North, and South

directions, see Figure 2. Unless hE = hW and hN = hS , the
local truncation error is O(h) here. The coefficient matrix is
in general not symmetric. The second method uses a linear
combination of weighted linear interpolations in the x- and
Lh uh :=

{h u (W ) + hW uh (E) − (hE + hW )uh (P )}
hE E h

{h u (S) + hS uh (N ) − (hN + hS )uh (P )}
hN N h

h + hW
h + hN
hW E
f (P ), P ∈ h
+ hS S
Here, the coefficient matrix is symmetric.

Remark 1. To provide an alternative to the Shortley–Weller approximation for treating curved boundaries,
much work has been devoted to numerical grid generation
during the past 20 to 30 years. For instance, curvilinear
boundary-fitted finite difference methods became popular and applied extensively in numerical fluid dynamics
problems (see Thompson, Warsi and Mastin, 1985). More
recently, much effort has been devoted to variational grid
generation methods, which can provide more robust methods applicable also for very complicated geometries; see,
for example, Garanzha (2004) and the references there in.

3.2 Higher-order schemes
3.2.1 The nine-point difference scheme
Let u = f , and consider first the cross-directed five-point
difference scheme on a local equidistant square submesh,
= 12 h−2 [uh (x − h, y − h) + uh (x + h, y − h)
+ uh (x − h, y + h) + uh (x + h, y + h)
− 4uh (x, y)],

x, y ∈ h


It is readily seen that, for a sufficiently smooth function u,
h = u +


2 2  (4)
+ h4 u(6)
h ux + u(4)
x + uy

+ O(h6 )

= u + h2 u(4)
+ h4
x + 6uxxyy + uy

× ux + 15uxxxxyy + 15uxxyyyy + u(6)
+ O(h6 )

Let (9)
h be the nine-point difference scheme defined by
h =

2 (5) 1 (5,×)
 + h
3 h


Finite Difference Methods

The coefficients in this stencil equal 1/6 for the corner
vertex points in the square with edges 2h, equal 2/3 for the
midedge points, and equal −10/3 for the center point.
A computation shows that for a uniform rectangular
h uh

1 4 2
= f + h2 f +
h  f + 2fxxyy + O(h6 )

where 2 = (f ). Using a modified right-hand side
in the difference formula, it follows that the difference

h2 (5)

12 h

(x, y) ∈ h


has truncation error O(h4 ).
Further, it follows from (24) that for a sufficiently smooth
2 2
function f , f = (9)
h f − (1/12)h  f + O(h ). A
f −  f]+
computation shows that h fxxyy = 2[
O(h ) and, therefore, the nine-point stencil with the next
modified right-hand side
h uh = f +

1 2 (9)
1 4 (5)
h h f −
h h f

1 4
h fxxyy ,

(x, y) ∈ h

has a truncation error O(h6 ).
The implementation of this scheme is simplified if f
is given analytically so that f and so on can be computed explicitly. If f ≡ 0, then (9)
h uh ≡ 0 has an order
of approximation O(h6 ). Hence, this scheme provides very
accurate approximation, for instance, for far-field equations,
where frequently u = 0.
The above is an example of a compact difference scheme
(for further references on such schemes, see Houstis and
Rice, 1982).

3.2.2 Difference methods for anisotropic problems
and problems with a mixed derivative
Consider first the anisotropic differential equation auxx +
buyy = f (x, y), (x, y) ∈ , where u = g(x, y), (x, y) ∈
∂ and f and g are given, sufficiently smooth functions.
Let a > 0 and b > 0.
Here, the nine-point difference approximation has a stencil, as shown in (25) with c = 0. If we modify the righthand side to be f + 1/12h2 f , it can be seen that in this
case the local truncation error becomes O(h4 ).
Consider next the differential equation with a mixed
auxx + 2cuxy + buyy = f (x, y),

(x, y) ∈ 

with given boundary conditions. We assume that a > 0,
b > 0, and c2 < ab, which are the conditions for ellipticity
of the operator. For the mixed derivative, we use the central
difference approximation uxy ≈ Dx0 Dy0 u, that is,
uxy ≈

[u (x − h, y − h) − uh (x + h, y − h)
4h2 h
− uh (x − h, y + h) + uh (x + h, y + h)]

Combined with the nine-point difference stencil, it becomes

1 −2 
h 

− 3c

5a − b
+ 3c

5b − a
−10(a + b)
5b − a

+ 3c
5a − b  (25)
− 3c

3.2.3 Difference schemes for other regular
Finite differences can be extended to nonrectangular
For a regular (isosceles) triangular mesh, one can form
the obvious seven-point difference stencil. For a hexagonal
(‘honeycomb’) mesh, one finds a four-point stencil




The symmetrically located nodepoints in the seven-point
scheme allows one to readily approximate second-order
cross derivatives. Similarly, in 3D problems, a cuboctahedral stencil involves 13 nodepoints. If a Cartesian grid
is used, approximations of the second-order cross derivatives require at least nine points in 2D and 19 points in
3D, that is, two and six more than for the hexagonal and
cuboctahedral stencils.
The biharmonic operator 2 u = f , (x, y) ∈  with
boundary conditions such as u = g(x, y), ∂u/∂n = q(x, y)
on ∂ give rise to a 12-point stencil for a regular equidistant
2 −8
−4 
h  1 −8 20 −8 1 
2 −8

Finite Difference Methods 15
which has truncation error O(h2 ). The biharmonic problem
can, however, more efficiently be solved as a coupled


ξ − u = 0
ξ = f


using a variational formulation hereof; see, for example,
Axelsson and Gustafsson (1979) and Ciarlet and Raviart



3.3 Approximating boundary conditions
So far, we have only considered Dirichlet type boundary conditions u = g on ∂. For a Neumann (∂u/∂n :=
∇u · n = g) or the more general Robin (∂u/∂n + σu = g)
(cf. Gustafson and Abe, 1995) type boundary conditions,
one must use special approximation methods. The regular
difference mesh is then extended around the boundary line.
If the normal to the boundary goes through meshpoints,
we can use a central difference approximation for ∂u/∂n,
using an interior meshpoint and its symmetrically reflected
point in the extended domain (see Figure 3(a)). In other
cases, we can still use central difference approximations
for ∂u/∂n (at U, R in Figure 3(b)) but we must interpolate the function value in the symmetrically located points
in the interior. This can be done using linear interpolation
from the surrounding points (P, N, NW, W) in Figure 3(b)
to find the value at T , or using biquadratic interpolation,
involving some additional points. The local truncation error
becomes O(h) or O(h2 ), respectively. It can be seen that
one can always get a positive-type scheme in the first case,
but not in the second case.
For discretization errors for problems with curved boundaries, see Section 6.

3.4 Monotonicity and discretization error
Monotone operators provide a basis for pointwise bounds
of the solution and the discretization errors corresponding
to various difference approximations.
The general form of a linear difference operator Lh
depending on some discretization parameter h is
Lh uh (Pi ) ≡


aij (Pi )uh (Pi ) = f(Pi ),

i = 1, 2, . . . , n

j =1

Here, the function f includes the given source function and
the given boundary data.
The operator is said to be monotone if Lh v ≥ 0 implies
v ≥ 0, where v is a function defined on h . Note that

















Figure 3. Neumann boundary conditions.

a monotone operator is nonsingular because if Lh is
monotone and Lh v ≤ 0, then Lh (−v) ≥ 0, so −v ≥ 0, that
is, v ≤ 0. Hence, if Lh v = 0, then both v ≥ 0 and v ≤ 0
hold, so v = 0.
It is also readily seen that a nonsingular matrix A is
monotone if and only if A−1 ≥ 0, where the inequality
holds componentwise. Further, for any monotone operator
A, there exists a positive vector v > 0 such that Av > 0.
Take, for example, v = A−1 e, where e = [1, 1, . . . , 1]T . If
some component vi = 0, it would follow that all entries of
the ith row of A−1 were zero, which is clearly not possible.
Consider now a monotone difference operator Lh and
a normalized function w with maxi w(xi ) = 1, called a
barrier function, such that Lh w(x) ≥ α for all x ∈ h and
for some positive constant α. Then,





Finite Difference Methods

Lh ∞ = supv=0 Lh v ∞ / v ∞ and
v ∞=
max |v(xi )|, xi ∈ h is the supremum norm.
As is readily seen, this yields a discretization error


≡ u − uh


L u − fh
α h


where u is the exact solution to a given differential
equation Lu = f , Lh is a discrete approximation to L,
uh is the solution to the corresponding discrete equation
Lh uh = fh , and Lh u − fh = Lh u − Lh uh is the truncation
The barrier function thus leads to an easily computable
discretization error estimate if the discrete operator is monotone. In addition, as pointed out by Collatz (1986), the
monotonicity property is of great practical value because
it can readily be used to get lower and upper bounds for
the solution.
Monotone operators arise naturally for positive-type difference schemes. The difference operator Lh is of positive
type if the coefficients satisfy the sign pattern aii > 0 and
aij ≤ 0, j = i.

If, in addition, aii ≥ jn=1 |aij |, i = 1, 2, . . . , n, j = i
with strong inequality for at least one index i, then the
operator is monotone. This is equivalent with that the matrix
A = [aij ]ni,j =1 is an M-matrix.
Stieltjes proved that a positive definite operator of positive type is monotone.
However, even if the operator is not of positive type, it
might nevertheless be monotone. For instance, a familiar
result states that if A = M − N is a weak regular splitting
(i.e. M is monotone and M −1 N ≥ 0), then A is monotone
if and only if the splitting is convergent, namely, there
holds ρ(M −1 N ) ≤ 1 with ρ being the spectral radius of
M −1 N .
Hence, monotonicity of a given matrix A, respectively,
a linear discrete operator Lh , can be proven by constructing convergent weak regular splittings A = M − N . As is
shown below, this result can be extended to matrix splittings
of a more general form.

3.4.1 Bounding inverses of monotone matrices
To bound the supremum norm of the inverse of a nonsingular monotone matrix A, we consider the following Barrier
Lemma both in its general and sharp forms.
Lemma 4 (The Barrier Lemma) Let A be a monotone
matrix of order n and let v be a normalized vector, v ∞ =
1, such that mini (Av)i ≥ α for some positive scalar α.



≤ 1/α;

A−1 ∞ = 1/ max{mini (Av)i , v ∈ VA }, where VA =
{v ∈ Rn , v ∞ = 1, Av > 0};
(c) A−1 ∞ = x ∞ , where x is the solution of Ax = e.


Proof. For a proof, see Axelsson and Kolotilina (1990).

For later use, note that for a strictly diagonally dominant
matrix A = [aij ], it holds the inequality





min |aii | −
i 
We now extend
positive vector v
This result can be
conditions hold at

|aij |

j =i


the barrier lemma to the case where the
satisfies the weaker condition Av ≥ 0.
particularly useful if Dirichlet boundary
some part of the domain .

Lemma 5. Let

 A be monotone of the form A =
A11 −A12
, where A11 is monotone, A12 and A21
are nonnegative and A22 has no zero rows. Further, let
v = [v1 ; v2 ]T be a positive vector such that v ∞ = 1 and
A[v1 ; v2 ] = [p; q]T , p > 0, q ≥ 0. Then there holds
12 ∞  
21 ∞ 
A−1 ∞ ≤ 1 +
min(A11 v1 )i
min(A11 v1 )i
× max
 min(A11 v1 )i min(q + A21 A11 p)i 


A proof can be found in Axelsson and Kolotilina (1990).

3.4.2 Proving matrix monotonicity
Here, we summarize some classical results on weak regular
splittings, convergent splittings, and Schur complements of
matrices partitioned into a two-by-two block form, which
can be used to ascertain that a given matrix is monotone.
Let M, N ∈ Rn×n ; then, A = M − N is called a weak
regular splitting if M is monotone and M −1 N is
nonnegative. The splitting is convergent if ρ(M −1 N ) < 1.
Theorem 6. A weak regular splitting P A = M − N ,
where P is nonsingular and nonnegative, is convergent if
and only if A is monotone.
Proof. See Axelsson and Kolotilina (1990).

In practical applications, it can be more convenient
to use, instead of Theorem 6, the following sufficient

Finite Difference Methods 17
Theorem 7. Let PAQ = M − N be a weak regular splitting with P and Q nonsingular and nonnegative. Then, A is
monotone if there exists a positive vector v such that either
M −1 PAQv > 0 or vT M −1 PAQ > 0.
Proof. Since by assumption M − N is a weak regular
splitting, it follows by Theorem 6 that PAQ is nonsingular
and monotone if ρ(M −1 N ) < 1. But M −1 PAQv = (I −
M −1 N )v or, since M −1 N is nonsingular, 0 ≤ M −1 N v =
(I − M −1 PAQ)v < v if M −1 PAQv > 0. Hence, with
D = diag(v1 , v2 , . . . , vn ), that is, De = v, 0 ≤ D −1 M −1
N De < e or D −1 M −1 N D ∞ < 1, so ρ(M −1 N ) < 1.
In a similar way, if vT M −1 PAQ > 0, it follows that

DM −1 N D −1 1 < 1, where · 1 is the 1 -norm.
Remark 2. Theorem 7 can be particularly useful when
we have scaled A with diagonal matrices P and Q.

where A11 and A22 are square submatrices and A11 is
(a) If A11 is monotone and A−1
11 A12 and A21 A11 are nonnegative, then A is monotone if S = A22 − A21 A−1
A12 is monotone.
(b) Conversely, assume that A is monotone. Then S is

Proof. To prove the existence of the inverse in part (a),
use the block matrix factorization of A, which shows that
A is invertible if and only if A11 and S are invertible.
The monotonicity statements in both parts follow from
the explicit form of the inverse of A,



11 + A11 A12 S A21 A11

11 A12 S

S −1 A21 A−1

S −1


From Theorem 7, one can deduce the following important
monotonicity comparison condition.
Corollary 1. Let B1 ≤ A ≤ B2 , where B1 and B2 are
monotone matrices. Then, A is monotone and B2−1 ≤ A−1 ≤
B1−1 .
Proof. See Axelsson and Kolotilina (1990).

Theorem 8. Let A = [Aij ] be an m × m block matrix
satisfying the following properties
(i) Aii are nonsingular and A−1
ii > 0, i = 1, 2, . . . , m.
(ii) For i = j , there exist matrices Pij ≥ 0 such that
Aij ≤ −Pij Ajj .
(iii) There exists a positive vector v such that either the
block components (AT v)i are nonzero and nonnegative for i = 1, 2, . . . , m or Au > 0, where ui = A−1
ii vi ,
i = 1, 2, . . . , m.

Example 1.

Consider the mixed derivative elliptic

−auxx − 2cuxy − buyy = f in  = [0, 1]2


u = 0 on ∂, with variable coefficients a(x, y), b(x, y) >
0, and c(x, y) ≥ 0. After elimination of boundary conditions, the standard nine-point second-order accurate finite
difference approximation of this problem on a uniform
mesh, yields the block tridiagonal n2 × n2 matrix
A = 2 block tridiag T


, −bi , i ,

T (−ai , 2(ai + bi ), −ai ), T

− i , −bi , − i


Then A is monotone.
Proof. Let D = blockdiag(A11 , . . . , Amm ). By (i), A−1
> 0 and by (ii), Aij Ajj
≤ 0, i = j , implying that AD −1 ≤
I . Now Theorem 7 with P = M = I , Q = D −1 can be
applied to prove monotonicity of A if vT AD −1 > 0 or
AD −1 v > 0, which, however, holds by (iii). (Note that we

have assumed strict inequality, A−1
ii > 0.)
The following theorem shows that monotonicity of twoby-two block matrices holds if and only if its Schur complement is monotone.
Theorem 9. Let A be a two-by-two block matrix


A11 A12
A21 A22

where T (ai , bi , ci ) stands for a tridiagonal matrix with
diagonal coefficients bi and off-diagonal ai and ci . Let
B = 2 block tridiag T 0, −bi , i ,

T (−ai , 2(ai + bi ), −ai ), T i , −bi , 0


Clearly, A ≤ B and by Theorem 7 with P = Q = I , M =
B, monotonicity of A follows if the inequality B −1 Ae > 0
and monotonicity of B hold.
The diagonal blocks of B are clearly monotone, and since
the block components (Be)i , i = 1, 2, . . . , n are nonzero
and nonnegative, applying Corollary 1 to B T , we conclude


Finite Difference Methods

that B is monotone if
ci ≤

ai bi+1
ai bi−1
, ci ≤
ai+1 + bi+1
ai−1 + bi−1

By a symmetry argument, B is also monotone if
ci ≤

ai−1 bi
ai+1 bi
, c ≤
ai−1 + bi−1 i
ai+1 + bi+1

To prove the monotonicity of A, it remains to be shown that
B −1 Ae > 0. Clearly, (Ae)i is nonzero and nonnegative for
i = 1, 2, . . . , n, and since the diagonal blocks of B −1 are
positive and B −1 ≥ 0, the required inequality follows.
Note that in the constant coefficient case, the conditions
above take the form



which is stronger than the ellipticity condition c < ( ab).
Note also that in the matrix in (25), where a nine-point
stencil has been used for the approximation of auxx + buyy ,
the difference stencil is of positive type only if |c| ≤
(a + b)/6.

The equation can be solved by a semidiscretization
method, such as the method of lines as it has also been
called. In such a method, one usually begins with discretization of the space derivatives in the equation, which leaves us
with a system of ordinary differential equations in variable
t, that is, an initial value problem. The system is stiff, and
to enable choosing the time-steps solely based on approximation properties one uses stable implicit methods. In
particular, a simple method called the θ-method can be used.
Alternatively, we may begin with discretization in time
using, for instance, the θ-method. This results in an elliptic boundary value problem to be solved at each timestep, which can be done using the methods presented in
Section 3.
Usually, the order in which we perform the discretizations, first in space and then in time, or vice versa, is
irrelevant in the respect that the same algebraic equations,
and hence the same numerical solution, result at each timestep. However, the analysis of the methods may differ. Also,
if we intend to use a variable (adaptive) mesh in space, it is
more natural to begin with the discretization in time. At various time-steps, we can then use different approximations
in space.

4.1 Properties of the solution



In this section, we discuss the numerical solution of
parabolic problems of the form
= Lu + f (x, t),

x ∈ , t > 0


with initial condition u(x, 0) = u0 (x) and given boundary
conditions on ∂, valid for all t > 0. Here, L is an elliptic
operator and f is a given, sufficiently smooth function.
The boundary conditions can be of general, Robin type,
(∂u/∂n) + σ[u − g(x, t)] = 0, where σ ≥ 0, g is given and
n is the outer unit vector normal to ∂. Here, σ = ∞
corresponds to Dirichlet boundary conditions. Frequently,
in applications we have L = , the Laplace operator. As
is well known, this equation is a model for the temperature
distribution in the body , for instance. The equation is
called the heat conduction or diffusion equation.
Stability and uniqueness of the solution of problem (30)
can be shown using a maximum principle, which holds for
nonpositive f , or decay of energy for the homogeneous
problem. Such and other properties of the solution can
be important for the evaluation of the numerical solution

4.1.1 A maximum principle
For ease of exposition, we describe the maximum principle
for an equation of the form (30) in one space variable. On
the other hand, we allow for a domain whose boundary
may depend on time. (Such problems arise in so-called
free boundary value problems, where the boundary between
two matters, such as ice and water, may vary with time.
Frequently, the temperature of ice is assumed to be constant
and it remains to compute the temperature distribution in
the water. Such a problem also arises in connection with
Hence, let the domain D be defined by the indicated parts
of the boundary lines
L0 = {(x, 0) | φ1 (0) ≤ x ≤ φ2 (0)}
L1 = {(x, T ) | φ1 (T ) ≤ x ≤ φ2 (T )},

T >0

and the curves
K1 = {(φ1 (t), t) | 0 ≤ t ≤ T }
K2 = {(φ2 (t), t) | 0 ≤ t ≤ T }
where φ1 (t) < φ2 (t) are continuous functions. Let 0 =
L0 ∪ K1 ∪ K2 and  = 0 ∪ L1 (see Figure 4).

Finite Difference Methods 19
where u(0, t) = u(1, t) = 0 and u(x, 0) = u0 (x) is a given
sufficiently smooth function. We assume that a is a
constant, satisfying a ≤ K − c, where K = ( ux / u )2 ,
c is positive constant and u 2 = 0 u2 dx. Letting E(t) =
1 2
(1/2) 0 u (·, t) dx (the square L2 norm, a measure
of energy) and using (32), we find that 0 ut u dx =
1 2
1 2
− 0 ux dx + a 0 u dx, that is, E  (t) = u 2 (a − K) ≤
−c u 2 = −2cE(t), or



ϕ1(t )


ϕ2(t )

E(t) ≤ e−2ct E(0),

t >0



Figure 4. Domain of definition for a parabolic problem.

Theorem 10. If u is the solution of
∂ 2u
= 2 + f (x, t),

(x, t) ∈ D


where f ≤ 0, and if u is sufficiently smooth in D, then
max(x,t)∈D u(x, t) ≤ max(x,t)∈0 u(x, t), that is, u takes its
maximum value on 0 .

Proof. By contradiction.

The same maximum principle can be proven for problem (30). From this maximum principle, uniqueness and
stability of solutions of (30) follow.
Theorem 11. Let D =  × [0, T ], t > 0.
(a) Problem (30) has at most one solution in D, which
takes prescribed boundary and initial conditions on 0 .
(b) If the boundary and initial conditions are perturbed by
an amount of at most ε, then the solution is perturbed
by at most this amount.
(c) If f in (30) is perturbed by a nonpositive function (but
the boundary and initial conditions are unchanged),
the solution u is also perturbed by a nonpositive function.
As an application of the maximum principle, we consider
the derivation of two-sided bounds of the solution.
Corollary 2. Let Ku ≡ u (t) + Lu = f (t), t > 0 and let
u and u be two sufficiently smooth functions that satisfy
Ku ≤ f ≤ Ku in D, u ≤ u on 0 . Then u ≤ u ≤ u.

4.1.2 Exponential decay of energy
Consider a heat equation with a reaction term
ut = uxx + au,

0 < x < 1, t > 0


with E(0) = (1/2) 0 u20 dx. Hence, E(t) → 0 exponentially, as t → ∞.
The constant K can take arbitrary large values. For
example, for u(x, t) = sin kπx g(t) and g = 0, K = (kπ)2 ,
and here k can be arbitrary large. By the classical Sobolev
inequality, there holds that K ≥ π2 . Hence, if a ≤ 0, then
we can take any c ≤ K, or c = π2 and E(t) ≤ e−2π t E(0),
t > 0.
A similar result holds also for more general parabolic
problems ut + Lu = au, where the operator L is coercive,
for example,  Lu u dx ≥ ρ  u2 dx for some positive
constant ρ.

4.1.3 Exponential decay of the solution
The solution of a parabolic problem depends on initial data
on the whole space  at all times. However, an exponential
decay holds not only for the energy but also for the solution,
away from the support of the initial function, assuming it
has a bounded support. This observation can be based on
the explicit form of the solution of the pure initial value
problem (referred to as the Cauchy problem),
ut = uxx in − ∞ < x < ∞, t > 0
 ∞ (x−y)2
e− 4t u0 (y) dy
u(x, t) = √
4πt −∞


where u(x, 0) = u0 (x), u0 ≡ 0 outside a bounded domain
. From (33), it readily follows that

|u(x, t)| = O |x|−1 e− 4t

as |x| → ∞

which shows a rapid decay away from the region of
compact support of the initial function. Formula (33) also
shows that the solution u(x, t) is an infinitely differentiable
function of x and t for any positive t and that a similar
decay holds for its derivatives.
Remark 3. Since, by a partitioning of unity, any initial
data can be partitioned in packets of initial data with compact support, it can hence be efficient for linear problems


Finite Difference Methods

to compute the solution for each such data packet, even
in parallel, at least for the smallest values of t. Since the
domain of influence widens with O( t), as t increases, we
must, however, add an increasing number of contributions
to the solution from several such subdomains to get the
solution at each point in the solution domain.


Finite difference schemes: the method of

For the numerical solution of parabolic problems, one
commonly uses a semidiscretization method.
Consider first a discretization of space derivatives in (30),
for instance by use of central difference approximations.
Then, for each nodepoint in the mesh h ⊂ , including
the boundary mesh points, where the solution u is not
prescribed, we get an approximation U(t), which is a vector
function of t, and whose ith component approximates u at
xi ∈ h at time t. The vector function U satisfies the system
of ordinary differential equations we get when the operator
L in (30) is replaced by a matrix Ah , corresponding to a
difference approximation of . Hence,
= Ah U(t) + b(t),

exp(tA) −−−→ 0

Let α(A) = maxi Reλi (A) be the so-called spectral abscissa of A. Analysis shows that the system is asymptotically stable for any perturbation of initial data if and
only if α(A) < 0. Similarly, for the solution of (35), there
holds |u(t)| → 0 if α(A) < 0 and |b(t)| → 0 and |u(t)| is
bounded if α(A) < 0 and |b(t)| is bounded.
Similar to the energy estimates, an alternative analysis
can be based on the spectral abscissa β(A) of the symmetric
part of A, that is, β(A) = maxi λi [(1/2)(A + AT )]. We
assume that β(A) < 0.
Then, let E(t) = (1/2) u(t) 2 . It follows from (35) that
E  (t) = (Au, u) + (b(t), u) = 12 [(A + AT )u, u] + (b(t), u)
where (u, v) = uT v. Hence,
E  (t) ≤ β(A)E(t) + 12 |β(A)|E(t) + |β(A)|−1 b(t)

t >0


4.2.1 Stability of the method of lines
We comment first on the stability of linear systems of
ordinary differential equations,
t > 0, u(0) = u0

where A is a n × n matrix. Its solution is
exp[(t − s)A]b(s) ds,
u(t) = exp(tA)u0 +


t >0


that is

E(t) ≤ e 2 β(A)t E(0)
|β(A)|−1 e 2 β(A)(t−s) b(s)




and, since β(A) < 0, if b(t) → 0, t → ∞ then exponential decay of energy follows.
The following relation between α(A) and β(A) holds.
Lemma 6 (Dahlquist (1959), Lozinskij (1958)) −β(−A)
≤ α(A) ≤ β(A).
Proof. If Ax = λx, where Re(λ) = α(A) and x = 1,
then x∗ A∗ = λx∗ , and hence (1/2)x∗ (A∗ +A)x = (1/2)(λ+
λ) = α(A). But, β(A) = max x =1 (1/2)x∗ (A∗ + A)x, by
the Rayleigh quotient theory. Hence α(A) ≤ β(A). Similarly, α(−A) ≤ β(−A). But α(A) ≥ −α(−A), as an elementary consideration shows. Hence α(A) ≥ −β(−A). 
Unlike the scalar case, supt≥0 exp(tA) may be strictly
greater than 1 even though α(A) is negative. The next
relation illustrates further how exp(tA) may grow and
shows exactly when supt≥0 exp(tA) = 1.


The analysis of stability of such systems can be based on
the Jordan canonical form of A; see, for example, Iserles
(1996). Without going into full details, we consider only


E  (t) ≤ 12 β(A)E(t) + |β(A)|−1 b(t)

where b(t) contains the values of f (x, t) at x = xi and any
nonhomogeneous boundary condition. For t = 0, we have
that the ith component of U(0) satisfies Ui (0) = u0 (xi ) at
the meshpoints xi of h .
In general, for elliptic operators, Ah has a complete
eigenvector space with eigenvalues with negative real parts
(see Section 3), and (34) is stable if b is bounded. This
method has been called method of lines because we approximate u along each line perpendicular to  in the time-space
domain and beginning in xi .

= Au + b(t),

the case where the homogeneous system is asymptotically
stable, that is,

Theorem 12.
(a) eα(A)t ≤ exp(tA) ≤ eβ(A)t , t ≥ 0.
(b) supt≥0 exp(tA) = 1 ⇔ β(A) ≤ 0.

Finite Difference Methods 21
Proof. Let v be an arbitrary unit vector, v = 1, and let
φ(t) = v∗ exp(tA∗ ) exp(tA)v. Then,
φ (t) = (exp(tA)v)∗ (A∗ + A) exp(tA)v

sup φ(t) = exp(tA)
v =1


exp(tA) ≥ max |eλi t | = max eReλi t = eα(A)t

which proves (a). By (a), the sufficiency part of (b) is
already proven. The converse follows by
sup exp(tA) = 1 ⇒ φ(t) ≤ 1,

t ≥ 0 ∀v ∈ C n , v = 1


Hence, since φ(0) = 1, we have φ (0) ≤ 0 ∀v, v = 1 or
by (36), v∗ (A∗ + A)v ≤ 0 ∀v, v = 1. Hence, β(A) ≤ 0.

Corollary 3. If A is a normal matrix (i.e. A∗ A = AA∗
that is, A is diagonalizable), then β(A) = α(A) and the
inequalities in Theorem 12 (a) are sharp: exp(tA) =
eα(A)t .

4.3 The θ-method
For the numerical solution of the system of ordinary differential equations, arising from the method of lines, there
exist a plethora of methods. We consider here in more detail
only one such method, the so-called θ-method. It is presented for linear problems but it is also readily applicable
for certain nonlinear problems.
For (34), the θ-method takes the following form
[I − (1 − θ)kA]V(t + k) = (I + θkA)V(t)
+ k[(1 − θ)b(t + k) + θb(t)]


t = 0, k, 2k, . . ., where V(0) = U0 , I is the identity operator and V(t) is the corresponding approximation of U(t).
Here, θ is a method parameter.
The θ-method takes familiar forms for particular values
of θ. For example, θ = 1 yields the familiar Euler forward
V(t + k) = V(t) + k[AV(t) + b(t)]

while θ = (1/2) determines the trapezoidal (or Crank–
Nicolson) rule
A [V(t) + V(t + k)]
+ b(t) + b(t + k)

V(t + k) = V(t) +


This implies the rightmost inequality of (a). To prove the
leftmost inequality, note that for any matrix B, B ≥
maxi |λi (B)|.

V(t + k) = V(t) + k [AV(t + k) + b(t + k)]


that is, by Rayleigh quotient theory, φ (t) ≤ 2β(A)φ(t).
Hence, φ(t) ≤ e2β(A)t and for every t,

θ = 0 yields the backward Euler (or Laasonen) method

We see from (37) that the θ-method is explicit only for θ =
1; for all other values of θ, the method is implicit, requiring
the solution of a linear algebraic system of equations at
each step. The extra computational labor caused by implicit
methods have to be accepted, however, for reasons of
Observe further that I − (1 − θ)kA is nonsingular if k is
small enough or if Reλi ≤ 0 ∀i and θ ≤ 1, where λi is an
eigenvalue of A.
For the stability analysis of time-stepping methods, one
can use the Jordan canonical form of the corresponding
system matrix. For a general linear system of first-order
difference equations
V(t + k) = BV(t) + c(t),

t = 0, k, 2k, . . . , V(0) = U0
It shows that the homogeneous system is asymptotically
stable, that is, B r → 0, r → ∞ if and only if µ(B) < 1,
where µ(B) = maxi |λi (B)|. Further, the solutions of the
inhomogeneous system (38) are bounded if µ(B) < 1 and
|c(rk)| is bounded for r → ∞ and satisfy x(rk) →
0, r → ∞ if µ(B) < 1 and |c(rk)| → 0, r → ∞.
However, even if µ(B) < 1, for nonnormal (i.e. nondiagonalizable) matrices, it may happen that B r takes huge
values for finite values of r. Hence, in practice, we may
have to require that B < 1.
An alternative approach is to use the numerical radius
r(B) of B, where r(B) = max{|x∗ Bx|; x ∈ Cn , (x, x) = 1}.
It is seen that r(B) ≤ B , and it can be shown (see
e.g. Axelsson, 1994) that B ≤ 2r(B). It can further be
shown that if r(B) ≤ 1, then r(B k ) ≤ r(B)k ≤ 1, k = 1, . . .
for any square matrix B. Hence, in general, the stronger
stability condition B < 1 can be replaced by r(B) ≤ 1.
Unfortunately, the computation of the numerical radius is,
in general, complicated. In the case when B is nonnegative,
it can be shown that r(B) = ρ[(1/2)(B + B T )], which can
be used to compute r(B).

4.3.1 Preservation of stability
We show now the conditions when the θ-method preserves
the stability of the system (35), that is, for which it holds


Finite Difference Methods

α(A) < 0 ⇒ µ(B) < 1. We then first put (37) in the form
of (38). Thus, (38) holds with
B = [I − (1 − θ)kA]−1 [I + θkA]
c(t) = k[I − (1 − θ)kA]−1 [(1 − θ)b(t + k) + θb(t)]

where c1 , . . . , cn are determined by the expansion u0 =

j =1 cj vj and {λj , vj } are the eigenpairs of A. The corresponding solution of the difference equation is
V(t) =


cj µtk vj ,

t = 0, k, 2k, . . .

j =1

Let (λj , vj ), j = 1, 2, . . . , n, denote the eigensolutions of
A with the ordering
Re(λ1 ) ≥ Re(λ2 ) ≥ · · · ≥ Re(λn )

µ(β; θ) =

Then, the eigensolutions of B are (µj , vj ), j = 1, 2, . . . , n,
1 + θkλj
µj =
1 − (1 − θ)kλj
Theorem 13. Assume that (35) is asymptotically stable
(Reλ1 < 0) and that b(t) is bounded, all t ≥ 0. Then
(a) (38), with 0 ≤ θ ≤ (1/2), is asymptotically stable ∀k
> 0. (Unconditional stability.)
(b) (38), with (1/2) < θ ≤ 1, is asymptotically stable if
and only if


2Re(λj )
(2θ − 1)|λj |2

j = 1, 2, . . . ,

where µj = µ(kλj ; θ) and where we have introduced the

(Conditional stability) (40)

1 + θβ
1 − (1 − θ)β

−∞ < β < ∞

Hence, µj is the damping factor corresponding to eλj t (note
Reλj < 0), and it is important that each factor is sufficiently
small, even for the large eigenvalues. There holds

lim µ(β; θ) =




Hence, |µ| is less than one only if θ < 1/2. In particular, for
θ = 1/2, it holds that µ → −1, as β → −∞. Therefore, the
Crank–Nicolson method can have an insufficient damping,
especially at the initial, transient, phase (small t) where all
eigenvector components may have a significant value. It is
then better to choose θ < 1/2, for instance, θ = (1/2) − ζk
for some ζ > 0. (The latter choice will preserve the second
order of discretization error, as is shown below.) Otherwise,
for a fixed value of θ < 1/2, the time-discretization error
is only O(k).

Proof. An easy calculation establishes that
|µj |2 = µj µj =

4.3.2 Discretization error estimates for the θ-method

1 + 2θk Re(λj ) + θ2 k 2 |λj |2
1 − 2(1 − θ)k Re(λj ) + (1 − θ)2 k 2 |λj |2

where µj is given by (39). Hence,
|µj |2 < 1 ⇔ k(2θ − 1)|λj |2 < −2Re(λj )


The assumption that (35) is asymptotically stable, implies
that Re(λj ) < 0 for j = 1, 2, . . . , n. It is easy to see that
the inequality on the right side of (41) is satisfied for j =
1, 2, . . . , n, precisely under the conditions on k presented

The solution of the homogeneous system u (t) = Au(t),
t > 0, u(0) = u0 , where A has a complete eigenvector
space, can be written in the form
u(t) =


j =1

cj etλj vj ,

t >0

The discretization error estimates for the θ-method can
be readily derived. The fully discretized scheme takes the
[I − k(1 − θ)Lh ]v(x, t + k) = [I + kθLh ]v(x, t) + k(1 − θ)
× f (x, t + k) + kθf (x, t),

t = 0, k, 2k, . . .


We have f = ut − Lu and


u(x, t + k) − u(x, t) =

ut (x, s) ds

Substituting the differential equation solution in (42), we
then obtain that
[I −k(1−θ)Lh ]u(x, t+k) = [I +kθLh ]u(x, t)+k(1−θ)
ut (x, s) ds−k(1−θ)
× f (x, t+k) + kθf (x, t) +

Finite Difference Methods 23

× Lh u(x, t+k)−kθLh u(x, t) − k(1−θ)[ut (x, t+k)

− Lu(x, t+k)] − kθ[ut (x, t)−Lu(x, t)]
Letting e(x, t) = u(x, t) − v(x, t), we then find
[I − k(1 − θ)Lh ]e(x, t + k) = [I + kθLh ]e(x, t)
+ τ(x, t; h, k)
where τ is the truncation error
ut (x, s) ds
τ(x, t; h, k) =

− [k(1 − θ)ut (x, t + k) + kθut (x, t)]

In applying this method to a nonstationary partial differential equation problem, we first keep k constant and run the
problem over a time-step with a mesh in space with the
following mesh sizes: say h and (1/2)h (independently of
each other). Then, as usual, (1/3)[4vh/2 (xi , t) − vh (xi , t)]
provides an estimate of the space-discretization error, provided that this is O(h2 ), h → 0.
Similarly, we can let h be fixed and run the numerical
method with stepsize k and two steps with k/2. In the
same manner as above, we can then estimate the timediscretization error.
If the solution is smooth, these estimates are frequently
accurate. However, in an initial transient phase, for instance,
when the solution is less smooth, the estimates are less

+ k(1 − θ)(L − Lh )u(x, t + k) + kθ(L − Lh )u(x, t)
Note that the truncation error consists of two parts: (1)
the time-discretization error and (2) the space-discretization
error. The stability analysis shows that they can be estimated independently of each other and the discretization
error in L2 norm in space becomes
u(·, t) − v(·, t)

≤ ) θ − 12 O(k)) sup utt + O(k 2 ) sup uttt

+ O(h )(



4.4 Nonlinear parabolic problems
We consider nonlinear evolution equations of parabolic
type, that is, for which an asymptotic stability property holds.
The stability and discretization error for infinite time for
the θ-method can be analyzed for nonlinear problems
+ F (t, u) = 0,




h, k → 0


If we choose θ = (1/2) − ζk for some positive constant ζ,
then the total discretization error becomes O(k 2 ) + O(h2 ).
Full details of such an analysis and applications for nonlinear problems of parabolic type can be found in Axelsson
Remark 4. Note that we were able to derive this estimate without the use of the monotonicity of (−h ) (cf.
Section 3). Hence, (43) is valid even for nonmonotone
approximations and for more general operators of second
order. In particular, it is valid for central difference approximations Lh of the convection–diffusion operators (6, 8),
as long as h is small enough so that the real part of the
eigenvalues of Lh is negative. This is due to the property
of the θ-method for 0 ≤ θ < (1/2) that it is stable for operators with arbitrary spectrum in the left-half plane of the
complex plane. Methods with such a property are called
A-stable methods; see Dahlquist (1963). It is known that
multistep time-discretization methods that are A-stable can
be at most of second-order of approximation. On the other
hand, certain implicit Runge–Kutta methods can have arbitrary high order and still be A-stable; see Axelsson (1964).
Remark 5. A simple method to estimate the discretization
errors numerically is provided by Richardson extrapolation.

t > 0, u(0) = u0 ∈ V


where V is a reflexive Banach space and F : V → V  , where
V  denotes the space dual to V . We have V → H → V 
for some Hilbert space H , where → denotes continuous
injections. Let (·, ·) be the scalar product in H and ·
the associated norm. We assume that F (·, t) is strongly
monotone (or dissipative), that is,
(F (t, u) − F (t, v), u − v) ≥ ρ(t) u − v



for all u, v ∈ V , t > 0, where ρ ≥ ρ0 > 0. It follows then
1 d
( u − v 2 ) = −(F (t, u) − F (t, v), u − v)
2 dt
≤ −ρ0 u − v 2 ,

t >0

u(t) − v(t) ≤ exp(−ρ0 t) u(0) − v(0)
For such nonlinear problems of parabolic type, the stability
and discretization error for the θ-method can be analyzed
(cf. Axelsson, 1984). Here, only a one-sided Lipschitz
constant enters in the analysis. The accuracy of the θmethod, however, is limited to second order, at most.
Here, we consider an alternative technique, which allows


Finite Difference Methods

an arbitrary high order of approximation. It is based on a
boundary value technique.
The method is presented in a finite dimensional space
V = Rn and on a bounded time interval. Given a system
of differential equations
+ F (t, u) = 0,

0 < t ≤ T , u(0) = u0 prescribed

we first make a transformation of the equations to a more
suitable form. In many problems, there may exist positive
 and of
parameters εi , 0 < εi ≤ 1, such that parts of F
the corresponding parts of the Jacobian matrix ∂ F
are unbounded as O(ε−1
corresponding equation by this parameter to get

+ F (t, u) = 0,

0<t ≤T


where ε is a diagonal matrix with entries εi , and now it
is assumed that F and (∂F /∂u) are bounded with respect
to ε.
From the monotonicity of F for t ≥ t0 > 0 follows
1 d
(ε(u − v), u − v) = −(F (t, u) − F (t, v), u − v)
2 dt
≤ −ρ(t) u − v


≤ −ρ(t)(ε(u − v, u − v)),

t ≥ t0


u(t) − v(t)


≤ exp


−2ρ(s) ds


≤ u(t0 ) − v(t0 ) 2ε ,

u(t0 ) − v(t0 )


t0 ≤ t ≤ T

where v ε = (εv, v)1/2 . This means that the system is contractive for t ≥ t0 . In the initial phase t ∈ (0, t0 ), the system
does not have to be contractive, that is, the eigenvalues of
the Jacobian may have positive real parts there. In this interval, we may choose a step-by-step method with sufficiently
small step sizes, in order to preserve stability and to follow
the transients of the solution.

4.4.1 Boundary value techniques for initial value
The traditional numerical integration methods to solve initial value problems
= f (t, u(t)),

t > 0, u(0) = u0 given


are step-by-step methods. Some such familiar methods are
Runge–Kutta and linear multistep methods. In step-by-step
methods, the error at the current time-step depends on the

local error at this step and also on the errors at all previous
time-steps. In this way, the errors accumulate and the total
error is typically larger by a factor O(k −1 ) than the local
errors. A direct global error control cannot be justified since
the global error depends on the stability of the problem and
the errors at the previous time-points.
In boundary value methods, on the other hand, all approximations at all time-points are computed simultaneously
and such methods may be coined global integration methods. By its very nature, a boundary value method is better
adapted for global error control.
For problems where one envisions that the solution suddenly may become unsmooth, one can implement boundary
value methods as time-stepping methods but with much
larger time-steps than for standard step-by-step methods.
A boundary value method can be composed of a sequence
of forward step-by-step methods followed by a stabilizing backward-step method. As a simple example, let u0
be given and
un+1 − un−1 − 2kf (tn , un ) = 0,
u −u


n = 1, 2, . . . , N − 1

− kf (tN , u ) = 0

whose solution components u1 , u2 , . . . , uN must be computed simultaneously. Such a method was analyzed in
Axelsson and Verwer (1984). For a more recent presentation of boundary value and other methods, see Brugnano
and Trigiante (1998).
The Galerkin method. For the interval (t0 , T ), the following global Galerkin method can be used. The interval is divided into a number of subintervals (ti−1 , ti ), i =
1, 2, . . . , N , where tN = T . The length of the intervals,
ti − ti−1 , may vary smoothly with some function h(ti ),
but for ease of presentation, we assume that the intervals have the same length ti − ti−1 = h, i = 1, 2, . . . , N .
Consider each interval as an element on which we place
some extra nodal points, ti,j , j = 0, 1, . . . , p, such that
ti,j = ti + ξj h, where ξj are the Lobatto quadrature points
satisfying 0 = ξ0 < ξ1 < · · · < ξp = 1 and ξj + ξp−j = 1.
Hence, the endpoints of the interval are always nodal points
and (if p > 1) we also choose p − 1 disjoint nodal points
in the interior of each element.
To each nodal point, we associate a basis function φi,j .
The basis functions may be exponential or trigonometric
functions, and may also be discontinuous, but here we
consider the case that they are continuous and polynomial
over each element. Basis functions corresponding to interior nodes have support only in the element to which they
belong, and those corresponding to endpoints have a support over two adjacent elements (except those at t0 and tN ).
The number of nodal points in each closed interval then
equals the degree of the polynomial p plus one.

Finite Difference Methods 25

V = ϕi,j

Let S h be the subspace of test functions that are zero at
t0 , that is,

j = 0, i = 0, 1, . . . , N − 1,

r = 1, 2, . . . , m



S h = span{φi,j , i = 0, 1, 2, . . . , N − 1, j = 1, 2, . . . , p}


+ F (t, U ), V
U, V ∈ H 1 (t0 , T )

a(U ; V ) ≡



where H 1 (t0 , T ) is the first-order Sobolev space of functions with square-integrable derivatives. To get an approx of the solution of (46), we take a vectorial test
imation U
function V = ϕi,j
, multiply the equation, and after integration, we obtain
 ti+1 * 
; ϕ ) =
), ϕ
dt = 0,
+ F (t, U
i = 1, . . . , N − 1
j =0
 ti+1 * 
; ϕ ) =
), ϕ
+ F (t, U
dt = 0,
i = 0, 1, . . . , N − 1
j = 1, 2, . . . , p − 1
and at tN = T , we get
 tN * 
; ϕ[r] ) =
), ϕ[r] dt = 0
+ F (t, U





corresponding to S h , where
di,j φi,j ,

di,j ∈ Rm

i=0 j =1

that is, we have imposed the essential boundary conditions
at t0 . Clearly,
a(U ; V ) = 0

∀V ∈ [H 1 (t0 , T )]m

From (48), we obtain
; V ) =
a(U ; V ) − a(U



% *




) , V dt = 0
+ F (t, U ) − F (t, U




U − UI


dt ≤ C0 h2(p+1)

, dU
dUI ,
, dt − dt , dt ≤ C1 h



p+1 dt

, dU ,2
dt (52)
, dt ,

T p+1
= t0 k=0 [(∂ k U/∂t k ), (∂ k U/∂t k )] dt is the
Here, U p+1
norm in the Sobolev space H p+1 (t0 , T ).

(F (t, U ) − F (t, V ), U − V ) ≥ ρ(t) U − V

= ϕi,j er , where er is the rth coorWe choose in turn

dinate vector. This defines the Galerkin approximation U


Theorem 14. Let U be the solution of (46), and conditions


 = U (t )φ +
0 0,0

and similarly for (49) and (50).
To estimate the Galerkin discretization error U − U
we let UI ∈ Sh be the interpolant to U on {ti,j }, j =
 = (U −
0, 1, . . . , N − 1. From the representation U − U
) = η − θ, we see that θ = U − U
 ∈S .
UI ) + (UI − U
Assuming that the solution U is sufficiently smooth, from
the interpolation error expansion in integral form, we obtain
the usual Sobolev norm estimates for the interpolation error:


F (t, U ) − F (t, V ) ≤ C U − V , t > 0
, in the space
be satisfied. Then the Galerkin solution U
of piecewise polynomial continuous functions of degree p,
defined by (48–50) satisfies
 | = O(hp+ν )
|U − U




+ U





where ν = 1 if p = 1, 1 ≥ ν ≥ (1/2) if p = 3, 5, . . . and
ν = 0 if p is even, and
|V |2 =

(εV (T ), V (T )) +


ρ(t) V (t)




Proof. For a detailed proof, see the supplement of (Axels
son and Verwer, 1984).
We have assumed that F is Lipschitz-continuous, that is,
F (t, u) − F (t, v) ≤ C u − v for all u, v ∈ Rm . This is,
however, a severe limitation as it is a two-sided bound.
Difference schemes. In order to get a fully discretized
scheme, one has to use numerical quadrature, which results
in various difference schemes. We consider this only for
p = 1. Then, ϕi,p = φi are the usual hat functions and there

 = U (t )φ + N U φ ,
are no interior nodes. With U
0 0
i=1 i I


Finite Difference Methods

relations (48), (50) imply
 φ +U
F (t, U
i−1 i−1
i i
 φ )φ dt
i+1 i+1 i

 ε(U
F (t, U
N−1 N−1 + UN φN )φN dt


We call this the generalized
midpoint rule difference
 )) . If we use numerical intescheme. Let Fi = F (t, U
i t=ti
gration by the trapezoidal rule, that is,
F φi dt = [Fi−1 φi (ti−1 ) + Fi φi (ti )] = hFi
which is known to be of O(h2 ) accuracy, on the basis of
(53), we may derive a more accurate difference scheme.
For this purpose, let
F (t) 

(Fi−1 + Fi ) + t − ti +
ti−1 ≤ t ≤ ti

(F − Fi−1 ),
h i

except that for the last formula in (53), we use F (t) 
(1/2)(FN−1 + FN ), tN−1 ≤ t ≤ tN . Then,
F φi dt  (Fi−1 + Fi ) + (Fi − Fi−1 )

(F + 2Fi ),
6 i−1

i = 1, 2, . . . , N − 1

and similarly



F φi dt 

(F + 2Fi )
6 i+1

Hence, the generalized midpoint rule (53) takes the form

= Fi−1 + 4Fi + Fi+1 ,
i = 1, 2, . . . , N − 1,


+ FN
ε UN − U
N−1 =
2 N−1
We notice that this is a combination of the Simpson and
trapezoidal rules.
For this combination, numerical tests in Axelsson and
Verwer (1984) indicate very accurate results. It is found that
already on a very coarse mesh (h = 1/4), the accuracy is
high. For (du/dt) = δu and δ < 0, the order of convergence
seems to be  3.5.
The above method is related to certain types of implicit
Runge–Kutta methods; see, for example, Axelsson (1964)
and Butcher (1977). As such, they are A-stable. The global
coupling preserves this stability.

4.5 Computational aspects
The θ-method, θ = 1 and other unconditionally stable methods are implicit, that is, they give rise to a linear system of
algebraic equations to be solved at each time-step. For this
purpose, one can use direct or iterative solution methods.
Hereby, the associated cost is important, in particular as
there will be a large number (O(k −1 )) of such systems to
be solved.

4.5.1 Iterative solution methods. The conjugate
gradient method
The matrix in the linear system, which arises in the θmethod and has the form I − (1 − θ)kA. If A is symmetric
and negative definite, which is typical for parabolic problems, then the system has a condition number

1 + (1 − θ)kb
1 + (1 − θ)ka

where a = − maxj λj and b = − minj λj and λi are the
eigenvalues of A. Hence, the condition number is bounded
by κ < 1 + (1 − θ)kb, that is (if a ≤ 1), the condition
number of A is typically multiplied by (1 − θ)k/a, which
can significantly decrease its value. For difference methods
for second-order operators L, it holds b = O(h−2 ), so the
number of iterations using the standard unpreconditioned
conjugate gradient method grows only as O(h1/2 ), if k =
O(h). Hence, the arising systems can normally be solved
with an acceptable expense, in particular if one, in addition,
uses a proper preconditioning of the matrix.
For stability reasons, for an explicit method, one must
let k = O(h2 ). Hence, there are O(h−1 ) more time-steps,
so explicit methods are generally more costly except where
there is a need to choose small time-steps of O(h2 ) due to
approximation accuracy reasons. However, for k = O(h2 ),
this is balanced by the condition number κ = O(1), that
is, the number of iterations are bounded irrespective of h,
which make implicit methods still preferable because of
their robust stability.
For problems with a small condition number, approximate inverse preconditioners can be quite accurate; see,
for example, Axelsson (1994) for references to such methods. This can make implicit methods behave essentially as
explicit but with no time-step stability condition. An example is the Euler backward method, where the arising matrix,
I − kA, can be approximated by (I − kA)−1 ≈ I + kA, the
first term in the Neumann expansion. This is equivalent to
the Euler forward method. Adding more terms increases the
accuracy of this approximation if kA < 1.
Frequently, it can be efficient to reuse the same preconditioner for several time-steps, whence its construction cost
can be amortized.

Finite Difference Methods 27

4.5.2 Periodic forcing functions

4.5.4 Alternating direction implicit methods

For problems

The alternating direction implicit method (ADI); see, for
example, Peaceman and Rachford (1955) and the similar
fractional step method (Yanenko, 1971) can sometimes be
used to solve difference equations for elliptic problems
more efficiently than some standard methods can. However,
they can also be seen as special time-splitting methods using
a direct solution method for the arising systems, which is
the aspect of the ADI methods to be dealt with here.
For notational simplicity, consider a homogeneous

ut = u + f (x, t),

x ∈ , t > 0

where the forcing function is periodic in time, that is,
f (x, t) = f0 (x)eiωt , one can apply the Ansatz
u(x, t) = v(x, t)eiωt
where u and v are complex-valued functions. We find then,