Principal
Encyclopedia of Computational Mechanics (3 Volume Set)
Encyclopedia of Computational Mechanics (3 Volume Set)
Erwin Stein, René de Borst, Thomas J.R. Hughes, Editors
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, geoenvironmental modelling, biomechanics, electromagnetism, metal forming, to name but a few.This three volume set provides the most comprehensive and uptodate collection of knowledge about this increasingly important area of engineering science. The Encyclopedia provides a wellrounded 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 indepth 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, indepth 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.
Categories:
Año:
2007
Edición:
1
Editorial:
Wiley
Idioma:
english
Páginas:
2244
ISBN 10:
0470091355
ISBN 13:
9780470091357
File:
PDF, 57.87 MB
The file will be sent to your email address. It may take up to 15 minutes before you receive it.
The file will be sent to your Kindle account. It may takes up to 15 minutes before you received it.
Please note you need to add our email km@bookmail.org to approved email addresses. Read more.
Please note you need to add our email km@bookmail.org to approved email addresses. Read more.
You may be interested in


Most frequently terms
finite^{4571}
equations^{3252}
boundary^{3118}
equation^{3049}
finite element^{2637}
mesh^{2143}
numerical^{1998}
computational^{1842}
linear^{1817}
matrix^{1626}
eng^{1391}
domain^{1369}
comput^{1347}
approximation^{1332}
nonlinear^{1330}
formulation^{1318}
fluid^{1291}
dimensional^{1181}
mech^{1162}
elastic^{1154}
discretization^{1146}
mechanics^{1043}
discrete^{1035}
schemes^{1016}
vector^{1015}
velocity^{1010}
numer^{987}
appl^{971}
galerkin^{895}
corresponding^{886}
tensor^{873}
operator^{870}
stability^{863}
simulation^{859}
integral^{858}
solutions^{829}
deformation^{819}
integration^{801}
displacement^{787}
int^{785}
algorithm^{775}
flows^{754}
residual^{736}
modeling^{732}
dynamics^{730}
math^{726}
convergence^{714}
variables^{709}
constitutive^{659}
finite elements^{645}
applications^{643}
explicit^{641}
adaptive^{634}
finite element methods^{632}
parameter^{631}
methods appl^{619}
algorithms^{616}
interpolation^{616}
boundary conditions^{615}
stokes^{601}
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.
1

2

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 1 1 2 6 1 MOTIVATION AND SCOPE 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 firstclass uptodate representations of all major computeroriented 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 selfcontrolling 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: 0470846992. 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. 2 STAGES OF DEVELOPMENT AND FEATURES OF COMPUTATIONAL MECHANICS One can say that we are now in about the third period of development of computerbased 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 boundaryvalue 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, 2 Fundamentals: Introduction and Survey such as boundaryvalue 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 problemoriented discretization methods. These new integral methods of weighted residuals are characterized by two important properties: – – – – – – – – (a) (b) Methodical width: This is the intended simple logical and algorithmic structure (e.g. with symmetry properties and wellposedness 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 boxtype 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 technology. 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 follows: – Meshfree and particle methods, finite elements with discontinuities for damage and fracture. Errorcontrolled adaptive modeling and approximation of physical events near to nature, also scalebridging 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 waveletbased 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, 3Danimation 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 problemoriented iterative algebraic solvers and preconditioners with advanced data management for highend 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 considerably. 3 SURVEY OF THE CHAPTERS OF VOLUME 1 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 timedependent 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 hversion of the usually used 2D and 3Dfinite elements. The following chapter, Finite element methods, by Susanne Brenner and Carsten Carstensen, treats the displacement method (primal finite element method) for boundaryvalue problems of secondorder elliptic PDE’s as well as a priori and a posteriori error estimates of the weak solutions and related hadaptivity, including nonconforming elements and algorithmic aspects. This basic chapter is followed by The pversion 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 hrefinement in subdomains with large gradients of the solution are advantageous against the hversion. 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 porders, 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) hpadaptivity for 3Dsystems as well as (v) the efficient implementation of anisotropic pextensions 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 ptype 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 highorder 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 timedependent 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 3D theory of elasticity to 2D theories of thinwalled 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) pextensions 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 saddlepoint 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 problemdependent inf–sup condition. Chapter 10, this Volume, by Antonio Huerta, Ted Belytschko, Sonia FernándezMéndez, and Timon Rabczuk, provides an advanced and systematic representation of different versions and alternatives of the socalled meshfree and particle methods, known as moving least squares, partition of unity FEM, corrected gradient methods, particleincell methods 4 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–Galerkintype 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 errorcontrolled approximation of the essential boundary conditions of a boundaryvalue 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 boundaryvalue 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 GalerkinBEMs for elliptic boundaryvalue 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 problemdependent preconditioners. Special features are Signorinitype contact problems using both primal and dualmixed finite element approximations. Recent features are adaptive hpmethods. There seems to be a lack of available software for 2D – and even more for 3D – problems. 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, fluidsolid interaction – by suitable coordinates and metrics for each of the constituents – for example, Lagrangian coordinates for solids and Eulerian coordinates for fluids – using the wellknown tangential pushforward and pullback 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 ALEconcept, 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 vertexcentered 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 FranzErich 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 boundarybased (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 selfcontrolled 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 freeform 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 PDEsolutions. 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 problemoriented 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, BiCGSTAB, and GMRES. Special eigenvalue problems of a Hermitian matrix are analyzed mainly with iterative QRmethods, 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 Wcycle 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 sideproduct of adaptive mesh refinements. 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 farfield expansion and the panel cluster tree. A generalized variant is the construction of hierarchical matrices (Hmatrices) with different matrix–vector and matrix–matrix 6 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 Timedependent problems. 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 timestepping methods for some classical parabolic equations, like the instationary heat equation and the reactiondiffusion 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, Laplacetransform methods, and timestepping 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 timeharmonic Maxwell equations. With the stabilized variational formulations and Nedelec’s three fundamental elements, hpdiscretizations and hpadaptivity 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 nullspace, in addition to projectionbased 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. 4 WHAT WE DO EXPECT 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 knowledge. Chapter 2 Finite Difference Methods Owe Axelsson University of Nijmegen, Nijmegen, The Netherlands 1 Introduction 2 Twopoint Boundary Value Problems 3 Finite Difference Methods for Elliptic Problems 4 Finite Difference Methods for Parabolic Problems 5 Finite Difference Methods for Hyperbolic Problems 6 Convection–Diffusion Problems 7 A Summary of Difference Schemes References Further Reading 7 9 12 18 28 36 50 52 53 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 symmetric. 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 posed. Definition 1. A boundary value problem is well posed if 1 INTRODUCTION 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: 0470846992. (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 illposed 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. 8 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 ∂ (1) 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 outwardpointing normal vector to . For the Poisson problem −u = ρ(x), x∈R 3 where we ignore the influence of boundary data, the solution is 1 x − ξ−1 ρ(ξ) dξ (2) u(x) = 2π 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 firstorder hyperbolic problem aux + buy = 0, x > 0, y > 0 (3) 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 a a x≥ y h x − by b u(x, y) = b a g y − x x< y a b Therefore, the solution is constant along the lines bx − ay = const, called characteristic lines, so the solution at any point depends just on a singleboundary 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 0 1 vt v dx + or, by letting E(t) = we find 1 0 1 0 avx v dx = 0 v 2 dx (called the energy integral), a 1 d E(t) + v 2 (1, t) = 0, 2 dt 2 t>0 that is, E (t) ≤ 0. Hence, E(t) ≤ E(0), t > 0, where 1 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 wellposedness 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) := 1 [u(x + h) − u(x)], h the forward difference (4) 1 [u(x) − u(x − h)], h the backward difference (5) 1 [u(x + h) − u(x − h)], 2h the firstorder central difference (6) Finite Difference Methods 9 Note that Dx0 u = (1/2)(Dx+ + Dx− )u. More generally, one can use the socalled “θmethod” u (x) [θDx+ + (1 − θ)Dx− ]u(x) (7) where θ is a method parameter. Thus, u (x) = [θDx+ + (1 − θ)Dx− ]u(x) + − h (1 − 2θ)u (x) 2 h2 [(1 − θ)u(3) (x − η2 h) + θu(3) (x + η1 h)] 6 where 0 < ηi < 1, i = 1, 2, so if 1 − 2θ = O(h) otherwise O(h2 ) O(h) 1 [θu(x + h) + (1 − 2θ)u(x) − (1 − θ)u(x − h)] h h2 h + (1 − 2θ)u (x) − u(3) (x + η3 h), −1 < η3 < 1 2 6 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 as 0 u(x), u (x) Dx+ Dx− u(x) ≡ Dxx the central difference of second order (8) 0 u(x) = Dx+ Dx− u(x) = h−2 [u(x + h) − We have then Dxx 2u(x) + u(x − h)], thus 0 Dxx u(x) h2 − u(4) (x + η4 h), 12 x 0 Dxx u(x) + O(h ), u (x) = or u (x) = 2 + u(x − h, y) + u(x, y + h) + u(x, y − h) − 4u(x, y)] (5) is called the 5point 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) 1 (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 (10) 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. 2 TWOPOINT BOUNDARY VALUE PROBLEMS 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 twopoint 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) = α (12) h→0 r1 (u) ≡ γ1 u(b) + δ1 k(b)u (b) = β Similar expressions hold for Dy+ , Dy− , Dy0 , and so on. In (4) 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 (9) 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 realvalued functions and γ0 , δ0 , γ1 , δ1 , α, β are given real numbers. The b b operator L is selfadjoint, 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, α), 10 Finite Difference Methods a problem: b x 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 heatconducting 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 solution. 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. 2.1 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 differences −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) = β (14) 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 ) 0≤i≤N+1 i=0,N+1 that is, the largest value of vh is taken at one of the endpoints. 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 (15) 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 principle. In order to estimate eh , it is convenient to use the following notations: vh πh = max vh (xi ), 1≤i≤N i=0,N+1 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 (17) 1 2 (4) h u (ξi ), 12 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 b−a 2 (18) 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 8 (b − a)2 (b − a)2 Lh u − f πh = τh πh = 8 8 (b − a)2 h2 ≤ max u(4) (ξ) 8 12 u − uh πh ∪∂πh ≤ u − uh ∂πh + vh ∂πh = max vh (xi ) (16) 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 1 [(b − a)h]2 max u(4) (ξ) a<ξ<b 96 (19) 1 (4) u , 12 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 proof. 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 higherorder 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. 12 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 2.3 Computing highorder 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 ydirection 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 fivepoint difference approximation, uh (xi+1 ) − uh (xi−1 ) + O(h2 ) 2h 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 higherorder 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 hexpansion of the errors. 3 x, y h1 x − h1, y u(xi+1 ) − u(xi−1 ) uh (xi+1 ) − uh (xi−1 ) = 2h 2h ϕ(x ) − ϕ(x ) i+1 i−1 + O(h3 ) = u (xi ) − h2 2h h2 + u(3) (η1 ) − h2 ϕ (xi ) + O(h3 ) 6 where η1 ∈ (xi−1 , xi+1 ). h2 (a) N hN W FINITE DIFFERENCE METHODS FOR ELLIPTIC PROBLEMS 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 E P hE hW hS Γ S (b) Figure 2. Local difference meshes (a) rectangular domain, (b) curved boundary. Finite Difference Methods 13 + − + − (5) 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 1 + h−2 2 [uh (x, y − h2 ) − 2uh (x, y) + uh (x, y + h2 )] = f (x, y), x, y ∈ h (20) If we order the meshpoints in lexicographic (horizontal) order, the corresponding matrix takes the form Ah = block tridiag(Ii , A(i) , Ii ) −2 where A(i) = h−2 1 tridiag (1, −2, 1) − 2h2 I and Ii = (i) and the identity matrix I h−2 2 I , i = 1, 2, . . . , N2 . Here, A have order N1 × N1 and there are N2 blockrows 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) 1 h1 ux (x + ξ1 , y) + h22 u(4) Lh u − f = y (x, y + ξ2 ) 12 12 = O(h ), 2 −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 fivepoint difference approximation with, in general, nonequidistant meshlengths (see Shortley and Weller, 1938): 1 1 [u (E) − uh (P )] − [u (P ) − uh (W )] Lh uh := hE h hW h 2 1 × + [u (N ) − uh (P )] hW + hE hN h 2 1 = f (P ), P ∈ h − [uh (P ) − uh (S)] hS hS + hN (21) 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 ydirections, Lh uh := 1 {h u (W ) + hW uh (E) − (hE + hW )uh (P )} hE E h 1 {h u (S) + hS uh (N ) − (hN + hS )uh (P )} hN N h h + hW h + hN 1 hW E f (P ), P ∈ h + hS S = 2 2 2 (22) 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 boundaryfitted 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 Higherorder schemes 3.2.1 The ninepoint difference scheme Let u = f , and consider first the crossdirected fivepoint difference scheme on a local equidistant square submesh, = 12 h−2 [uh (x − h, y − h) + uh (x + h, y − h) (5,×) h + uh (x − h, y + h) + uh (x + h, y + h) − 4uh (x, y)], x, y ∈ h (23) It is readily seen that, for a sufficiently smooth function u, (5) h = u + (5,×) h 2 2 (4) 2 (6) + h4 u(6) h ux + u(4) y x + uy 4! 6! + O(h6 ) 2 2 (4) = u + h2 u(4) + h4 x + 6uxxyy + uy 4! 6! (6) × ux + 15uxxxxyy + 15uxxyyyy + u(6) + O(h6 ) y Let (9) h be the ninepoint difference scheme defined by (9) h = 2 (5) 1 (5,×) + h 3 h 3 14 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 mesh, (9) h uh 1 1 4 2 = f + h2 f + h f + 2fxxyy + O(h6 ) 12 360 where 2 = (f ). Using a modified righthand side in the difference formula, it follows that the difference approximation h2 (5) u = I + f, (9) h h 12 h (x, y) ∈ h (24) has truncation error O(h4 ). Further, it follows from (24) that for a sufficiently smooth 2 2 4 function f , f = (9) h f − (1/12)h f + O(h ). A 2 (5,×) (5) f − f]+ computation shows that h fxxyy = 2[ 4 O(h ) and, therefore, the ninepoint stencil with the next modified righthand side (9) h uh = f + + 1 2 (9) 1 4 (5) h h f − h h f 12 240 1 4 h fxxyy , 180 (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 farfield 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 ninepoint 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 derivative 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 ≈ 1 [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 ninepoint difference stencil, it becomes a+b 1 −2 h 6 − 3c 2 5a − b a+b + 3c 2 5b − a −10(a + b) 5b − a a+b + 3c 2 5a − b (25) a+b − 3c 2 3.2.3 Difference schemes for other regular tessellations Finite differences can be extended to nonrectangular meshes. For a regular (isosceles) triangular mesh, one can form the obvious sevenpoint difference stencil. For a hexagonal (‘honeycomb’) mesh, one finds a fourpoint stencil 4/3 4 4/3 4/3 The symmetrically located nodepoints in the sevenpoint scheme allows one to readily approximate secondorder cross derivatives. Similarly, in 3D problems, a cuboctahedral stencil involves 13 nodepoints. If a Cartesian grid is used, approximations of the secondorder 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 12point stencil for a regular equidistant mesh, 1 2 −8 2 −4 (26) h 1 −8 20 −8 1 2 −8 2 1 Finite Difference Methods 15 which has truncation error O(h2 ). The biharmonic problem can, however, more efficiently be solved as a coupled problem y u=g ξ − u = 0 ξ = f δu =0 δn using a variational formulation hereof; see, for example, Axelsson and Gustafsson (1979) and Ciarlet and Raviart (1974). u=g u=g 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 positivetype 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 estimates 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 ) ≡ n 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 u=g x (a) NW U T W N x x hN P αW xE βW R Qx hs Γ βE αE S SE (b) 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, L−1 h ∞ ≤ 1 α 16 Finite Difference Methods where 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 estimate, eh ∞ ≡ u − uh ∞ 1 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 error. 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 positivetype 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 Mmatrix. 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 α. Then, (a) A−1 ∞ ≤ 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. (b) 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 A−1 ∞ ≤ 1 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 −A21 A22 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 A A 12 ∞ 21 ∞ 1+ A−1 ∞ ≤ 1 + min(A11 v1 )i min(A11 v1 )i i i 1 1 × max , (27) min(A11 v1 )i min(q + A21 A11 p)i i 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 twobytwo 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 conditions. 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 nonsingular. −1 (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 11 A12 is monotone. (b) Conversely, assume that A is monotone. Then S is monotone. 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, % −1 A = −1 −1 −1 A−1 11 + A11 A12 S A21 A11 −1 A−1 11 A12 S S −1 A21 A−1 11 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. problem Consider the mixed derivative elliptic −auxx − 2cuxy − buyy = f in = [0, 1]2 (28) 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 ninepoint secondorder accurate finite difference approximation of this problem on a uniform mesh, yields the block tridiagonal n2 × n2 matrix % 1 A = 2 block tridiag T h − ci c , −bi , i , 2 2 T (−ai , 2(ai + bi ), −ai ), T c c − i , −bi , − i 2 2 & Then A is monotone. Proof. Let D = blockdiag(A11 , . . . , Amm ). By (i), A−1 ii −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 twobytwo block matrices holds if and only if its Schur complement is monotone. Theorem 9. Let A be a twobytwo block matrix A= A11 A12 A21 A22 where T (ai , bi , ci ) stands for a tridiagonal matrix with diagonal coefficients bi and offdiagonal ai and ci . Let % c 1 B = 2 block tridiag T 0, −bi , i , h 2 c T (−ai , 2(ai + bi ), −ai ), T i , −bi , 0 2 & 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 18 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 c≤ ab a+b (29) √ which is stronger than the ellipticity condition c < ( ab). Note also that in the matrix in (25), where a ninepoint 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 timesteps 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 timesteps, we can then use different approximations in space. 4.1 Properties of the solution 4 FINITE DIFFERENCE METHODS FOR PARABOLIC PROBLEMS In this section, we discuss the numerical solution of parabolic problems of the form ∂u = Lu + f (x, t), ∂t x ∈ , t > 0 (30) 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 methods. 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 socalled 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 permafrost.) 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 , 1 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 1 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 t T L1 ϕ1(t ) D ϕ2(t ) E(t) ≤ e−2ct E(0), t >0 1 L0 x Figure 4. Domain of definition for a parabolic problem. Theorem 10. If u is the solution of ∂ 2u ∂u = 2 + f (x, t), ∂t ∂x (x, t) ∈ D (31) 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 twosided 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 (32) 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 2 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 1 e− 4t u0 (y) dy u(x, t) = √ 4πt −∞ (33) where u(x, 0) = u0 (x), u0 ≡ 0 outside a bounded domain . From (33), it readily follows that x2 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 20 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. 4.2 Finite difference schemes: the method of lines 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, dU(t) = Ah U(t) + b(t), dt exp(tA) −−−→ 0 t→∞ Let α(A) = maxi Reλi (A) be the socalled 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 (34) 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 t exp[(t − s)A]b(s) ds, u(t) = exp(tA)u0 + (35) t >0 2 that is 1 E(t) ≤ e 2 β(A)t E(0) 1 1 β(A)−1 e 2 β(A)(t−s) b(s) + 2 ds 0 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. 0 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 2 or 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 timespace domain and beginning in xi . du = Au + b(t), dt 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 ≤e exp(tA) ≥ max eλi t  = max eReλi t = eα(A)t i 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 t≥0 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 socalled θ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)] (37) 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 method V(t + k) = V(t) + k[AV(t) + b(t)] while θ = (1/2) determines the trapezoidal (or Crank– Nicolson) rule k' A [V(t) + V(t + k)] 2 ( + b(t) + b(t + k) V(t + k) = V(t) + 2β(A)t This implies the rightmost inequality of (a). To prove the leftmost inequality, note that for any matrix B, B ≥ maxi λi (B). Hence, i V(t + k) = V(t) + k [AV(t + k) + b(t + k)] (36) that is, by Rayleigh quotient theory, φ (t) ≤ 2β(A)φ(t). Hence, φ(t) ≤ e2β(A)t and for every t, 2 θ = 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 stability. 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 timestepping methods, one can use the Jordan canonical form of the corresponding system matrix. For a general linear system of firstorder difference equations V(t + k) = BV(t) + c(t), t = 0, k, 2k, . . . , V(0) = U0 (38) 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 22 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 = n j =1 cj vj and {λj , vj } are the eigenpairs of A. The corresponding solution of the difference equation is V(t) = n 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, where 1 + θkλj (39) µ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 k<− 2Re(λj ) (2θ − 1)λj 2 j = 1, 2, . . . , where µj = µ(kλj ; θ) and where we have introduced the function , (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 µ(β; θ) = β→−∞ θ , 1−θ −∞, − 0≤θ<1 θ=1 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 timediscretization 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 ) (41) 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 above. 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) = n 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 form [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, . . . (42) We have f = ut − Lu and t+k u(x, t + k) − u(x, t) = 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−θ) t+k ut (x, s) ds−k(1−θ) × f (x, t+k) + kθf (x, t) + 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 t+k ut (x, s) ds τ(x, t; h, k) = t − [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 timestep 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 spacediscretization 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 accurate. + 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 timediscretization error and (2) the spacediscretization 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 t>0 + O(h )( 2 u(4) x + 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 du + F (t, u) = 0, dt t>0 u(4) y ), h, k → 0 (43) 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 (1984). 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 lefthalf plane of the complex plane. Methods with such a property are called Astable methods; see Dahlquist (1963). It is known that multistep timediscretization methods that are Astable can be at most of secondorder of approximation. On the other hand, certain implicit Runge–Kutta methods can have arbitrary high order and still be Astable; 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 (44) 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 2 (45) 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 or 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 onesided 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 24 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 du + F (t, u) = 0, dt 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 /∂u the corresponding parts of the Jacobian matrix ∂ F ), ε → 0. We then multiply the are unbounded as O(ε−1 i i corresponding equation by this parameter to get ε du + F (t, u) = 0, dt 0<t ≤T (46) 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 2 ≤ −ρ(t)(ε(u − v, u − v)), t ≥ t0 so u(t) − v(t) 2 ε ≤ exp t −2ρ(s) ds t0 ≤ u(t0 ) − v(t0 ) 2ε , u(t0 ) − v(t0 ) 2 ε 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 stepbystep 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 problems The traditional numerical integration methods to solve initial value problems du = f (t, u(t)), dt t > 0, u(0) = u0 given (47) are stepbystep methods. Some such familiar methods are Runge–Kutta and linear multistep methods. In stepbystep methods, the error at the current timestep depends on the local error at this step and also on the errors at all previous timesteps. 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 timepoints. In boundary value methods, on the other hand, all approximations at all timepoints 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 timestepping methods but with much larger timesteps than for standard stepbystep methods. A boundary value method can be composed of a sequence of forward stepbystep methods followed by a stabilizing backwardstep method. As a simple example, let u0 be given and un+1 − un−1 − 2kf (tn , un ) = 0, u −u N N−1 n = 1, 2, . . . , N − 1 − kf (tN , u ) = 0 N 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 o [r] 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 (51) o S h = span{φi,j , i = 0, 1, 2, . . . , N − 1, j = 1, 2, . . . , p} Let T dU + F (t, U ), V dt t0 m U, V ∈ H 1 (t0 , T ) a(U ; V ) ≡ ε dt, where H 1 (t0 , T ) is the firstorder Sobolev space of functions with squareintegrable derivatives. To get an approx of the solution of (46), we take a vectorial test imation U [r] function V = ϕi,j , multiply the equation, and after integration, we obtain + ti+1 * d U [r] [r] ; ϕ ) = ), ϕ ε dt = 0, + F (t, U a(U i,j i,j dt ti−1 i = 1, . . . , N − 1 j =0 + ti+1 * d U [r] [r] ; ϕ ) = ), ϕ + F (t, U ε dt = 0, a(U i,j i,j dt ti i = 0, 1, . . . , N − 1 j = 1, 2, . . . , p − 1 and at tN = T , we get + tN * d U ; ϕ[r] ) = ), ϕ[r] dt = 0 + F (t, U ε a(U N,0 N,0 dt tN−1 (48) (49) (50) o 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 ti+1 ti−1 % * ε dU dU − dt dt & ) , V dt = 0 + F (t, U ) − F (t, U + T t0 T t0 U − UI 2 dt ≤ C0 h2(p+1) , ,2 T , dU dUI , 2p , , , dt − dt , dt ≤ C1 h t0 T U t0 2 p+1 dt , , , dU ,2 , , dt (52) , dt , p+1 T p+1 2 = 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 p N−1 Theorem 14. Let U be the solution of (46), and conditions [r] ϕi,j = U (t )φ + U 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 o ) = η − θ, we see that θ = U − U ∈S . UI ) + (UI − U I h 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: 2 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 2 p+2 + U 2 p+1 .(1/2) h→0 , where ν = 1 if p = 1, 1 ≥ ν ≥ (1/2) if p = 3, 5, . . . and ν = 0 if p is even, and V 2 = 1 (εV (T ), V (T )) + 2 T ρ(t) V (t) 2 dt t0 Proof. For a detailed proof, see the supplement of (Axels son and Verwer, 1984). We have assumed that F is Lipschitzcontinuous, 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 twosided 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 26 Finite Difference Methods relations (48), (50) imply ti+1 −U )=2 φ +U φ ε( U F (t, U i+1 i−1 i−1 i−1 i i ti−1 φ )φ dt +U i+1 i+1 i tN −U )= φ ε(U F (t, U N N−1 N−1 N−1 + UN φN )φN dt tN−1 (53) 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, ti h 1 F φi dt = [Fi−1 φi (ti−1 ) + Fi φi (ti )] = hFi 2 2 ti−1 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) h 1 (Fi−1 + Fi ) + t − ti + 2 2 ti−1 ≤ t ≤ ti 1 (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, ti h h F φi dt (Fi−1 + Fi ) + (Fi − Fi−1 ) 4 12 ti−1 = h (F + 2Fi ), 6 i−1 i = 1, 2, . . . , N − 1 and similarly ti ti+1 F φi dt h (F + 2Fi ) 6 i+1 Hence, the generalized midpoint rule (53) takes the form h −U εU = Fi−1 + 4Fi + Fi+1 , i+1 i−1 3 (54) i = 1, 2, . . . , N − 1, h F + 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 Astable. 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 timestep. 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 secondorder 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 timesteps, so explicit methods are generally more costly except where there is a need to choose small timesteps 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 timestep 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 timesteps, 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 timesplitting 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 equation 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 complexvalued functions. We find then, vt