Dianne P O'Leary.

Capacitance matrix methods for the Helmholtz equation on general three-dimensional regions online

. (page 1 of 7)
Font size
QR-code for this ebook


Courant Mathematics and
Computing Laboratory

U. S. Department of Energy

Capacitance Matrix Methods

for the Helmholtz Equation on
General Three-Dimensional Regions

Dianne P. Oleary and Olof Widlund

U. S. Department of Energy Report

Prepared under Contract EY-76-C-02-3077
with the Office of Energy Research

Mathematics and Computing
October 1978

New York University


Courant Mathematics and Computing Laboratory-
New York University

Mathematics and Computing COO-3077-155


Dianne P. O'Leary and Olof Widlund

Department of Mathematics, University of Michigan, Ann Arbor,
Michigan. The work of this author was supported in part by
a National Science Foundation Grant MCS76-O6595, while at
the University of Michigan. Current address: Department of
Computer Science, Univ. of Maryland, College Park, Maryland.

Courant Institute of Mathematical Sciences, New York University,
The work of this author was supported by ERDA, Contract No.
EY-76-C- 02-3 077* 000.

U. S. Department of Energy
Contract EY-76-C-02-3077*000


1. Introduction

It is well known that highly structured systems of linear
algebraic equations arise when Helmholtz's equation

(l.l) -Au+cu = f , c = constant ,

is discretized "by finite difference or finite element methods using
uniform meshes. This is true, in particular, for problems on a
region O which permits the separation of the variables. Very fast
and highly accurate numerical methods are now readily available to
solve separable problems at an expense which is comparable to that
of a few steps of any simple iterative procedure applied to the
linear system; see Bank and Rose [2,3], Buneman [5], Buzbee, Golub
and Nielson [8], Fischer, Golub, Hald, Leiva and Widlund [l6],
Hockney [24,26], Swarztrauber [50,51], Swarztrauber and Sweet
[52^53] and Sweet [54]- Adopting common usage, we shall refer to
such methods as fast Poisson solvers.

The usefulness of these algorithms has been extended in
recent years to problems on general bounded regions by the develop-
ment of capacitance matrix, or imbedding, methods; see Buzbee and
Dorr [6], Buzbee, Dorr, George and Golub [7], George [19], Hockney
[25,27], Martin [35], Polozhii [4o], Proskurowski [41,42,43],
Proskurowski and Widlund [44,45], Shieh [46,47,48] and Widlund
[57] • We refer to Proskurowski and Widlund [44] for a discussion
of this development up to the beginning of 1976. All of the
numerical experiments reported in those papers were carried out
for regions in the plane. Strong results on the efficiency of

certain of these methods have been rigorously established through
the excellent work of Shieh [46,47,48]. Algorithms similar to
those which we shall describe have recently been implemented very
successfully for two-dimensional regions by Proskurowski [42,43]
and Proskurowski and Widlund [45]. In that work, a new fast
Poisson solver, developed by Banegas [1 ], has been used exten-
sively; see Section 5. We note that the performance of computer
programs implementing capacitance matrix algorithms depends very
heavily on the efficiency of the fast Poisson solver, and if
properly designed, they can be easily upgraded by replacing that
module when a better one becomes available.

In this paper, we shall extend the capacitance matrix method
to problems in three dimensions. The mathematical framework, using
discrete dipole layers in the Dirichlet case, is an extension of
the formal discrete potential theory developed in Proskurowski and
Widlund [44] . We note that these algorithms must be quite
differently designed in the three-dimensional case. As in two
dimensions the fast Poisson calculations strongly dominate the
work. The number of these calculations necessary to meet a given
tolerance remains virtually unchanged when the mesh size is
refined. We have developed a FORTRAN program for Cartesian co-
ordinates and the Dirichlet problem, which turns out to be techni-
cally more demanding than the Neumann case. This program has been
designed to keep storage requirements low. The number of storage
locations required is one or two times N, the number of mesh points
in a rectangular parallelepiped in which the region is imbedded,
and a modest multiple of p, the number of mesh points which belong

to the region O and are adjacent to its boundary. A further sub-
stantial reduction of storage can be accomplished for very large
problems by using the ideas of Banegas [1], see further Section 5.

In the second section, we discuss the imbedding idea.
Following a review of classical potential theory, we derive our
capacitance matrix methods in Section 3- Section 4 focuses on
algorithmic aspects which are of crucial importance in the develop-
ment of fast, reliable and modular computer code. We solve the
capacitance matrix equations by conjugate gradient methods. These
methods, originally used in a similar context by George [19], are
reviewed in that section. We also discuss how spectral information
and approximate inverses of the capacitance matrices can be
obtained and used at a moderate cost in computer time and storage.
The fast Poisson solver which is used in our program is described
in Section 5. It is numerically stable even for negative values of
the coefficient c of the Helmholtz operator. Finally, we give
details on the organization of our computer program and results
from numerical experiments. These tests were designed to be quite
severe and the method has proved efficient and reliable.

A listing of our program is provided as an appendix. It has
been checked by the CDC ANSI FORTRAN verifier at the Courant
Mathematics and Computing Laboratory of New York University. It
has been run successfully on the CDC 660O at the Courant Institute,
a CDC 7600 at the Lawrence Berkeley Laboratory and the Amdahl
47OV/6 at the University of Michigan.

Acknowledgements . The authors want to thank John G. Lewis,

wTodzimierz Proskurowski and Arthur Shieh for their interest and
help with this project.


2. Discrete Helmholtz Problems and. Imbedding

2.1. The Imbedding of the Discrete Problem

In this section, we shall discuss how discretizations of the

-Au+ cu = f on O ,

with a boundary condition and data given on hO, can be imbedded in
problems for which fast Poisson solvers can be used. In the second
subsection, we describe in detail how these ideas apply to the
finite difference scheme which we have used in our numerical
experiments .

The efficiency of capacitance matrix methods depends on the
choice of appropriate finite difference and finite element meshes.
Interior parts of the mesh should be made regular in the sense that
the linear equations at the corresponding mesh points match those
of a fast Poisson solver. We denote the set of these mesh points
t)y O, where h is a mesh width parameter. The set of the remaining,
irregular mesh points is denoted by ^O, • These points are typi-
cally located on or close to the boundary ^O and the discrete equa-.
tions associated with them are computed from local information on
the geometry of the region. For efficiency, the number of unknowns
associated with the points ia So, should be kept small, since the
equations and other information required at the regular mesh points
are inexpensive to generate and can be stored in a very compact

If we work in Cartesian coordinates it is natural to imbed our
open, bounded region O in a rectangular parallelepiped and to use

a rectangular mesh. Other choices which permit the separation of
the variables on the larger region, can equally well be chosen. On
the larger region a mesh suitable for a fast Poisson solver is intro-
duced which coincides with the regular part of the mesh previously
introduced for the region Q. The position of the larger region
relative to O is largely arbitrary but when using discrete dipoles
(see Section 3 )j we need a layer of exterior mesh points, one mesh
width thick, outside of Ov.^U^O, . We shall use some or all of the
discrete equations at exterior mesh points to expand our original
linear system into one which is of the same size as the one which
is solved by the fast Poisson solver. The set of mesh points
corresponding to these equations is denoted by CO-..

Before we describe how these larger systems of equations
are derived, we shall show by two examples how these sets of mesh
points can be constructed. We first consider a Dirichlet problem
solved by a classical finite difference scheme on a rectangular
mesh. The values of the approximate solution are sought at the
mesh points which belong to O. The discretization of the Helmholtz
operator on the larger region induces, for each mesh point, a
neighborhood of points used by its stencil. A mesh point in O be-
longs to n, if and only if all its relevant stencil neighbors are

in O.and hO, is the set of the remaining mesh points in O. The set
' h

CO, is the set of all mesh points which belong to the complement
of O- It thus includes any mesh point which is on the boundary hO.
As a second example, consider a Neumann problem for Laplace's
equation in two dimensions solved by a finite element method with
piecewise linear trial functions. The region is approximated by a


union of triangles using a regular triangulation, based on a uni-
form mesh, in the interior of the region. The set O, will then
correspond to the set of equations which are not affected by the
particular geometry of the region. Values of the discrete solution
are also sought at the vertices on the boundary. These points
normally fail to lie on a regular mesh. They belong to So, to-
gether with certain mesh points which are close to the boundary.
Each irregular point can be assigned to a close-by mesh point of
the regular mesh which covers the larger region and we then define
CO, as the set of remaining, exterior mesh points. There are a
number of permissible ways in which this assignment can be made.
Similar constructions can be carried out for higher order accurate
finite element methods; see Proskurowski and Widlund [45] for
further details.

Let us write the expanded linear system in the form

(2.1) Au = b

where u is the vector of values of the discrete solution at the
mesh points and the components of b are constructed from the func-
tion f and the data given on hO. By construction, our formulas for
the interior and irregular mesh points do not involve any coupling
to exterior mesh points, and the matrix is therefore reducible,
i.e. there exists a permutation matrix P such that

P AP = ^


■ A A
V 21 22


The block matrix A-,-, represents the approximation of the probl


on O, U^, . It is clear from the structure of this system that the
restriction of the solution of the system (2.1) to this set is
independent of the solution and the data at the exterior points.
Our methods also produce values of a mesh function for the points
of COr^ hut they are largely arbitrary and useless. Similarly, we
must provide some extension of the data to the set CO, , but the
performance of the algorithms is only marginally affected by this

Let B denote the matrix representation of the operator
obtained by using the basic discretization at all the mesh points.
Only those rows of A and E which correspond to the irregular mesh
points differ provided the equations and unknowns are ordered in
the same way. We can therefore write

A = B + UZ^ ,

where U and Z have p columns, with p equal to the number of ele-
ments of the set ^O,. It is convenient to choose the columns of U


to be unit vectors in the direction of the positive coordinate axes
corresponding to the points of ^O, . The operator U is then an
extension operator which maps any mesh function, defined only on
bo, , onto a function defined on all mesh points. The values on So,

are retained while all the rer^aining values are set equal to zero.

The transpose of U, U , is a restriction, or trace, operator which

maps any mesh function defined everywhere onto its restriction to

bo, . The matrix Z can, with this choice of U, be regarded as a

compact representation of A-B, obtained by deleting the zero rows

corresponding to the equations for the mesh points in O, U CO, • It
is important to note that Z and U are quite sparse, a reflection
of the sparsity of A and B.

In Sections 3 and h, we shall discuss efficient and stable
ways of solving the linear system (2.1).

2.2. The Shortley-Weller Scheme

We shall now discuss the finite difference scheme which has
been used in our numerical experiments to solve the Dirichlet
problem and also describe how the necessary information on the
geometry of the boundary is handled.

The second order accurate Shortley-Weller formula (see
Collatz [9], Chap. 5.1 or Forsythe and Wasow [I7], Sec. 20.7) can
be understood as the sum of three point difference approximations
for the second derivative with respect to each of the three inde-
pendent variables. The value at the nearest mesh neighbor in each
positive and negative coordinate direction is used unless this
neighbor belongs to the set CO, . In that case the Dirichlet data
at the point of intersection of the mesh line and the boundary is

As an example, suppose that the mesh spacings in the x, y and
z directions are all equal to h. Consider an irregular mesh point,
with indices (i,j,k), which has two exterior neighbors in the x
direction and one in the positive y direction. Let 5 , 5 and

— X "tX

5 be the distances to the boundary, in the respective coordinate
directions, measured in units of the mesh size h and let g ? g,
and g be the Dirichlet data at th_ corresponding points on the

boundary hO. Then our approximation to -Au +cu = f at this
irregular point is.

(2/(5^ 5 ) + 2/6 , + 2 + ch^)u. .,

- (2/(l+5^y))u,^j_i^i, - u. .^^^^ - u. .^^_^

h^f . ., + (2/(5^ + 5^ 5 ) )g^
ijk ^ '^ ^ +x +x -x' ^^+x

+ (2/(5^ + 5^ 5 ))g + (2/(5^ + 5^ ) )g^

^ ' ^ -X +x -X ^-x +y +y +y

At the regular points the formula reduces to a simple seven point

The Shortley-Weller formula has a matrix of positive type.
This permits the use of the classical error estimates based on a
discrete maximum principle, as in the references given above. The
only information required on the geometry of the region is the
coordinates of the irregular mesh points and the distances along
the mesh lines from each such point to the boundary. This appears
to be close to the minimum information required by any method with
more than first order accuracy. See Proskurowski and Widlund [44],
Pereyra, Proskurowski and Widlund [39] and Strang and Fix [49] for
more details. This geometrical information is also sufficient to
construct higher order accurate approximations to the Helmholtz
equation, as in Pereyra, Proskurowski and Widlund [39] where a
family of methods suggested by Kreiss is developed. These methods
have proven quite effective for two dimensional problems but their
usefulness is limited by the requirement that each irregular mesh


point must have several interior mesh neighbors along each mesh
line. This requirement is met by shifting the. region and refining
the mesh if necessary. Although this is practical in two dimen-
sions, it is much more difficult for three dimensional regions.

We are free to scale the rows of the matrix A which corre-
spond to the irregular mesh points. The choice of scaling is
important since it affects the rate of convergence of our iterative
method. Based on the analysis given in the next section, the
experience in the two dimensional case (see Proskurowski and
Widlund [-^"^ ] ) and our numerical experiments, we have chosen to make
all diagonal elements of A equal to one.


3 . Potential Theory and Discrete Dlpoles

3.1. The Continuous Case

In this section, we shall give a brief survey of certain
results of classical potential theory and also develop an analo-
gous, formal theory for the discrete case. We shall mainly follow
the presentation of Garabedian [18] when discussing the continuous
case, specializing to the case of c = 0. A discrete, formal theory
has previously been developed by Proskurowski and Widlund [44] but
our presentation in Sections 3-2 - '^ .'k will be more complete in
several respects.

We first introduce the volume, or Newton, potential

(3.1) u^(x) = (l/^ir) J f(0/r d|

where x - {x-^,x^,^), ^ = {^^,i^,^) and r = ((x^- ^^) + [x^ - i^)
+ (x^ -1^)^) • We note that (l/4Tr)(l/r) is a fundamental solu-
tion of the operator -A, i.e.,

-Au^ = f .

A single layer potential, with a charge density p, is given by,

(3.2) %^{x) = {l/2jr) J p(|)/r da

and a double layer potential, with a dipole moment density \x, by

(3.3) VM = {l/2Tr) J n(|)(V^v^)(l/r)da .

Here v denotes the normal of the bour.aary hO directed towards the


interior of O- By V and "V' , we denote the limits of ">' when
the boundary is approached from the outside and inside respectively
ajid similar notations are also used for the limits of '/r . The func-
tions -2^ and '^/^ are real analytic functions in the complement
of ho. By using a Green's formula one can establish that '2' and
h'^^y^v are continuous and that jiomp conditions hold for S'?^/Sv
and '%'' ; see Garabedian [18], Chapter 9' Thus, for a region with a
smooth boundary,

hV^'-^/hv - ( + ) p + il/2ir) J p(V^v^)(l/r)da ,


>r^^^ = (-) H + {l/2ir)J p(V^v^)(l/r)do ,


With the aid of these relations the Neumann and Dirichlet problems
can be reduced to Fredholm integral equations. For the interior
Neumann problem,

-Au = f in O ,

hu/hv = gp. on So ,

we make the Ansatz,

u(x) = u^(x) +'V{^) .
The boundary condition is satisfied by choosing p such that


h-ryhv = -p + (l/27r) j"p(V^v^)(l/r)da


(3.4) " % ~ (V^)u^


This equation can be written as (l-K)p = -g, where K is a compact
operator defined by the formula above. It is a Fredholm integral

equation of the second kind with a simple zero eigenvalue. Since

2 2

K is compact in L the integral operator I-K is bounded in L and

it has an inverse of the same form on a space of codimension one.

Equation (3''^) is solvable if g is orthogonal to the left eigen-

f unction of (l-K) corresponding to the zero eigenvalue. In this

case this simply means that g should have a zero mean value. By

using the same Ansatz for the exterior Neumann problem, we obtain

an integral equation with the operator I+K.

If we use the same single layer Ansatz for the interior

Dirichlet problem, with data g , we get an integral equation of the

first kind,

(l/27r) / p/r do = g^ - u,.

(l/27r) j p/r do = gj^ - u^


This operator does not have a bounded inverse in Lp. The use of
an analogous Ansatz for the discrete Dirichlet problem gives rise
to capacitance matrices which become increasingly ill-conditioned
as the mesh is refined.
The Ansatz

u(x) = u^(x) +Mx) ,
which employs a double layer potential, leads to a Fredholm


integral equation of the second kind?

'}f- = ix + {l/2Tr) J [iih/hv^Hl/r)da

(3.5) =^^^-^U*

The integral operator is now I+K , where K is the transpose of the

operator introduced when solving the Neumann problems. We shall

obtain well-conditioned capacitance matrices when using a discrete

analogue of this approach.

The close relationship between the integral equations for the
interior Dirichlet and exterior Neiimann problems is used to
establish the solvability of the Dirichlet problem; see Garabedian
[18], Chapter 10. A similar argument is given in Section 3.3 for
a discrete case.

The integral operator K is not symmetric except for very
special regions. Nevertheless it has real eigenvalues; see e.g.
Kellogg [32], p. 309. For future reference, we also note that
there exist variational formulations of the Fredholm integral equa-
tions given in this section; see Nedelec and Planchard [37] . It
can be shown that the mapping defined by the single layer potential
'Y is an isomorphism from H~ ''^ (^)/P to the subspace of ir"(o)/P
of weak solutions of Laplace's equation. Here H (o) is the space
of functions with square integrable first distributional deriva-
tives, H / (hn) the space of traces of H (o), H~ ^ Cdo) the space

dual to H ' (^O), and P the space of constants. By substituting

the single layer potential into the standard variational formula-
tion of the interior Neumann problem and using a Green's formula.


an alternative formulation is obtained. The resulting bilinear
form is coercive on H" "'"'^^ ( So ) /P^ and is equivalent to equation


Before we turn to the discrete problems, we note that, in the

theory just developed, the function (l/47r)(l/r) can be replaced by

other fundamental solutions of the Laplace operator. In particular,

we can use a Green's function for a rectangular parallelepiped in

which the region O is imbedded. The theory can also be extended,

in a straightforward way, to Helmholtz's equation with a nonzero

coefficient c.

3.2. Discrete Potential Theory

We now return to the solution of Au = b, (equation (2.1))
with A = B +\]Z^ . Guided by the theory for the continuous case, we
shall develop two algorithms, one suitable for the Neumann and the
other for the Dirichlet case.

We shall assume that B is invertible. This is not a very
restrictive assumption since we have a great deal of freedom to
choose the boundary conditions on the larger region.

We recall from Section 2.1 that the columns of U were chosen
to be unit vectors corresponding to the irregular mesh points. If
we order the points of Oj^ first, followed by those of ho^ and CO^,
we can obtain the representation.

U =




where I is a pxp identity matrix. Let us, in analogy to the con-
tinuous case, make the Ansatz

(3.6) u = Gb + GWs

where the vector s has p components, G is the inverse of B, and W
has the form

W =


The operator G plays a role very similar to that of ^a fundamental
solution for the continuous problem. The second term GWs corre-
sponds to a single or double layer potential. For additional
flexibility, we have introduced the mesh function b which coincides
with b except possibly at the irregular points of So, • In particu-
lar, if the Helmholtz equation has a zero right hand side, we can
often choose b = 0, eliminating the first term of the Ansatz. To
arrive at an equation for the vector s, we calculate the residual,

b-Au = b -(B+ UZ^ ) ( Gb + GWs )

= (b-b)-UZ^Gb- (l+UZ^G)Ws .

From the form of b, U, and W, we have the following result:

Lemma 3 ♦ 1 • The residuals for the system (2.1) corresponding to the

points of O, are zero for any choice of the vector s in (3 •6). If

the matrix W-;, is zero they also vanish at all points of GO. •

We now demand that the residuals vanish on the set ^Q, :


- U^(b -Au) = U^(b-b)- Z^Gb -U^AGWs .


This gives us a system of p equations:

(3.7) Cs = U^AGWs = (u\+z'^G¥)s = U^(b-b) - Z^Gb ,

where C is the capacitance matrix. We ignore the residuals

on the set CO, since the extension of the data to this set is

largely arbitrary. It follows from the reducible structure of A
that if the capacitance matrix C is nonsingular the restriction of
the mesh function u, given by formula {^.6), solves the discrete
Helmholtz equation. We shall now discuss two choices of the matrix
W and study the invertibility of the resulting matrices.

For a Neumann problem, our choice of W should correspond to

a single layer Ansatz. We therefore choose W = U and note that the

capacitance matrix G = U AGU is then the restriction of AG to the

subspace corresponding to the set SOj^. Using equations (3.6) and

(3.7), we find,

~ /T x-l/T~ T/ ~\x

u = Gb -GU(U AGU) (Z-^Gb - U (b-b ) ) .

This is, for b = b, the well known Woodbury formula; see
Householder [29]. For completeness, we give a proof of the
following result.

Theorem 3.1. The capacitance matrix C„ is singular if and only if
the matrix A is singular. For b = b the equation (3.7) fails to
have a solution if and only if b does not lie in the range of A.

Proof: Let cj) be a nontrivial element of the null space of C

Then since C„ = I + Z GU, the vector

Z^GUcj) = -cj)



is nonzero and therefore GUif) cannot vanish identically. But

AGUcj) = UC„(f) = and therefore A is singular. Let now ijj belong to

the null space of C^ and assume that

^^(Z^Gh) = (^^Z^G)b ^ .

Then b does not belong to the range of A since

1 3 4 5 6 7