Multigrid

By Ulrich Trottenberg, Cornelius W. Oosterlee and Anton Schuller

Multigrid presents both an elementary introduction to multigrid methods for solving partial differential equations and a contemporary survey of advanced multigrid techniques and real-life applications.Multigrid methods are invaluable to researchers in scientific disciplines including physics, chemistry, meteorology, fluid and continuum mechanics, geology, biology, and all engineering disciplines. They are also becoming increasingly important in economics and financial mathematics.Readers are presented with an invaluable summary covering 25 years of practical experience acquired by the multigrid research group at the Germany National Research Center for Information Technology. The book presents both practical and theoretical points of view.

* Covers the whole field of multigrid methods from its elements up to the most advanced applications* Style is essentially elementary but mathematically rigorous* No other book is so comprehensive and written for both practitioners and students

• Language: English
• Publisher: Academic Press
• Release date: Nov 20, 2000
• ISBN: 9780080479569

Book preview

Table of Contents

Cover image

Title page

Dedication

Copyright

PREFACE

Chapter 1: INTRODUCTION

1.1 TYPES OF PDEs

1.2 GRIDS AND DISCRETIZATION APPROACHES

1.3 SOME NOTATION

1.4 POISSON’S EQUATION AND MODEL PROBLEM 1

1.5 A FIRST GLANCE AT MULTIGRID

1.6 INTERMEZZO: SOME BASIC FACTS AND METHODS

Chapter 2: BASIC MULTIGRID I

2.1 ERROR SMOOTHING PROCEDURES

2.2 INTRODUCING THE TWO-GRID CYCLE

2.3 MULTIGRID COMPONENTS

2.4 THE MULTIGRID CYCLE

2.5 MULTIGRID CONVERGENCE AND EFFICIENCY

2.6 FULL MULTIGRID

2.7 FURTHER REMARKS ON TRANSFER OPERATORS

2.8 FIRST GENERALIZATIONS

2.9 MULTIGRID IN 3D

Chapter 3: ELEMENTARY MULTIGRID THEORY

3.1 SURVEY

3.2 WHY IT IS SUFFICIENT TO DERIVE TWO-GRID CONVERGENCE FACTORS

3.3 HOW TO DERIVE TWO-GRID CONVERGENCE FACTORS BY RIGOROUS FOURIER ANALYSIS

3.4 RANGE OF APPLICABILITY OF THE RIGOROUS FOURIER ANALYSIS, OTHER APPROACHES

Chapter 4: LOCAL FOURIER ANALYSIS

4.1 BACKGROUND

4.2 TERMINOLOGY

4.3 SMOOTHING ANALYSIS I

4.4 TWO-GRID ANALYSIS

4.5 SMOOTHING ANALYSIS II

4.6 SOME RESULTS, REMARKS AND EXTENSIONS

4.7 h-ELLIPTICITY

Chapter 5: BASIC MULTIGRID II

5.1 ANISOTROPIC EQUATIONS IN 2D

5.2 ANISOTROPIC EQUATIONS IN 3D

5.3 NONLINEAR PROBLEMS, THE FULL APPROXIMATION SCHEME

5.4 HIGHER ORDER DISCRETIZATIONS

5.5 DOMAINS WITH GEOMETRIC SINGULARITIES

5.6 BOUNDARY CONDITIONS AND SINGULAR SYSTEMS

5.7 FINITE VOLUME DISCRETIZATION AND CURVILINEAR GRIDS

5.8 GENERAL GRID STRUCTURES

Chapter 6: PARALLEL MULTIGRID IN PRACTICE

6.1 PARALLELISM OF MULTIGRID COMPONENTS

6.1.1 Parallel Components for Poisson’s Equation

6.2 GRID PARTITIONING

6.3 GRID PARTITIONING AND MULTIGRID

6.4 PARALLEL LINE SMOOTHERS

6.5 MODIFICATIONS OF MULTIGRID AND RELATED APPROACHES

Chapter 7: MORE ADVANCED MULTIGRID

7.1 THE CONVECTION-DIFFUSION EQUATION: DISCRETIZATION I

7.2 THE CONVECTION-DIFFUSION EQUATION: MULTIGRID I

7.3 THE CONVECTION-DIFFUSION EQUATION: DISCRETIZATION II

7.4 THE CONVECTION–DIFFUSION EQUATION: MULTIGRID II

7.5 ILU SMOOTHING METHODS

7.6 PROBLEMS WITH MIXED DERIVATIVES

7.7 PROBLEMS WITH JUMPING COEFFICIENTS AND GALERKIN COARSE GRID OPERATORS

7.8 MULTIGRID AS A PRECONDITIONER (ACCELERATION OF MULTIGRID BY ITERANT RECOMBINATION)

Chapter 8: MULTIGRID FOR SYSTEMS OF EQUATIONS

8.1 NOTATION AND INTRODUCTORY REMARKS

8.2 MULTIGRID COMPONENTS

8.3 LFA FOR SYSTEMS OF PDEs

8.4 THE BIHARMONIC SYSTEM

8.5 A LINEAR SHELL PROBLEM

8.6 INTRODUCTION TO INCOMPRESSIBLE NAVIER–STOKES EQUATIONS

8.7 INCOMPRESSIBLE NAVIER-STOKES EQUATIONS: STAGGERED DISCRETIZATIONS

8.8 INCOMPRESSIBLE NAVIER-STOKES EQUATIONS: NONSTAGGERED DISCRETIZATIONS

8.9 COMPRESSIBLE EULER EQUATIONS

Chapter 9: ADAPTIVE MULTIGRID

9.1 A SIMPLE EXAMPLE AND SOME NOTATION

9.1.2 Hierarchy of Grids

9.2 THE IDEA OF ADAPTIVE MULTIGRID

9.3 ADAPTIVE MULTIGRID AND THE COMPOSITE GRID

9.4 REFINEMENT CRITERIA AND OPTIMAL GRIDS

9.5 PARALLEL ADAPTIVE MULTIGRID

9.6 SOME PRACTICAL RESULTS

Chapter 10: SOME MORE MULTIGRID APPLICATIONS

10.1 MULTIGRID FOR POISSON-TYPE EQUATIONS ON THE SURFACE OF THE SPHERE

10.2 MULTIGRID AND CONTINUATION METHODS

10.3 GENERATION OF BOUNDARY FITTED GRIDS

10.4 L, SS: A GENERIC MULTIGRID SOFTWARE PACKAGE

10.5 MULTIGRID IN THE AERODYNAMIC INDUSTRY

10.6 HOW TO CONTINUE WITH MULTIGRID

Appendixes

Appendix A: AN INTRODUCTION TO ALGEBRAIC MULTIGRID

Appendix B: SUBSPACE CORRECTION METHODS AND MULTIGRID THEORY

Appendix C: RECENT DEVELOPMENTS IN MULTIGRID EFFICIENCY IN COMPUTATIONAL FLUID DYNAMICS

REFERENCES

INDEX

Dedication

Dedicated to Linde, Lukas, Sophia, Kristian, Eva, Filipp, Katharina, Anasja, Wim, Agnes, Annette, Sonja and Timo

Copyright

This book is printed on acid-free paper.

Copyright © 2001 by ACADEMIC PRESS

All Rights Reserved.

No part of this publication may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without the prior permission in writing from the publisher.

Academic Press

A Harcourt Science and Technology Company Harcourt Place, 32 Jamestown Road, London NW1 7BY, UK http://www.academicpress.com

Academic Press

A Harcourt Science and Technology Company 525 B Street, Suite 1900, San Diego, California 92101-4495, USA http://www.academicpress.com

ISBN 0-12-701070-X

Library of Congress Catalog Number: 00-103940

A catalogue record of this book is available from the British Library

Typeset by Newgen Imaging Systems (P) Ltd., Chennai, India Printed and bound in Great Britain by MPG Books, Cornwall

01 02 03 04 05 MP 8 7 6 5 4 3 2 1

PREFACE

This book is intended for multigrid practitioners and multigrid students: for engineers, mathematicians, physicists, chemists etc.

We will give a systematic introduction to basic and advanced multigrid. Our intention is to give sufficient detail on multigrid for model problems so that the reader can solve further applications and new problems with appropriate multigrid techniques. In addition, we will cover advanced techniques (adaptive and parallel multigrid) and we will present applications based on these new developments.

Clearly, we would not have written this book if we had thought that there was a better book that fulfilled the same purpose. No doubt, there are a number of excellent multigrid books available. However, in our opinion all the other books emphasize multigrid aspects different from those that we are interested in and which we want to emphasize.

Mathematical results that can be rigorously proved may be formulated in different ways. Practically oriented mathematicians often prefer a presentation in which the assumptions and results are motivated by some typical examples and applications. These assumptions are usually not the most general ones, and thus the results may not be the strongest that can be obtained. We prefer such a presentation of results, as we want to provide as much motivation and direct understanding as possible. However, in many cases, we add some remarks about generalizations and extensions, and about weaker assumptions.

With respect to multigrid theory, we give (elementary) proofs, wherever we regard them as helpful for multigrid practice. All the theoretical tools which we use should be understood by mathematicians, scientists and engineers who have a general background in analysis, linear algebra, numerical analysis and partial differential equations (PDEs). Wherever more sophisticated mathematical tools are needed to derive practically relevant theoretical results, we cite the appropriate literature. However, we are not interested in theory as an end in itself.

This book has three authors and, in addition, three guest authors. While the guest contributions are supposed to be fully consistent with the contents of the book and fit with its general philosophy and message, they are still independent and self-contained. The guest authors are experts in the fields they present here and they use their own style to express their views and approaches to multigrid.

The three main authors, on the other hand, have agreed on all the material that is presented in this book. They did not distribute the chapters among themselves and did not distribute the responsibility. They also agreed on the way the material is presented.

Multigrid methods are generally accepted as being the fastest numerical methods for the solution of elliptic partial differential equations. Furthermore, they are regarded as among the fastest methods for many other problems, like other types of partial differential equations, integral equations etc. If the multigrid idea is generalized to structures other than grids, one obtains multilevel, multiscale or multiresolution methods, which can also be used successfully for very different types of problems, e.g. problems which are characterized by matrix structures, particle structures, lattice structures etc. However, the literature does not have a uniform definition of the terms multigrid, multilevel etc.

This book is devoted to PDEs and to the "algebraic multigrid approach" for matrix problems.

We assume that the reader has some basic knowledge of numerical methods for PDEs. This includes fundamental discretization approaches and solution methods for linear and nonlinear systems of algebraic equations. Implicitly, this also means that the reader is familiar with basics of PDEs (type classification, characteristics, separation of variables etc. see, for example, [395]) and of direct and iterative solvers for linear systems.

We will not, however, assume detailed knowledge about existence theories for PDEs, Sobolev spaces and the like. In this respect, the book is addressed to students and practitioners from different disciplines. On the other hand, in some sections, advanced applications are treated, in particular from computational fluid dynamics. For a full understanding of these applications, a basic knowledge of general PDEs may not be sufficient. In this respect, we will assume additional knowledge in these sections and we will give references to background material in the corresponding fields.

We do not assume that the reader works linearly with this book from the beginning to the end though this is suitable to obtain a good insight into multigrid and its relation to similar approaches. The multigrid beginner may well skip certain sections. We will lead the reader in this respect through the book, pointing out what can be skipped and what is needed.

The overall structure of the book is determined by its chapters. The first half of the book (Chapters 1–6) discusses standard multigrid, the second half (Chapters 7–10) deals with advanced approaches up to real-life applications. Accordingly, the style and presentation in the first half is more detailed. In addition to the basic techniques introduced in the first six chapters, we add many more specific remarks and algorithmical details. These may not be so interesting for beginners but should be helpful for practitioners who want to write efficient multigrid programs. Mistakes that are easily made are mentioned in several places.

The second part of the book (Chapters 7–10) is presented in a more condensed form, i.e. in a more research oriented style.

This structure of the book is also reflected by the nature of the equations and applications we deal with. There is no doubt about the fact that multigrid methods work excellently for nicely elliptic PDEs. This is confirmed by rigorous mathematical theory.

For typical real-life applications (PDE systems with nonelliptic features and nonlinear terms), however, such a theory is generally not available. Nevertheless, as we will see in this book, multigrid can be applied to such problems although they may not be nicely elliptic or even not elliptic at all. In answering the question when does multigrid work?, we will give insight, based on 20 years of multigrid practice and multigrid software development.

ACKNOWLEDGMENTS

We would like to thank the many friends, colleagues, collaborators and students who have helped us to write this book. First, we thank the three guest authors, Achi Brandt, Peter Oswald and Klaus Stüben for their contributions. Klaus Stüben was also the chief reader of the manuscript. His criticism and comments were extremely helpful.

Rudolph Lorentz also commented extensively on the manuscript and checked our English.

Achi Brandt, who, to our regret, never wrote a multigrid book himself, closely followed our progress and made many helpful comments.

Others who commented on our manuscript included: Tanja Füllenbach, Jürgen Jakumeit, Rene Redler, Roman Wienands and Kristian Witsch. We would like to thank them for their efforts and suggestions.

All colleagues working in the GMD–Institute for Algorithms and Scientific Computing (SCAI) who supported us are gratefully acknowledged. In particular, we thank Ingrid Filter, Wolfgang Joppich, Axel Klawonn, Johannes Linden, Hans-Günther Reschke, Hubert Ritzdorf, Horst Schwichtenberg, Frauke Sprengel, Barbara Steckel, Clemens-August Thole and Monika Wendel.

Some of the colleagues and friends whom we thankfor clarifying discussions, interesting cooperations and for pointers to the literature include: Wolfgang Dahmen, Francisco Gaspar, Michael Griebel, Norbert Kroll, Oliver McBryan, Takumi Washio, Pieter Wesseling, Olof Widlund, Jinchao Xu, Irad Yavneh and Christoph Zenger.

We would like to thank Eric Peters for the cover design.

The initiative to write this book was strongly influenced by Olof Widlund, who invited Ulrich Trottenberg to the Courant Institute of Mathematical Sciences of the New York University, New York City as a visiting professor to teach a multigrid course for graduate students.

Of course, this book could not have been written without the support and patience of our families, to whom we owe most.

INTRODUCTION

We start this chapter with a short introduction of some of the equations that we will treat in this book in Section 1.1 and with some information on grids and discretization approaches in Section 1.2. In Section 1.3 we will introduce some of our terminology. The 2D Poisson equation with Dirichlet boundary conditions is the prototype of an elliptic boundary value problem. It is introduced and discussed in Section 1.4. In Section 1.5 we will take a first glance at multigrid and obtain an impression of the multigrid idea. Some facts and methods on basic numerics are listed in Section 1.6.

1.1 TYPES OF PDEs

As we will see in this book, elliptic boundary value problems are the type of problem to which multigrid methods can be applied very efficiently. However, multigrid or multigrid-like methods have also been developed for many PDEs with nonelliptic features.

We will start with the usual classification of second-order scalar 2D PDEs. Generalizations of this classification to 3D, higher order equations or systems of PDEs can be found [150]. We consider equations Lu = f where

     (1.1.1)

with coefficients aij, ai, ao and a right-hand side f which, in general, may depend on x, y, u, ux, uy (the quasilinear case). In most parts of this book, Lu = f is assumed to be a linear differential equation, which means that the coefficients and the right-hand side f only depend on (x, y). L is called

.

In general, this classification depends on (x, y) and, in the nonlinear case, also on the solution u. Prototypes of the above equation are

• the Poisson equation −Δu = −uxx − uyy = f,

• the wave equation uxx − uyy = 0,

• the heat equation uxx − uy = 0.

Since multigrid methods work excellently for nicely elliptic problems, most of our presentation in the first chapters is oriented to Poisson’s and Poisson-like equations. Other important model equations that we will treat in this book include

• the anisotropic model equation −εuxx − uyy = f,

• the convection-diffusion equation −εΔu + a1ux + a2uy = f,

• the equation with mixed derivatives −Δu + τuxy = f.

All these equations will serve as model equations for special features and complications and are thus representative of a larger class of problems with similar features. These model equations depend crucially on a parameter ε or τ. For certain parameter values we have a singular perturbation: the type of the equation changes and the solution behaves qualitatively different (if it exists at all). For instance, the anisotropic equation becomes parabolic for ε → 0, the equation with mixed derivatives is elliptic for |τ| < 2, parabolic for |τ| = 2 and hyperbolic for |τ| > 2. All the model equations represent classes of problems which are of practical relevance.

In this book, the applicability of multigrid is connected to a quantity, the "h-ellipticity measure Eh", that we will introduce in Section 4.7. This h-ellipticity measure is not applied to the differential operator itself, but to the corresponding discrete operator. It can be used to analyze whether or not the discretization is appropriate for a straightforward multigrid treatment. Nonelliptic problems can also have some h-ellipticity if discretized accordingly.

The above model equations, except the wave and the heat conduction equations, are typically connected with pure boundary conditions. The wave and the heat conduction equations are typically characterized by initial conditions with respect to one of the variables (y) and by boundary conditions with respect to the other (x).

We will call problems which are characterized by pure boundary conditions space-type problems. For problems with initial conditions, we interpret the variable for which the initial condition is stated as the time variable t (= y), and call these problems time-type. Usually these problems exhibit a marching or evolution behavior with respect to t. Space-type equations, on the other hand, usually describe stationary situations. Note that this distinction is different from the standard classification of elliptic, hyperbolic and parabolic. Elliptic problems are usually space-type while hyperbolic and parabolic problems are often time-type. However, the stationary supersonic full potential equation, the convection equation (see Chapter 7) and the stationary Euler equations (see Section 8.9) are examples of hyperbolic equations with space-type behavior. (In certain situations, the same equation can be interpreted as space- or as time-type, and each interpretation may have its specific meaning. An example is the convection equation.)

In this book, we will present and discuss multigrid methods for space-type problems.

For time-type problems, there is usually the option to discretize the time variable explicitly, implicitly or in a hybrid manner (i.e. semi-implicitly). The implicit or semi-implicit discretization yields discrete space-type problems which have to be solved in each time step. We will briefly consider multigrid approaches for the solution of such space-type problems which arise in each time step in Section 2.8. Typically, these problems have similar but more convenient numerical properties than the corresponding stationary problems with respect to multigrid.

Remark 1.1.1

Some authors have proposed multigrid approaches which include the time direction directly. These approaches consider time to be just another space-type direction. Such approaches are discussed in [78, 175, 198, 199, 396] for so-called parabolic multigrid and multigrid-like methods based on waveform relaxation.

1.2 GRIDS AND DISCRETIZATION APPROACHES

In this book, we assume that the differential equations to be solved are discretized on a suitable grid (or, synonymously, mesh). Here, we give a rough survey (with some examples) of those types of grids which are treated systematically in this book and those which are only touched on. We also make some remarks about discretization approaches.

1.2.1 Grids

The general remarks in this section may not be so interesting to multigrid beginners. They may start with Section 1.4 on Poisson’s equation and return to this section later.

Most parts of this book refer to: Cartesian grids, boundary-fitted logically rectangular grids and block-structured boundary-fitted grids. Figure 1.1 is an example of a Cartesian grid. For simple domains with simple boundaries, Cartesian grids are numerically convenient. We will use them for Model Problem 1 (see Section 1.3.2) and several other cases.

Figure 1.1 A square Cartesian grid (left) in a domain Ω and a boundary-fitted curvilinear (logically rectangular) grid (right).

Figure 1.1 also gives an example of a boundary-fitted grid. Boundary-fitted grids will be used in more advanced examples in this book. A systematic introduction into boundary-fitted grids is given in [391]. We will discuss the generation of boundary-fitted grids in Section 10.3.

In the context of boundary-fitted grids, there are two different approaches. In the first, coordinate transformations are used to obtain simple (for example rectangular) domains, and correspondingly simple (rectangular) grids. Here the differential (and/or the discrete) equations are transformed to the new curvilinear coordinates. In the second approach, the computations are performed in the physical domain with the original (nontransformed) equations. In this book we concentrate on the second approach.

Block-structured boundary-fitted grids are used if the given domain cannot (or cannot reasonably) be mapped to a rectangular domain, but can be decomposed into a finite number of subdomains each of which can be covered with a boundary-fitted grid. Quite complicated domains may be treated with this approach, as can be seen in Fig. 1.2. Block-structured boundary-fitted grids are discussed in Chapters 6 and chapter 8–10.

Figure 1.2 Boundary-fitted block-structured grid around a car.

In many software packages, unstructured, irregular grids are used today. These grids have become quite popular because unstructured automatic mesh generation is much easier than the generation of block-structured grids for very complicated 2D and 3D domains.

An example of an unstructured grid is the grid around a car during a crash simulation. A part of this grid is presented in Fig. 1.3.

Figure 1.3 Unstructured grid around a car in a crash simulation application.

From the multigrid point of view, unstructured grids are a complication. For a given unstructured grid, it is usually not difficult to define a sequence of finer grids, but it may be difficult to define a sequence of reasonable coarser grids. (A hierarchy of coarse grids is needed for multigrid.)

The algebraic multigrid (AMG) method presented in Appendix A constructs a hierarchy of coarse grids automatically and is thus particularly well-suited for problems on unstructured grids.

Although unstructured grids are widely used, even complicated domains allow other than purely unstructured grid approaches. Often a hybrid approach, an unstructured grid close to the boundary and a structured grid in the interior part of the domain (or vice versa) is suitable for the treatment of such complicated domains.

More generally, all types of grids mentioned above can be used in the context of overlapping grids. A typical situation for the use of overlapping grids is when an overall Cartesian grid is combined with a local boundary-fitted grid (Chimera technique). An example of this approach is shown in Fig. 1.4. Such overlapping grids are also called composite grids.

Figure 1.4 An overlapping grid around an object.

Finally, we mention self-adaptive grids. They are constructed automatically during the solution process according to an error estimator that takes the behavior of the solution into account. Self-adaptive grids are very natural in the multigrid context. We will treat the self-adaptive multigrid approach systematically in Chapter 9. An example is given in Fig. 1.5.

Figure 1.5 Adaptively refined grid around an airfoil.

1.2.2 Discretization Approaches

In principle, any type of grid can be used with any type of discretization approach. In practice, however, finite difference and finite volume methods are traditionally used in the context of Cartesian, logically rectangular and block-structured boundary-fitted grids, whereas finite elements and finite volumes are widely used in the context of unstructured grids.

In this book, we will focus on finite difference and finite volume discretizations on structured and block-structured grids.

Finally, another important choice to be made in discretization is the arrangement of the unknowns within a grid. The unknowns can be defined at the vertices of a grid (vertex-centered location of unknowns). Another option is the choice of the unknowns at cell centers (cell-centered location of unknowns).

For systems of PDEs it is possible to choose different locations for different types of unknowns. A well-known example is the staggered grid for the system of the incompressible Navier–Stokes equations discussed in Chapter 8, in which pressure unknowns are placed at the cell centers and velocity components at cell boundaries. Examples of a vertex-centered, a cell-centered and a staggered Cartesian grid are sketched in Fig. 1.6. Often, the discrete equations are defined at the same locations as the unknowns. It is hard to say which location of unknowns and which location of the discrete equations is best in general. Often these choices depend on the type of boundary conditions and on the application. In the following chapters we mainly present results for vertex-centered locations of unknowns. Multigrid components for cell-centered arrangements of unknowns are presented in Section 2.8.4.

Figure 1.6 Three arrangements of unknowns in a Cartesian grid: (a) a vertex-centered grid; (b) a cell-centered grid; (c) a staggered grid.

1.3 SOME NOTATION

We start with some general notation. In this book, we use the classical formulation of differential equations with differential (and boundary) operators rather than a weak formulation. For discretization, we use the terminology of discrete differential operators. This also means that we use grid functions rather than vectors and grid operators rather than matrices. In our opinion, the grid-oriented notation emphasizes more clearly the correspondence between the discrete and continuous problems, and between the discrete formulations on different grids of the multigrid structure. In that respect, we regard it as a natural formulation for the multigrid approach. If, for example, the discrete differential operator can be described by a difference stencil (see Section 1.3.4) this is clearly a more condensed formulation than a matrix notation. On the other hand, there are situations in which the matrix notation has some advantages and is more general. Then we will not hesitate to use it.

1.3.1 Continuous Boundary Value Problems

We denote scalar linear boundary value problems by

     (1.3.1)

Here x = (x1,…, xd)T is a given open bounded domain with boundary Γ. LΩ is a linear (elliptic) differential operator on Ω and Lγ represents one or several linear boundary operators. fω denotes a given function on ω and fγ one or several functions on γ. We always denote solutions of (1.3.1) by u = u(x). For brevity, we also simply write Lu = f instead of (1.3.1). Most concrete considerations refer to the cases d = 2 and d = 3. (Multigrid is usually not needed for d = 1, in which case it degenerates to direct or other simple well-known solvers, see Section 6.4.1.)

In the case d = 2 or d = 3, we will usually write (x, y) instead of (x1, x2) and (x, y, z) instead of (x1, x2, x3).

This and the following chapters essentially refer to scalar equations. Chapter 8 and other parts of the book, however, will refer to systems of PDEs. The above notation is also used in that case. LΩ then stands for a vector of differential operators and u, f etc. are vector-valued functions.

1.3.2 Discrete Boundary Value Problems

The following considerations are formulated for the 2D case. They can, of course, be directly generalized to higher dimensions. The discrete analog of (1.3.1) is denoted by

     (1.3.2)

Here h is a (formal) discretization parameter. Using the infinite grid

     (1.3.3)

where h = (hx, hy) is a vector of fixed mesh sizes, we define Ωh = Ω ∩ Gh and Γh as the set of discrete intersection points of the grid lines with Γ. In the special case of square grids, we simply identify h = hx = hy.

The discrete solution uh is a function defined on the grid Ωh ∪ Γhare discrete analogs of fΩ and fΓ, respectively, where fΩ, fΓ are restricted to the grid. Instead of uh(x, y)=uh(xi, yj)=uh(ihx, jhy), we sometimes simply write ui, j.

is also called a discrete (differential) or difference operatora discrete boundary operator.)

Clearly the concrete definitions of Ωh, Γh etc. also depend on the given PDE, the domain Ω, the boundary conditions, the grid approach and the discretization.

For ease of presentation, we will first assume that discrete boundary equations are eliminated from (1.3.2). (In general, a proper multigrid treatment of boundary conditions needs an explicit consideration of noneliminated boundary conditions. This is discussed in more detail in Section 5.6.)

eliminated boundary conditions can be considered as a natural approach. We then simply write

     (1.3.4)

Here uh and fh are grid functions on Ωh and Lh is a linear operator

     (1.3.5)

denotes the linear space of grid functions on Ωh. Clearly, (1.3.4) can be represented as a system of linear algebraic equations. However, we usually consider it as one grid equation.

Even if the boundary conditions are not eliminated, we may use the notation (1.3.4) for the discrete problem. The notation then represents an abstract equation (in a finite dimensional space).

1.3.3 Inner Products and Norms

. Most of our considerations (and many of our measurements) will be based on the Euclidean inner product

     (1.3.6)

where #Ωh is the number of grid points of Ωh. The scaling factor (#ΩhThe corresponding operator norm for discrete operators Lh is the spectral norm

where Bh , and where ρ is the spectral radius.

For Lh, symmetric and positive definite, we also consider the energy inner product

     (1.3.7)

, which is given by

     (1.3.8)

For practical purposes the infinity norm

     (1.3.9)

are commonly used.

1.3.4 Stencil Notation

(on a Cartesian or a logically rectangular grid) the stencil terminology is convenient. We restrict ourselves to the case d = 2. The extension to d = 3 (and to general d) is straightforward (see Section 2.9). We first introduce the stencil notation for the infinite grid Gh (defined in Section 1.3.2). On Gh, we consider grid functions

defines an operator on the set of grid functions by

     (1.3.10)

Here we assume that only a finite number of coefficients Sk1k2 are nonzero.

Many of the stencils we will consider will be five-point or compact nine-point stencils

     (1.3.11)

are usually given only on a finite grid ωh, we usually have to restrict the stencil to ωh (instead of Gh). Near boundary points the stencils may have to be modified. In square or rectangular domains ω, which are the basis for the examples in to ωh will depend on (x, y):

if the coefficients of LΩ and/or LhΩ depend on (x, y).

1.4 POISSON’S EQUATION AND MODEL PROBLEM 1

In this section, we introduce Model Problem 1, the classical model for a discrete elliptic boundary value problem. Every numerical solver has been applied to this problem for comparison.

Model Problem 1

At many places in this book, we will study in detail the discrete Poisson equation with Dirichlet boundary conditions

     (1.4.1)

in the unit square is the standard five-point O(h²) approximation (explained below) of the partial differential operator L,

     (1.4.2)

on the square grid Gh.

O(h²) here means that one can derive consistency relations of the form

, for example).

To illustrate exactly what the elimination of boundary conditions means, we consider the following example.

Example 1.4.1

The discrete Poisson equation with eliminated Dirichlet boundary conditions can formally be written in the form (1.3.4). For (x, y) ∈ Ωh not adjacent to a boundary this means:

The notation

is a first example of the stencil notation (1.3.11).

For (x, y) ∈ Ωh adjacent to a (here: west) boundary, Lh (1.3.4) reads

For (x, y) ∈ ωh in a (here: the north-west) corner we have

In this example, elimination of the boundary conditions is simple. Often, elimination of boundary conditions may be complicated and is not preferred. In such cases, a particular treatment of the boundary conditions in the multigrid process may be needed (see Section 5.6).

1.4.1 Matrix Terminology

Discrete operators Lh are often represented by matrices Ah. Each matrix row then represents connections of one unknown in the discretization of a PDE to its neighbor unknowns. Which of the matrix entries is different from 0, depends on the ordering of grid points, i.e., the ordering of the components of the vector of unknowns.

As an example we consider Model Problem 1. For a column- or row-wise ordering of grid points (also called lexicographical ordering, see Fig. 1.7(a), starting with points at the left lower corner) and eliminated Dirichlet boundary conditions, the resulting matrix is a block tridiagonal matrix with a regular sparsity pattern:

Figure 1.7 (a) Lexicographic; (b) red–black ordering of grid points.

     (1.4.3)

where

Due to the elimination of Dirichlet boundary points, every matrix row corresponding to a grid point near a boundary has only three or two entries connecting them to neighbor grid points. This can be seen, for example, in the first rows of the matrix, where only two or three –1 entries can be found instead of four for interior points.

The dependence of the matrix structure on the ordering of the grid points can be seen when we write the matrix down for a red–black ordering of the grid points (see Fig. 1.7(b)). First all unknowns at odd (red) grid points are considered, then the unknowns at the even (black) points. The corresponding matrix Ah is now a block matrix with blocks Arr, representing the connections of the red grid points to red grid points, Arb the connections of the red points to the black points, Abr the connections of black points to red points, and Abb the connections of black points to black points. So:

     (1.4.4)

For Model Problem 1, the blocks Arr and Abb are diagonal matrices with 4/h² of the above example in Fig. 1.7(b) is

     (1.4.5)

Remark 1.4.1

As we will see below (Section 1.5), in a multigrid algorithm it is usually not necessary to build up the matrix Ah coming from the discretization. The multigrid components are based on local operations; multiplications and additions are carried out grid point by grid point. The storage that is needed in a multigrid code mainly consists of solution vectors, defects and right-hand sides on all grid levels.

1.4.2 Poisson Solvers

Table 1.1 gives an overview on the complexity of different solution methods (including fast Poisson solvers) applied to Model Problem 1. Here direct and iterative solvers are listed. For the iterative solvers, we assume an accuracy (stopping criterion) in the range of the discretization accuracy. This is reflected by the log ε term. The full multigrid (FMG) variant of multigrid which we will introduce in Section 2.6 is a solver up to discretization accuracy.

Table 1.1 Complexity of different solvers for the 2D Poisson problem (N denotes the total number of unknowns).

It is generally expected that the more general a solution method is, the less efficient it is and vice versa. Multigrid is, however, a counter example for this pattern–indeed multigrid will turn out to be a general principle. On the other hand, Table 1.1 shows that the most efficient version of multigrid yields an algorithm that is at least as efficient as the highly efficient and highly specialized fast Poisson solvers. Here, the total reduction method and the Buneman algorithm are typical fast Poisson solvers, very efficient but essentially designed and tailored exclusively for elliptic model problems.

1.5 A FIRST GLANCE AT MULTIGRID

1.5.1 The Two Ingredients of Multigrid

In this section, we introduce the multigrid idea heuristically for Model Problem 1. We use the grid function oriented notation introduced in Section 1.3.2.

The iteration formula of the classical lexicographical Gauss–Seidel method (GS-LEX) for Poisson’s equation reads

     (1.5.1)

where (xi, yj) ∈ Ωhbefore and after an iteration, respectively.

If we apply (1.5.1) to Poisson’s equation, we recognize the following phenomenon. After a few iteration steps, the error of the approximation becomes smooth. It doesn’t necessarily become small, but it does become smooth. See Fig. 1.8 for an illustration of this error smoothing effect. Looking at the error

Figure 1.8 Influence of lexicographic Gauss–Seidel iteration (1.5.1) on the error (Model Problem 1).

Formula (1.5.1) means

     (1.5.2)

Obviously, the iteration formula can be interpreted as an error averaging process. Error smoothing is one of the two basic principles of the multigrid approach.

Smoothing principle

Many classical iterative methods (Gauss–Seidel etc.) if appropriately applied to discrete elliptic problems have a strong smoothing effect on the error of any approximation.

The other basic principle is the following: a quantity that is smooth on a certain grid can, without any essential loss of information, also be approximated on a coarser grid, say a grid with double the mesh size. In other words: if we are sure that the error of our approximation has become smooth after some iteration steps, we may approximate this error by a suitable procedure on a (much) coarser grid.

Qualitatively, this is the coarse grid approximation principle.

Coarse grid principle

A smooth error term is well approximated on a coarse grid. A coarse grid procedure is substantially less expensive (substantially fewer grid points) than a fine grid procedure.

As this principle holds for error or correction quantities, we also speak of the coarse grid correction (CGC) principle.

considered as a function of the discrete variables x and y can be written as

     (1.5.3)

For (x, y) ∈ Ωh, the functions

     (1.5.4)

are the discrete eigenfunctions of the discrete operator Δh. The fact that this error becomes smooth after some iteration steps means that the high frequency components in (1.5.3), i.e.

     (1.5.5)

become small after a few iterations whereas the low frequency components, i.e.

     (1.5.6)

hardly change. The distinction between high and low frequencies is important in the multigrid context. In the following subsection, we will give a first definition of high and low frequencies and show how this concept is related to the coarse grid under consideration.

1.5.2 High and Low Frequencies, and Coarse Meshes

We again consider Model Problem 1 on a grid Ωh with mesh size h = 1/n. Additionally, we consider Model Problem 1 on a coarser grid ΩH with mesh size H > h. Assuming that n is an even number, we may, for instance, choose H = 2h, which is a very natural choice in the multigrid context. This choice of the coarse grid is, therefore, called standard coarsening.

in (1.5.4). For given (k, l), we consider the (four) eigenfunctions

and observe that they coincide on Ω2h in the following sense:

This means that these four eigenfunctions cannot be distinguished on Ω2h. (For k or l = n/2, the ϕk, l vanish on Ω2h.) This gives rise to the following definition of low and high frequencies:

Definition

(in the context of Model Problem 1) For k, l ∈ {1,…,n – 1}, we denote ϕk, l to be an eigenfunction (or a component) of

low frequency if max (k, l) < n/2,

high frequency if n/2 ≤ max (k, l) < n.

Obviously, only the low frequencies are visible on ω2h since all high frequencies coincide with a low frequency on ω2h (or vanish on ω2h). The fact that high frequencies coincide with low ones is also called aliasing of frequencies. For the 1D case with n = 8, we illustrate the above definition in Fig. 1.9. We summarize this consideration:

Figure 1.9 Low and high frequency components for a ID example (n = 8).

The low frequency components also represent meaningful grid functions on a grid ω2h with double the mesh size, whereas the high frequency components do not. Their high frequencies are not visible on the ω2h grid.

If we apply this distinction of low and high frequencies to the representation of vh (x, y) in (1.5.3), we can decompose the sum in (1.5.3) into corresponding partial sums:

where

and

Remark 1.5.1 (H = 4h and other choices of H)

From our definition, it immediately becomes clear that the terms high frequency and low frequency are related to both the fine grid ωh and the coarse grid ωH that we consider. (If we want to emphasize this dependence on the grids ωh and ωH, we will speak of (h, H)-low and (h, H)-high frequencies.) If, for example, we would choose

(assuming that n is a multiple of 4), our definition of high and low frequencies would have to be modified in the following way:

The choice H = 4h is not very practical in the multigrid context since it usually does not lead to the most efficient algorithms, but it is not out of the question, either.

Other, more practical choices of coarse grids, different from standard coarsening, are introduced in Section 2.3.1.

1.5.3 From Two Grids to Multigrid

In the following, we will continue our heuristic considerations a little further extending them from two-grid levels Ωh, Ω2h to a sequence of levels. In this setting, we will also be able to explain the multigrid (not only the two-grid) idea. Fig. 1.10 shows such a sequence of grids.

Figure 1.10 A sequence of coarse grids starting with h = 1/16.

We assume additionally that n is a power of 2 meaning that h = 2−p. Then we can form the grid sequence

     (1.5.7)

just by doubling the mesh size successively. We assume that this sequence ends with a coarsest grid Ωh0 (which may well be the grid consisting of only one interior grid point (x, y) = (1/2, 1/2), i.e. h0 = 1/2, see Fig. 1.10).

We now consider a decomposition of the error representation (1.5.3) into partial sums which correspond to the grid sequence (1.5.7), having the following idea in mind. In the same way as we have distinguished low and high frequencies with respect to the pair of grids Ωh and Ω2h in the previous section, we now make an additional distinction between high and low frequencies with respect to the pair Ω2h and Ω4h. We continue with the pair Ω4h and Ω8h and so on.

As discussed above, by using Gauss-Seidel type iterations on the original Ωh grid, we cause the (h, 2h)-high frequency error components to become small rapidly. The remaining low frequency error components are visible and can thus be approximated on Ω2h. (Of course, an equation which determines the low frequency error has to be defined in a suitable way on Ω2h.) Performing Gauss-Seidel-type iteration not only on the original Ωh grid, but also on Ω2h and correspondingly on Ω4h etc., causes the (2h, 4h)-high frequency, the (4h, 8h)-high frequency components etc. to decrease rapidly. Only on the coarsest grid might it be necessary to do something special, which is trivial if Ωh0 consists of only one (or few) point(s). Altogether, this leads to a very fast overall reduction of the error.

Summary

What we have described and explained so far is the basic idea of multigrid. Our description is not at all an algorithm. It is not even a clear verbal description of an algorithmic principle. For example, we have neither said precisely what a Gauss-Seidel iteration on Ω2h, Ω4h etc. means algorithmically and to which grid functions it is to be applied, nor how to go from one level to the other etc. However, the flavor of multigrid is in it.

Flavor of multigrid

Gauss–Seidel iteration (or more generally, an appropriate iterative scheme) on different grid levels gives rapid reduction of the corresponding high frequency components and as this process covers all frequencies, a rapid reduction of the overall error can be achieved.

1.5.4 Multigrid Features

The multigrid idea is very fundamental and can also be applied in other contexts. Accordingly, there are also several different ways to view multigrid. In order to clarify our understanding of multigrid, we want to briefly discuss different multigrid aspects and to point out those multigrid features which we regard as the most significant ones.

Multigrid as an iterative linear solver

The most basic way to view multigrid is to consider it as an iterative linear solver for a discrete elliptic boundary value problem. Here we assume the problem, the grid and the discretization to be given and fixed. A characteristic feature of the iterative multigrid approach is that the multigrid convergence speed is independent of the discretization mesh size h and that the number of arithmetic operations per iteration step is proportional to the number of grid points. The multigrid principle allows us to construct very efficient linear solvers, and, in that respect, the iterative approach is important and fundamental. This approach is described in detail in Chapter 2 of this book. On the other hand, this view is somewhat restricted and does not exploit the full potential of the multigrid idea. For example, even the direct application of multigrid to nonlinear problems (discussed in detail in Section 5.3) is not covered by this approach.

Multigrid as a solver for the differential problem

From a more sophisticated point of view, multigrid can be regarded as a solution method for the (continuous) differential problem. In this view it is not appropriate to separate the discretization and the solution of the discrete problem, but rather to regard both processes as interdependent: the solution process can, according to this view, be performed the more efficiently the more the continuous differential background is exploited.

One simple, but very natural way of looking at multigrid as a differential solver is represented by the full multigrid method (FMG, see Section 2.6). Here, the method is oriented to minimizing the differential error rather than to minimizing the algebraic error (corresponding to the linear system). Self-adaptive grid refinements and related approaches, which are very natural in the multigrid context, also belong to the differential view of multigrid (and will be presented in Chapter 9).

Efficiency

Multigrid methods are highly efficient. Efficiency here relates to both a theoretical feature and a practical one. The theoretical efficiency is more precisely expressed by the term optimality in a complexity theory sense. If interpreted appropriately, (full) multigrid methods are typically optimal in the sense that the number of arithmetic operations needed to solve a (discrete) problem is proportional to the number N of unknowns in the problem considered.

Efficiency in the practical sense means that the proportionality constants in this O(N) statement are small or moderate. This is indeed the case for multigrid: if designed well, the h-independent convergence factors can be made very small (in the range 0.1–0.2 or even less) and the operation count per unknown per iteration step is also small. If a concrete multigrid algorithm has this property, we speak of the typical multigrid efficiency (sometimes also called top multigrid efficiency).

Generality

The second striking property of multigrid is its generality. Multigrid methods are as efficient as the so-called fast elliptic solvers [95, 342] but are less limited in their range of application. Multigrid methods can be applied with full efficiency

– to general elliptic equations with variable coefficients,

– in general domains,

– for general boundary conditions,

– in 2D, 3D and higher dimensions (trivially also for 1D problems),

– to scalar equations and to systems of equations.

Very importantly, multigrid can also be applied directly, i.e. without global linearization, to nonlinear elliptic problems.

Furthermore, multigrid is not restricted to a certain discretization approach. In principle, multigrid can be used in connection with any type of grid-based discretization: finite differences, finite volumes and finite elements. (Collocation, variational, spectral or particle-based discretization methods can also be combined with the multigrid principle. In this book, however, we will concentrate on grid-type approaches.)

Optimality versus robustness

Multigrid methods are characterized by their so-called components. We will introduce them systematically in Sections 2.1–2.4. The components are the smoothing procedure, the coarsening strategy, the coarse grid operators, the transfer operators from fine grids to coarse and from coarse to fine and the cycle type. These components have to be specified for each concrete problem. Although it is well known how to choose suitable multigrid components for large classes of problems, it may be (very) difficult to define the right or reasonable ones in complicated new applications. This is still an art requiring theoretical insight, experience and numerical experiments. There are two trends with respect to the choice of multigrid components.

In optimized multigrid algorithms, one tries to tailor the components to the problem at hand in order to obtain the highest possible efficiency for the solution process. This optimized approach clearly makes sense if a very large scale problem has to be (repeatedly) solved or if a smaller core problem needs to be solved many times every day, like, for example, the 3D Helmholtz equation in daily weather prediction [251].

On the other hand, the idea of robust multigrid algorithms is to choose the components independently of the given problem, uniformly for as large a class of problems as possible. The robust approach is often used in software packages where the highest efficiency for a single problem is not so important. The AMG method in Appendix A is an example of a robust multigrid method.

Optimization of multigrid is in many cases limited by practical circumstances. For instance, important large software packages are based on certain grid structures, discretization techniques and solution methods, which are not oriented to multigrid requirements. In such situations, the users and developers of the software packages may be very interested in accelerating the program by introducing some multigrid features. It may, however, be impossible or too costly to change the data structures in such a way that optimal multigrid features can be integrated. From a puristic multigrid point of view, such an accelerating approach may be very unsatisfactory because much higher efficiency is achievable in principle. From a practical point of view, this approach can, however, be a good compromise. A simple modification of the program may give a significant reduction in computing time.

Adaptivity

Defining a global grid for the discretization of a given problem independently of the solution process is often insufficient. Only during the solution process may certain features of the solution, like shocks, singularities, oscillations, turbulent behavior and the like be recognized. In such cases, the local discretization error is of different size in different parts of the domain and therefore it would be natural to adapt the grids (and perhaps also the discretization) to the behavior of the solution. This is one of the main reasons for using adaptive grids that are dynamically constructed within the solution process. For 3D problems, this is particularly necessary since a grid that is globally as fine as is needed at some crucial parts of the domain is not affordable.

For these reasons, adaptivity of grids is one of the major trends in today’s numerical simulation and scientific computing.

Adaptivity can be combined with the multigrid principle in a very natural way. In the adaptive multigrid process finer and finer grids are not constructed globally. They are only constructed in those parts of the domain where the current discretization error is significantly large. Essentially, all other multigrid components are maintained as usual. What is specifically needed for this approach, are criteria to measure or estimate the current local discretization error. Adaptive multigrid will be discussed in Chapter 9.

Parallel features

A promising and challenging trend in numerical simulation and scientific computing is the use of parallelism in numerical algorithms. The background to this trend is the fact that most high performance computers are now parallel systems.

In order to solve any problem on a parallel computer, an algorithm with an appropriate degree of parallelism has to be made available. In addition, the data being processed in the algorithm should be suitably organized. This is particularly important if the parallel computer has distributed memory. In this case, the communication overhead has to be limited, i.e. the internal transfer of data should cost only a (small) fraction of the overall effort for solving the given problem.

In Chapter 6, we deal mainly with the practical aspects of parallel multigrid, focusing on the grid partitioning approach, in which the given grid is partitioned into a number of subgrids. Each processor works on its subgrid but has to communicate with the other processors. In this approach all multigrid components that have to be specified should be parallel, or as parallel as possible, but at the same time as efficient as possible.

1D problems

The multigrid idea leads to optimal algorithms for many PDEs. Since ordinary differential equations (ODEs) are a special case of PDEs, the question arises whether multigrid methods are also useful for ODEs. Indeed, the multigrid principle can also be applied to ODEs and, of course, leads to optimal (i.e. O(N)) algorithms in 1D. However, in 1D other optimal methods are available and multigrid methods typically coincide with (and degenerate to) well-known optimal solvers. In that respect, the multigrid principle is also applicable to the 1D case, but not really needed there, at least not for standard problems.

For example, if linear ordinary boundary value problems are discretized with standard discretization methods (for example, finite differences) band matrices are obtained, where the bandwidth is independent of the mesh size h. In particular, if a three-point discretization is used for the differential operator, typically the discrete problem is characterized by a linear system with a tridiagonal (N × N) matrix. Such tridiagonal systems, and more general band matrices with a fixed small bandwidth, can be solved in O(N) operations by Gaussian elimination-type methods.

1.5.5 Multigrid History

The forerunners of multigrid are the ideas of nested iteration, error smoothing by relaxation and total reduction. Nested iteration (see also Section 2.6) has been used for a long time to obtain first approximations (initial guesses) on fine grids from coarse grids, for instance in the context of Newton’s method. Also, the fact that in many cases relaxation processes have an error smoothing property has been known for a long time [367, 368, 374]. The total reduction method by Schröder and Trottenberg [342343–344] has the complete multigrid structure. The main difference from standard multigrid is, however, that on the coarse grids equations are constructed, which are equivalent to the fine grid equations for the respective unknowns; smoothing plays no role then.

The first studies investigating multigrid methods in a strict sense were those by Fedorenko [139, 140] (1962, 1964) and that of Bakhvalov [16] (1966). While Fedorenko [140] restricted the convergence investigation to a discrete boundary value problem of second order with variable coefficients (in the unit square), Bakhvalov also indicated the possibility of combining multigrid methods with nested iteration.

Though the studies published by Fedorenko and Bakhavalov had shown the asymptotic optimality of the multigrid approach (and to a certain extent its generality), their actual efficiency was first recognized by Brandt in the early 1970s. Studying adaptive grid refinements and their relation to fast solvers, Brandt had been led to the papers of Fedorenko and Bakhvalov by Widlund. In his first two papers [56, 57] (1973, 1976) and then summarized in the systematic 1977 work [58], Brandt showed the actual efficiency of multigrid methods. His essential contributions (in the early studies) include the introduction of a nonlinear multigrid method (FAS) and of adaptive techniques (MLAT), the discussion of general domains, the systematic application of the nested iteration idea (FMG) and, last but not least, the provision of the local Fourier analysis tool for theoretical investigation and method design.

The following papers, which we would like to mention as being historically relevant contributions, are representative of the early development of multigrid.

In 1971 Astrakhantsev [6] generalized Bakhvalov’s convergence result to general boundary conditions. Like Bakhvalov he used a variational formulation in his theoretical approach. After a first study of multigrid methods for Poisson’s equation in a square, Nicolaides [278] (1975) discussed multigrid ideas in connection with finite element discretizations systematically in [279] (1977). In 1975, Frederickson [142] introduced an approximate multigrid-like solver which can be regarded as a forerunner of the MGR methods [141, 319].

In 1975–1976, Hackbusch developed the fundamental elements of multigrid methods anew. It was again Widlund who informed Hackbusch about the studies which were already available. Hackbusch’s first systematic report [169] (1976) contained many theoretical and practical investigations which have been taken up and developed further by several authors. So one finds Fourier analysis considerations, the use of red–black and four–colour relaxation methods for smoothing, the treatment of nonrectangular domains and of nonlinear problems etc. Hackbusch then presented a general convergence theory of multigrid methods [170–172].

Since the early 1980s, the field has been exploding and many researchers have contributed to the field. Two series of conferences dedicated to multigrid were set up. The European Multigrid Conferences (EMG) have been held at Cologne (1981, 1985), Bonn (1991), Amsterdam (1993), Stuttgart (1996) and Ghent (1999). And in the US the Copper Mountain Conferences on multigrid have been held bi-annually since 1983. Proceedings of the European meetings have appeared in [129, 174, 177, 178, 181, 189] and of the Copper Mountain Conferences in special issues of journals: Applied Numerical Mathematics (Vol. 13, 1983; Vol. 19, 1986), Communications in Applied Numerical Methods (Vol. 8, 1992), SIAM Journal of Numerical Analysis (Vol. 30, 1993), Electronic Transactions on Numerical Analysis (Vol. 6, 1996). Another rich source of information on multigrid is the MGNet website maintained by C.C. Douglas: http://www.mgnet.org. This website includes a large multigrid bibliography with more than 3000 entries. Some multigrid textbooks and monographies are [54, 66, 91, 176, 206, 262, 332, 351, 378, 415].

1.6 INTERMEZZO: SOME BASIC FACTS AND METHODS

This section contains some general numerical material which is needed at several places of the book. One reason for presenting this material here is to clarify terminology. The reader who has a general idea of the material presented here may have a quick look at it and return to it for details later. Those readers for whom the material is new may read it now or postpone its study. (We will refer these readers back to these sections later.)

1.6.1 Iterative Solvers, Splittings and Preconditioners

Since the facts listed in this section are valid for general matrices and are not restricted to discrete differential operators, we use linear algebra notation in this section, i.e. matrices Ah or A instead of discrete operators Lh.

Let

     (1.6.1)

be a linear system with an invertible matrix A. The simplest iterative scheme for this equation is the Richardson iteration

     (1.6.2)

with some acceleration parameter τ ≠ 0.

A more general iteration is

     (1.6.3)

Here, M is the iteration matrix. We assume that the original equation Au = f is equivalent to the fixed point equation u = Mu + s. For Richardson’s iteration, we have M = I − τA, s = τf.

The convergence (and the asymptotic convergence speed) of Richardson’s iteration and of the general iteration are characterized by the spectral radii ρ(I − τA) and ρ(M), respectively. The spectral radius of a matrix M is defined as

     (1.6.4)

The spectral radius is the asymptotic convergence factor of the iteration. Asymptotically