COO-3077-155

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

UNCLASSIFIED

Courant Mathematics and Computing Laboratory-

New York University

Mathematics and Computing COO-3077-155

CAPACITANCE MATRIX METHODS FOR THE HELMHOLTZ EQUATION

ON GENERAL THREE-DIMENSIONAL REGIONS

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,

');
}
elseif (getClientWidth() > 430)
{
document.write('');
}
else
{
document.write('');
}
//-->

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

UNCLASSIFIED

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.

4

2. Discrete Helmholtz Problems and. Imbedding

2.1. The Imbedding of the Discrete Problem

In this section, we shall discuss how discretizations of the

problem

-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

form.

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

6

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

T

P AP = ^

f^ll

■ A A

V 21 22

7

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

em

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

choice.

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

h

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.

T

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

used.

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

approximation.

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

10

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.

11

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

12

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 ,

bo

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

ho

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

13

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

an

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

ao

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^

ha

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

14

integral equation of the second kind?

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

(3.5) =^^^-^U*

T T

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

1/2

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.

15

an alternative formulation is obtained. The resulting bilinear

form is coercive on H" "'"'^^ ( So ) /P^ and is equivalent to equation

(3.4).

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 =

^0^

^°y

16

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 =

^0^

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, :

h

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

17

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

h

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

T

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

T

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

Z^GUcj) = -cj)

N'

18

is nonzero and therefore GUif) cannot vanish identically. But

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

T

the null space of C^ and assume that

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

Then b does not belong to the range of A since