# A Fourier Continuation Method for the Solution of Elliptic Eigenvalue Problems in General Domains

Authors: Oscar P. Bruno; Timothy Elling; Ayon Sen
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

We present a new computational method for the solution of elliptic eigenvalue problems with variable coefficients in general two-dimensional domains. The proposed approach is based on use of the novel Fourier continuation method (which enables fast and highly accurate Fourier approximation of nonperiodic functions in equispaced grids without the limitations arising from the Gibbs phenomenon) in conjunction with an overlapping patch domain decomposition strategy and Arnoldi iteration. A variety of examples demonstrate the versatility, accuracy, and generality of the proposed methodology.


---

## Body

## 1. Introduction

We present a new computational method for the solution of variable-coefficient elliptic eigenvalue problems in general two-dimensional domains. This method uses Arnoldi iteration to approximate the eigenvalues of numerical differential operators in such a way that the bulk of the associated computational cost arises from the needed repeated evaluations of the differential operator for given input functions (vectors). For accuracy and efficiency our algorithm produces such evaluations on the basis of the Fourier continuation (FC) method [1] which, using equispaced meshes, resolves the inaccuracies due to the Gibbs phenomenon by computing smooth, periodic extensions of arbitrary smooth functions, therefore allowing for accurate evaluation of derivatives in the frequency domain. The efficiency of the method thus results from its reliance on the Fast Fourier Transform (FFT) for numerical differentiation.Differential eigenvalue problems arise in many areas of physical science; a prominent example is the Laplace eigenvalue problem:(1)
∇
2
u
x
,
y
=
λ
u
x
,
y
.In the simplest physics perspective the eigenvalues correspond to the fundamental modes of vibration of a thin membrane with geometry corresponding to the domain of (1). Solutions to (1) completely characterize membrane motions since the eigenfunctions form a complete orthonormal basis for the space of membrane vibrations. Other physically relevant applications of the Laplace eigenproblem include the description of quantized energy states of particles (under various potentials) and the dynamics of electromagnetic waves traveling through waveguides [2]. Chain reactions in nuclear reactors are described by a more complicated eigenvalue problem, in which the smallest-magnitude eigenvalue describes whether the reaction is subcritical, critical, or supercritical [3]. Because natural processes do not only occur in simple separable geometric regions such as squares or disks, the only types of regions on which (1) can be solved analytically, much effort has been invested in the solution of eigenvalue problems on general domains.Existing approaches for differential eigenvalue problems include collocation methods as well as finite element or finite difference methods. Collocation methods utilize expressions for the eigenfunction in terms of adequately chosen bases. A classical example is the method of particular solutions, in which the coefficients in the expansion are obtained via consideration of a homogeneous system of linear equations that enforces the required (homogeneous) boundary conditions at points placed along the domain boundary [4, 5]. A noteworthy aspect of this method is its use of singular (Fourier-Bessel) functions in the expansion to incorporate, for example, reentrant corners. Improvements have been made by including points on the interior of the domain where the condition that the eigenfunction is nonzero is enforced [5]. Similar approaches were put forth in [6, 7] with different choices for the basis functions. Other generalizations include the use of conformal mappings to recast the problem in a circular domain as a generalized eigenvalue problem, as is done in [8] for domains with approximately fractal boundaries.Other eigenvalue solvers utilize finite difference or finite element approximations which reduce differential eigenvalue problem to a corresponding finite-dimensional algebraic eigenvalue problem. Finite element methods rely on unstructured meshes that can be generated and adaptively refined to account for complicated domain features [9]. Domain decomposition techniques are often used in conjunction with finite elements in order to deal with discontinuous coefficients or to reduce the problem to a number of smaller problems [10], as are iterative eigenvalue algorithms [11]; the algorithm proposed in this paper utilizes both domain decomposition and iterative methods. An extension of these ideas can be found in the method given by [12], which uses the basis from the method of particular solutions and domain decomposition, allowing for local usage of the basis near singularities. Results from a wide array of these types of methods are summarized in [2, 13].Even though collocation methods have demonstrated highly accurate evaluation of eigenvalues for some singular geometries (such as the L-shaped membrane), the applicability of these methods is limited to locally separable geometries. As indicated above, use of conformal maps has been suggested in the literature as a remedy for this difficulty [6], but determining a mapping from a simple shape to an arbitrary region is a nontrivial task. Finite element methods can also be effective, and they are applicable to rather general configurations. Finite difference methods, finally, are efficient by virtue of their finite bandwidths, but this advantage turns into a challenge near domain boundaries. Even high-order centered finite difference schemes must be reformulated near boundaries, resulting in an asymmetric stencil and generally a lower order of accuracy at the boundary [14].To address the difficulties arising in the treatment of general elliptic eigenvalue problems we propose a methodology wherein derivatives are produced by means of the FC(Gram) algorithm [1]. Differentiation of Fourier series is, of course, straightforward and, in view of the Fast Fourier Transform, can be produced in O
(
N
log
⁡
N
) operations, where N is the number of points in the mesh. Because the FC(Gram) algorithm produces smooth periodic extensions for given functions, Fourier series truncated to a small number of terms are still highly accurate representations of the continued function: the difficulties arising from the Gibbs phenomenon are eliminated. As a result of the proposed FC-based domain decomposition approach, eigenvalue problems for complicated geometries and general smooth coefficients can be tackled with high accuracy within the FC framework. Additionally, the proposed domain decomposition approach can be used to produce local mesh refinements as well as efficient parallel implementations.The capabilities afforded by the proposed methodologies are demonstrated by means of a variety of examples in two-dimensional space. Extensions of this framework to higher dimensions, which are straightforward, are not otherwise pursued in this paper.


## 2. Definitions

### 2.1. Overall Methodology

Our general approach relies on an decomposition of the domain into overlapping patches that are expressed in terms of curvilinear coordinates (i.e., mapped to rectangles with uniform grids), using parametrizations to describe the differential operator in terms of Cartesian reference coordinates, which in what follows are denoted bys and t. On each patch, boundary conditions and continuity are imposed, and then the operator is evaluated using Fast Fourier Transforms. The numerical differential operator thus defined on the whole domain can then be used in conjunction with an iterative method for eigenvalue evaluation. The basic concepts and notations used throughout this process are discussed in what follows.


### 2.2. Eigenvalue Problem

LetΩ
⊂
R
2 be a bounded domain and let L be a linear elliptic operator of the form L
=
∑
|
α
|
≤
m
‍
a
α
(
x
,
y
)
∂
α, where α is a two-dimensional multi-index. Throughout this paper, we search for a subset of the scalar values λ and corresponding functions u
∈
C
∞
(
Ω
,
R
) such that(2)
L
u
=
λ
u
,
u
∂
Ω
1
=
0
,
∂
u
∂
n
^
∂
Ω
2
=
0
,where λ are the eigenvalues and u
(
x
,
y
) are the eigenfunctions. The boundary conditions are defined so that ∂
Ω
1
∪
∂
Ω
2 is a disjoint union that partitions ∂
Ω and n
^ denotes the outward unit normal. For simplicity we restrict the scope of our study to only consider fixed (Dirichlet) and free (Neumann) boundary conditions.


### 2.3. Domain Decomposition

LetP
i
i
=
1
n be open sets relative to Ω such that there exist invertible smooth functions f
i
:
R
→
P
i, where R
=
[
0,1
]
×
[
0,1
]. We let (
s
,
t
) denote coordinates on R, which correspond to (
x
,
y
) on P
i under the mapping f
i. Using the chain rule to translate derivatives in domain coordinates to derivatives in reference coordinates (e.g., u
s
=
x
s
u
x
+
y
s
u
y) and inverting the resulting system yields the desired expressions for derivatives in x and y in terms of s and t. The derivatives in s and t are calculated numerically on the uniform reference grid using the FC(Gram) algorithm described in the following section. For the examples considered in this paper, the domain decompositions were achieved by inspection of the domain and trial-and-error manipulations, but a variety of algorithms are available that perform overlapping patch decompositions (e.g., [15]).


### 2.4. FC(Gram)

Evaluation of derivatives in the Fourier basis offers two main advantages: under the Fourier Transform, differentiation corresponds to scalar multiplication, and, in view of the Fast Fourier Transform algorithm, numerical evaluation of the Fourier Transform scales well to large problems, since the cost of a Fourier Transform onn modes is O
(
n
log
⁡
n
). As is known, however, truncated Fourier series provide poor representations for nonperiodic functions, an inaccuracy which is only compounded by differentiation. Specifically, the Fourier series oscillates rapidly at the endpoints of the domain whenever the periodic extension is discontinuous (Gibbs phenomenon). Since general functions to which the differential operator must be applied are not periodic, high orders of accuracy cannot be expected from straightforward uses of Fourier expansions.This problem may be addressed by means of the Fourier continuation algorithm, which, given a functionf, produces Fourier series for smooth periodic extensions of f to a suitably larger interval. Thus, given a function f defined on (say) [
0,1
] which is sampled on the discrete grid given by x
j
=
(
j
-
1
)
h, where h
=
1
/
(
n
-
1
) and j
=
1
,
…
,
n, the algorithm constructs a new function f
c on [
0
,
b
] for some b
>
1 such that f
c
(
x
)
=
f
(
x
),  ∀
x
∈
[
0,1
], and f
c
m
b
=
f
c
m
0, ∀
m
≤
k. In view of these conditions the function f
c is smooth and periodic, and therefore its b-periodic trigonometric polynomial,(3)
f
c
x
=
∑
k
∈
t
F
a
k
e
2
π
i
/
b
k
x
,(where t
(
F
) is the set of n
+
n
c modes distributed symmetrically about 0) is highly accurate. To enforce the approximation of the function f, one must compute the coefficients a
k in (3) such that f
c and f agree on the n original sample points. This could, for example, be accomplished by solving a (possibly oversampled) minimization problem:(4)
min
a
k
∣
k
∈
t
F
⁡
∑
j
=
1
n
∑
k
∈
t
F
a
k
e
2
π
i
/
b
k
x
j
-
f
x
j
2
.Written as a matrix equation, this amounts to solving(5)
M
a
=
f
,in the least-squares sense, where a is a column vector whose entries are the coefficients a
k, f is a column vector whose entries are f
(
x
j
), and M is a matrix with M
j
k
=
e
(
2
π
i
/
b
)
k
x
j.By itself, this procedure constructs the desired extension but is not desirable for high-efficiency computation, which is needed in the context of an iterative solver. This method can be modified, however, in such a way that only a constant number of data points are used in the calculation of the continuation for arbitrarily large values ofn. A brief summary of the resulting approach is presented in what follows; full details can be found in [1, 16]. Consider the set of points {
(
x
n
-
n
l
-
1
,
f
(
x
n
-
n
l
-
1
)
)
,
…
,
(
x
n
,
f
(
x
n
)
)
}
∪
{
(
x
1
+
1
+
d
,
f
(
x
1
)
)
,
…
,
(
x
n
r
+
1
+
d
,
f
(
x
n
r
)
)
}, where n
l and n
r are fixed integers chosen independently of n. By computing the trigonometric polynomial fit via (4) on this set of points, we can construct a function f
m
a
t
c
h such that(6)
f
c
x
=
f
x
,
x
∈
0,1
,
f
match
x
,
x
∈
1,1
+
dis (
1
+
d
)-periodic. A Fourier approximation of the smooth and periodic function f
c in the interval [
0,1
+
d
] (which also provides an approximation for f in the interval [
0,1
]) is highly accurate, and it provides, in particular, the needed Fourier expansion of the function f in the interval [
0,1
].It is possible to selectn
l and n
r small while maintaining high accuracy in the continuation by expressing the aforementioned least-squares problem in terms of the orthonormal (Gram) polynomial basis for the intervals [
x
n
-
n
l
-
1
,
1
] and [
1
+
d
,
x
n
r
+
1
+
d
]. The algorithm is completed by relying on highly accurate precomputed extensions for pairs of Gram polynomials on the left and right subintervals, projection of a given function onto the Gram basis, and a subsequent use of the Fast Fourier Transform; see [1, 16] for details. An example of an extension computed via FC(Gram) is shown for a quadratic function on a closed interval in Figure 1. All FC(Gram) continuations in our numerical examples use n
l
=
n
r
=
10 and a Gram polynomial basis up to degree 5.Figure 1
The functionf
(
x
)
=
x
2 on [
0,1
] extended via the FC(Gram) algorithm to f
c
(
x
). The function has been extended so that it is smoothly periodic over a slightly larger domain.


### 2.5. Numerical Differentiation

Iff
(
x
) is a 1-dimensional function whose continuation (computed via FC(Gram)) is f
c
(
x
), then the derivative of the Fourier series of f
c is computed as follows:(7)
f
c
x
=
∑
l
=
-
∞
∞
f
l
e
i
ω
l
x
,
d
k
f
c
d
x
=
∑
l
=
-
∞
∞
i
ω
l
k
f
l
e
i
ω
l
x
.The derivative of f is recovered by restricting the domain of f
c
(
k
). This method easily generalizes to multivariate functions defined on rectangular domains. In the discrete setting, we compute Fourier Transforms using the FFTW implementation [17] of the Fast Fourier Transform. The one-dimensional FC(Gram) differentiation scheme is summarized in Algorithm 1. Partial derivatives in two dimensions are computed numerically by applying the one-dimensional algorithm to each row (or column) of the equispaced grid in reference coordinates. This method is indeed accurate to high order, as shown in previous applications to ODE and PDE problems [1, 16, 18]. Since we use a 5th-order Gram basis to construct the fits, we expect Algorithm 1 to compute derivatives with fourth-order accuracy.Algorithm 1:FC(Gram) 1D kth-order differentiation.
(1)    Procedure  FCDIFF(u
,
k
,
h)            ▹ returns kth derivative of u, where h is the mesh step size
(2)     u
c
←
c
o
n
t(u)                       ▹ perform FC(Gram) continuation
(3)     u
c
^
←
f
f
t(u
c)
(4)     for  j
←
1
,
n
c  do                           ▹  n
c is the length of u
c
(5)      u
c
^
[
j
]
←
u
c
^
[
j
]
∗
(
2
π
i
/
(
n
c
h
)
)
k
(6)     end for
(7)     u
c
←
i
f
f
t(u
c
^)
(8)     u
←
r
e
s
t
r
i
c
t(u
c)                   ▹ store the portion corresponding to the original domain
(9)     return  u
(10) end procedure


## 2.1. Overall Methodology

Our general approach relies on an decomposition of the domain into overlapping patches that are expressed in terms of curvilinear coordinates (i.e., mapped to rectangles with uniform grids), using parametrizations to describe the differential operator in terms of Cartesian reference coordinates, which in what follows are denoted bys and t. On each patch, boundary conditions and continuity are imposed, and then the operator is evaluated using Fast Fourier Transforms. The numerical differential operator thus defined on the whole domain can then be used in conjunction with an iterative method for eigenvalue evaluation. The basic concepts and notations used throughout this process are discussed in what follows.


## 2.2. Eigenvalue Problem

LetΩ
⊂
R
2 be a bounded domain and let L be a linear elliptic operator of the form L
=
∑
|
α
|
≤
m
‍
a
α
(
x
,
y
)
∂
α, where α is a two-dimensional multi-index. Throughout this paper, we search for a subset of the scalar values λ and corresponding functions u
∈
C
∞
(
Ω
,
R
) such that(2)
L
u
=
λ
u
,
u
∂
Ω
1
=
0
,
∂
u
∂
n
^
∂
Ω
2
=
0
,where λ are the eigenvalues and u
(
x
,
y
) are the eigenfunctions. The boundary conditions are defined so that ∂
Ω
1
∪
∂
Ω
2 is a disjoint union that partitions ∂
Ω and n
^ denotes the outward unit normal. For simplicity we restrict the scope of our study to only consider fixed (Dirichlet) and free (Neumann) boundary conditions.


## 2.3. Domain Decomposition

LetP
i
i
=
1
n be open sets relative to Ω such that there exist invertible smooth functions f
i
:
R
→
P
i, where R
=
[
0,1
]
×
[
0,1
]. We let (
s
,
t
) denote coordinates on R, which correspond to (
x
,
y
) on P
i under the mapping f
i. Using the chain rule to translate derivatives in domain coordinates to derivatives in reference coordinates (e.g., u
s
=
x
s
u
x
+
y
s
u
y) and inverting the resulting system yields the desired expressions for derivatives in x and y in terms of s and t. The derivatives in s and t are calculated numerically on the uniform reference grid using the FC(Gram) algorithm described in the following section. For the examples considered in this paper, the domain decompositions were achieved by inspection of the domain and trial-and-error manipulations, but a variety of algorithms are available that perform overlapping patch decompositions (e.g., [15]).


## 2.4. FC(Gram)

Evaluation of derivatives in the Fourier basis offers two main advantages: under the Fourier Transform, differentiation corresponds to scalar multiplication, and, in view of the Fast Fourier Transform algorithm, numerical evaluation of the Fourier Transform scales well to large problems, since the cost of a Fourier Transform onn modes is O
(
n
log
⁡
n
). As is known, however, truncated Fourier series provide poor representations for nonperiodic functions, an inaccuracy which is only compounded by differentiation. Specifically, the Fourier series oscillates rapidly at the endpoints of the domain whenever the periodic extension is discontinuous (Gibbs phenomenon). Since general functions to which the differential operator must be applied are not periodic, high orders of accuracy cannot be expected from straightforward uses of Fourier expansions.This problem may be addressed by means of the Fourier continuation algorithm, which, given a functionf, produces Fourier series for smooth periodic extensions of f to a suitably larger interval. Thus, given a function f defined on (say) [
0,1
] which is sampled on the discrete grid given by x
j
=
(
j
-
1
)
h, where h
=
1
/
(
n
-
1
) and j
=
1
,
…
,
n, the algorithm constructs a new function f
c on [
0
,
b
] for some b
>
1 such that f
c
(
x
)
=
f
(
x
),  ∀
x
∈
[
0,1
], and f
c
m
b
=
f
c
m
0, ∀
m
≤
k. In view of these conditions the function f
c is smooth and periodic, and therefore its b-periodic trigonometric polynomial,(3)
f
c
x
=
∑
k
∈
t
F
a
k
e
2
π
i
/
b
k
x
,(where t
(
F
) is the set of n
+
n
c modes distributed symmetrically about 0) is highly accurate. To enforce the approximation of the function f, one must compute the coefficients a
k in (3) such that f
c and f agree on the n original sample points. This could, for example, be accomplished by solving a (possibly oversampled) minimization problem:(4)
min
a
k
∣
k
∈
t
F
⁡
∑
j
=
1
n
∑
k
∈
t
F
a
k
e
2
π
i
/
b
k
x
j
-
f
x
j
2
.Written as a matrix equation, this amounts to solving(5)
M
a
=
f
,in the least-squares sense, where a is a column vector whose entries are the coefficients a
k, f is a column vector whose entries are f
(
x
j
), and M is a matrix with M
j
k
=
e
(
2
π
i
/
b
)
k
x
j.By itself, this procedure constructs the desired extension but is not desirable for high-efficiency computation, which is needed in the context of an iterative solver. This method can be modified, however, in such a way that only a constant number of data points are used in the calculation of the continuation for arbitrarily large values ofn. A brief summary of the resulting approach is presented in what follows; full details can be found in [1, 16]. Consider the set of points {
(
x
n
-
n
l
-
1
,
f
(
x
n
-
n
l
-
1
)
)
,
…
,
(
x
n
,
f
(
x
n
)
)
}
∪
{
(
x
1
+
1
+
d
,
f
(
x
1
)
)
,
…
,
(
x
n
r
+
1
+
d
,
f
(
x
n
r
)
)
}, where n
l and n
r are fixed integers chosen independently of n. By computing the trigonometric polynomial fit via (4) on this set of points, we can construct a function f
m
a
t
c
h such that(6)
f
c
x
=
f
x
,
x
∈
0,1
,
f
match
x
,
x
∈
1,1
+
dis (
1
+
d
)-periodic. A Fourier approximation of the smooth and periodic function f
c in the interval [
0,1
+
d
] (which also provides an approximation for f in the interval [
0,1
]) is highly accurate, and it provides, in particular, the needed Fourier expansion of the function f in the interval [
0,1
].It is possible to selectn
l and n
r small while maintaining high accuracy in the continuation by expressing the aforementioned least-squares problem in terms of the orthonormal (Gram) polynomial basis for the intervals [
x
n
-
n
l
-
1
,
1
] and [
1
+
d
,
x
n
r
+
1
+
d
]. The algorithm is completed by relying on highly accurate precomputed extensions for pairs of Gram polynomials on the left and right subintervals, projection of a given function onto the Gram basis, and a subsequent use of the Fast Fourier Transform; see [1, 16] for details. An example of an extension computed via FC(Gram) is shown for a quadratic function on a closed interval in Figure 1. All FC(Gram) continuations in our numerical examples use n
l
=
n
r
=
10 and a Gram polynomial basis up to degree 5.Figure 1
The functionf
(
x
)
=
x
2 on [
0,1
] extended via the FC(Gram) algorithm to f
c
(
x
). The function has been extended so that it is smoothly periodic over a slightly larger domain.


## 2.5. Numerical Differentiation

Iff
(
x
) is a 1-dimensional function whose continuation (computed via FC(Gram)) is f
c
(
x
), then the derivative of the Fourier series of f
c is computed as follows:(7)
f
c
x
=
∑
l
=
-
∞
∞
f
l
e
i
ω
l
x
,
d
k
f
c
d
x
=
∑
l
=
-
∞
∞
i
ω
l
k
f
l
e
i
ω
l
x
.The derivative of f is recovered by restricting the domain of f
c
(
k
). This method easily generalizes to multivariate functions defined on rectangular domains. In the discrete setting, we compute Fourier Transforms using the FFTW implementation [17] of the Fast Fourier Transform. The one-dimensional FC(Gram) differentiation scheme is summarized in Algorithm 1. Partial derivatives in two dimensions are computed numerically by applying the one-dimensional algorithm to each row (or column) of the equispaced grid in reference coordinates. This method is indeed accurate to high order, as shown in previous applications to ODE and PDE problems [1, 16, 18]. Since we use a 5th-order Gram basis to construct the fits, we expect Algorithm 1 to compute derivatives with fourth-order accuracy.Algorithm 1:FC(Gram) 1D kth-order differentiation.
(1)    Procedure  FCDIFF(u
,
k
,
h)            ▹ returns kth derivative of u, where h is the mesh step size
(2)     u
c
←
c
o
n
t(u)                       ▹ perform FC(Gram) continuation
(3)     u
c
^
←
f
f
t(u
c)
(4)     for  j
←
1
,
n
c  do                           ▹  n
c is the length of u
c
(5)      u
c
^
[
j
]
←
u
c
^
[
j
]
∗
(
2
π
i
/
(
n
c
h
)
)
k
(6)     end for
(7)     u
c
←
i
f
f
t(u
c
^)
(8)     u
←
r
e
s
t
r
i
c
t(u
c)                   ▹ store the portion corresponding to the original domain
(9)     return  u
(10) end procedure


## 3. Solving the Eigenvalue Problem

### 3.1. Boundary Conditions

Enforcing a Dirichlet boundary condition is as simple as prescribing the value on a line of fringe points along the boundary. Due to the mappingf
i, the boundaries of R correspond exactly to the boundaries of P
i, so Dirichlet boundary conditions can be easily translated to reference coordinates. However, imposing a Neumann boundary condition is more complicated due to the fact that normal vectors are not held invariant by f
i. To translate the condition to reference coordinates, we again use the chain rule to write the normal derivative as a combination of change of variables terms and derivatives in parameter space (by inverting the Jacobian of the mapping f
i):(8)
0
=
∂
u
∂
n
^
=
∇
u
·
n
^
=
u
x
n
1
+
u
y
n
2
=
x
s
y
t
-
y
s
x
t
-
1
u
s
n
1
y
t
-
n
2
x
t
+
u
t
-
n
1
y
s
+
n
2
x
s
≔
γ
1
u
s
+
γ
2
u
t
.The values γ
1 and γ
2 depend only on the unit outward normal n
^
=
(
n
1
,
n
2
)
T and the partial derivatives of components of f
i (for a given patch P
i).By treating the points on the boundary, where a Neumann condition is prescribed as unknowns, we will already have data from the candidate functionu along the boundary. As a result, either u
s or u
t can be numerically computed along the boundary, which in turn allows us to describe the Neumann condition solely in terms of the other partial (i.e., by solving for it in (8)). This means that the condition on the normal derivative in domain space can be written as a condition on the normal derivative in reference coordinates. Knowing this, we use a least-squares system to fit a fifth-order 1D polynomial to each set of 10 grid points along the normal direction in parameter space. This polynomial is constructed with the condition that its derivative has the value derived from (8). Finally, the values at these grid points are replaced by the values from the polynomial fit when computing the Fourier continuation, but the original values are used for the calculation of the derivative (after continuation). This way, if a given function already obeys the Neumann condition, nothing is changed and the operator is computed properly. However, if the condition does not hold, this process introduces the desired discontinuity in the function which increases in magnitude when derivatives are computed.


### 3.2. Continuity

A similar method is used to enforce continuity across patches in overlap regions. SupposeP
i overlaps with P
j near one of its boundaries. A layer of fringe points is added to the patch on the overlap region. Then the function values for P
i on the fringe points are prescribed by interpolating from the function values in the interior of P
j (avoiding values in P
j that are themselves interpolated from other patches). We use Neville’s algorithm for cubic polynomial interpolation and perform all interpolation in reference coordinates. A simple implementation of Newton’s method is used to compute f
j
-
1
(
x
,
y
) (when f
j
-
1 is not explicitly known) for a given point (
x
,
y
) at the boundary of P
i.


### 3.3. Eigenvalue Algorithm

The methodologies described above in this paper give rise to a numerical algorithm for evaluation of differential operators for a given geometry with a provided domain decomposition and boundary conditions. To compute the eigenvalues in (2), we combine this numerical differential operator with the ARPACK++ [19] implementation of the implicitly restarted Arnoldi method. The Arnoldi method iteratively builds a lower dimensional approximation of the operator L as an upper Hessenberg matrix H. The eigenvalues of H (“Ritz values”) are, in practice, good estimates of the eigenvalues of L, as are the eigenvectors (up to an orthogonal transformation). The special structure of H allows for efficient calculation of the Ritz values as well. Furthermore, we find that the implicitly restarted Arnoldi method performs better than the full Arnoldi method (i.e., without restarts); the restarted method is dominated mostly by the cost of evaluating the operator L rather than the cost of computing the eigenvalues of a large H matrix. Combined with acceleration techniques that damp out components of unwanted eigenvectors at each restart, the former method is faster in practice, so ARPACK++ provides the necessary computational framework for computing eigenvalues of L using our FC algorithm. All eigenvalues in the following numerical experiments are computed to a tolerance of 1
0
-
7 unless otherwise stated.


## 3.1. Boundary Conditions

Enforcing a Dirichlet boundary condition is as simple as prescribing the value on a line of fringe points along the boundary. Due to the mappingf
i, the boundaries of R correspond exactly to the boundaries of P
i, so Dirichlet boundary conditions can be easily translated to reference coordinates. However, imposing a Neumann boundary condition is more complicated due to the fact that normal vectors are not held invariant by f
i. To translate the condition to reference coordinates, we again use the chain rule to write the normal derivative as a combination of change of variables terms and derivatives in parameter space (by inverting the Jacobian of the mapping f
i):(8)
0
=
∂
u
∂
n
^
=
∇
u
·
n
^
=
u
x
n
1
+
u
y
n
2
=
x
s
y
t
-
y
s
x
t
-
1
u
s
n
1
y
t
-
n
2
x
t
+
u
t
-
n
1
y
s
+
n
2
x
s
≔
γ
1
u
s
+
γ
2
u
t
.The values γ
1 and γ
2 depend only on the unit outward normal n
^
=
(
n
1
,
n
2
)
T and the partial derivatives of components of f
i (for a given patch P
i).By treating the points on the boundary, where a Neumann condition is prescribed as unknowns, we will already have data from the candidate functionu along the boundary. As a result, either u
s or u
t can be numerically computed along the boundary, which in turn allows us to describe the Neumann condition solely in terms of the other partial (i.e., by solving for it in (8)). This means that the condition on the normal derivative in domain space can be written as a condition on the normal derivative in reference coordinates. Knowing this, we use a least-squares system to fit a fifth-order 1D polynomial to each set of 10 grid points along the normal direction in parameter space. This polynomial is constructed with the condition that its derivative has the value derived from (8). Finally, the values at these grid points are replaced by the values from the polynomial fit when computing the Fourier continuation, but the original values are used for the calculation of the derivative (after continuation). This way, if a given function already obeys the Neumann condition, nothing is changed and the operator is computed properly. However, if the condition does not hold, this process introduces the desired discontinuity in the function which increases in magnitude when derivatives are computed.


## 3.2. Continuity

A similar method is used to enforce continuity across patches in overlap regions. SupposeP
i overlaps with P
j near one of its boundaries. A layer of fringe points is added to the patch on the overlap region. Then the function values for P
i on the fringe points are prescribed by interpolating from the function values in the interior of P
j (avoiding values in P
j that are themselves interpolated from other patches). We use Neville’s algorithm for cubic polynomial interpolation and perform all interpolation in reference coordinates. A simple implementation of Newton’s method is used to compute f
j
-
1
(
x
,
y
) (when f
j
-
1 is not explicitly known) for a given point (
x
,
y
) at the boundary of P
i.


## 3.3. Eigenvalue Algorithm

The methodologies described above in this paper give rise to a numerical algorithm for evaluation of differential operators for a given geometry with a provided domain decomposition and boundary conditions. To compute the eigenvalues in (2), we combine this numerical differential operator with the ARPACK++ [19] implementation of the implicitly restarted Arnoldi method. The Arnoldi method iteratively builds a lower dimensional approximation of the operator L as an upper Hessenberg matrix H. The eigenvalues of H (“Ritz values”) are, in practice, good estimates of the eigenvalues of L, as are the eigenvectors (up to an orthogonal transformation). The special structure of H allows for efficient calculation of the Ritz values as well. Furthermore, we find that the implicitly restarted Arnoldi method performs better than the full Arnoldi method (i.e., without restarts); the restarted method is dominated mostly by the cost of evaluating the operator L rather than the cost of computing the eigenvalues of a large H matrix. Combined with acceleration techniques that damp out components of unwanted eigenvectors at each restart, the former method is faster in practice, so ARPACK++ provides the necessary computational framework for computing eigenvalues of L using our FC algorithm. All eigenvalues in the following numerical experiments are computed to a tolerance of 1
0
-
7 unless otherwise stated.


## 4. Numerical Results

### 4.1. Simple Geometries

To test the efficiency and accuracy of the eigenvalue solver as well as the domain decomposition approach we first present numerical results for two simple geometries, the unit square and the unit disk, usingL
=
∇
2. For each of these geometries, we consider both Dirichlet and Neumann homogeneous boundary conditions. The exact eigenvalues and eigenfunctions are well known and easily derived. We report them here for reference:(1)
Unit square (Ω
=
[
0,1
]
×
[
0,1
]) with all sides fixed: λ
m
n
=
π
2
(
m
2
+
n
2
) for m
,
n
∈
Z
+.
(2)
Unit square with one side (x
=
0) free: λ
m
n
=
π
2
(
(
m
+
1
/
2
)
2
+
n
2
) for m
∈
Z
+
∪
{
0
} and n
∈
Z
+.
(3)
Unit disk (Ω
=
{
(
x
,
y
)
∈
R
2
∣
x
2
+
y
2
≤
1
}) with fixed boundary: λ
m
n
=
j
n
m
2, where j
n
m is the nth root of the mth Bessel function.
(4)
Unit disk with free boundary:λ
m
n
=
j
′
2
n
m, where j
′
n
m is the nth root of the derivative of the mth Bessel function.For these cases, the domain decompositions we use are not optimal (in the sense of minimizing the number of patches), as they are meant to test how well continuity conditions are enforced. In addition, we intentionally discretize neighboring patches with slightly different numbers of grid points (to prevent grid lines from overlapping and artificially increasing the accuracy of interpolation). The decompositions used for the square and disk cases are shown in Figure2.Figure 2
Decomposition of simple geometries: (a) the unit square, decomposed into two overlapping rectangles; (b) the unit disk, decomposed into eight overlapping sections of an annulus and a square centered at the origin. These examples demonstrate the usefulness of domain decomposition to avoid difficulties in mapping entire domains.
(a)
(b)The error in the first 10 computed eigenvalues for all the cases and for increasingly finer discretizations is shown in Figure3. In each case, we observe fourth-order convergence (e.g., doubling the number of points along one direction of each patch decreases the error in the eigenvalues by no less than 2
-
4). In addition, for the finest discretization in each case, we obtain at least 8 digits of accuracy in the principal eigenvalue. Finally, the number of calls to the numerical operator (denoted by L
u) and the time taken per case are reported in Table 1.Table 1
Computational time required for evaluation of the differential operatorL
u for four test cases, with breakdown of computing time per L
u calculation, and the number of L
u calculations needed for the eigenvalue solver to converge to the required tolerance. All computations were performed on a 2.67 GHz Intel Xeon processor. (a) Square with fixed edges; (b) square with one free edge; (c) disk with fixed boundary; (d) disk with free boundary.
(a)
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
1525
0.001691
728
1.2313
6100
0.006096
2179
13.2835
24400
0.020961
6003
125.8296
97600
0.103420
22861
2364.2878
(b)
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
1525
0.001891
923
1.7452
6100
0.005459
2702
14.7507
24400
0.019813
8328
165.0070
97600
0.094591
23440
2217.2179
(c)
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
3950
0.006348
1065
6.7605
15800
0.017500
2037
35.6459
63200
0.072966
4981
363.4445
252800
0.257892
12338
3181.8739
(d)
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
3950
0.006638
1584
10.5144
15800
0.012105
3336
40.3828
63200
0.061148
7027
429.6900
252800
0.222666
14820
3299.9156Figure 3
Convergence in the first ten eigenvalues for each of the four cases in Section4.1: (a) unit square and fixed edges; (b) unit square and one free edge; (c) unit disk and fixed boundary; (d) unit disk and free boundary (with zero eigenvalues omitted). In all cases, the dashed line for fourth-order convergence is shown for comparison. N denotes the number of points used to discretize the overall two-dimensional geometry in each case.
(a)
(b)
(c)
(d)


### 4.2. Comparison to Finite Differences

Because the Arnoldi iteration is a black box algorithm that depends only on the definition of the numerical operator, one might consider a very similar algorithm that replaces the FC(Gram) algorithm with another method for computing derivatives. In this section, we use the examples of the previous section to show that the choice of FC(Gram) is justified when compared to a standard finite difference scheme and that the advantages of FC(Gram) are quite apparent. (It is worth noting the recent contribution in [20], which relies on finite differences and Richardson extrapolation, can provide highly accurate approximations of Laplace eigenvalues provided that the domain can be discretized by means of grids which exactly sample the boundary of the domain.)In these tests, we compare the number of iterations, time taken, and number of grid points needed by the two methods to compute the principal eigenvalue of the Laplacian to a fixed amount of error. We consider both the unit disk and the unit square with Dirichlet boundary conditions, and a second-order accurate centered finite difference scheme is used in the comparison. As the results in Table2 indicate, the FC(Gram) algorithm requires significantly fewer iterations and grid points than the finite difference algorithm to reach the desired accuracy. This also results in a much lower overall run time, even though the time per operation is smaller for the finite difference scheme. In addition, the FC(Gram) algorithm scales well to higher accuracies, while the cost of the finite difference algorithm grows prohibitively large for the finer discretizations. These results are even more pronounced in the case of the unit disk. Finally, simple tests using higher order Gram bases for the FC(Gram) algorithm demonstrate that the order of accuracy can be increased without much change in the computational cost. At the same time, there is no significant alteration to the FC(Gram) algorithm at high orders, while high-order finite difference schemes must be specially reformulated at boundaries to take into account the extended stencils. Thus, FC(Gram) seems to be well suited to the elliptic eigenvalue problems we consider here and competitive with other methods.Table 2
Number of iterations and computational time required by the FC solver and the finite difference solver to computeλ
1 to a fixed accuracy for the unit square example (fixed boundaries). (The first 10 eigenvalues are calculated to the tolerance of 1
0
-
7, but we only show the error in the first eigenvalue for conciseness.) For each given accuracy, the FC(Gram) algorithm outperforms the finite difference algorithm by a large margin. All computations were performed on a 2.67 GHz Intel Xeon processor. (a) Results using centered finite difference; (b) results using FC(Gram) with 5th-order basis.
(a)
Absolute error
N
Time perL
u operation (s)
Number ofL
u operations
3.00 × 10−4
74725
0.009501
9693
7.57 × 10−5
298800
0.05130
37760
(b)
Absolute error
N
Time perL
u operation (s)
Number ofL
u operations
3.53 × 10−4
1525
0.001235
728
1.41 × 10−5
6100
0.004211
2245


### 4.3. Smoothed L-Shaped Membrane

Collocation methods were used to calculate the eigenvalues of the L-shaped membrane with Dirichlet boundary conditions in order to show how polygonal boundaries, including those with geometric singularities, could be handled. Improvements to Fox et al.’s work [4] on the method of particular solutions have yielded the first several eigenvalues to at least 10 digits of accuracy [5]. In its present form the FC method does not handle geometries containing corner singularities, but it can be applied to approximations thereof containing sharp but smooth corners. (In the presence of corners in the boundary of the domain, PDE solutions and eigenfunctions are singular and therefore cannot be approximated with high accuracy by means of the present version of the FC method. An extension of this method could be envisioned which takes into account such singularities explicitly, but the development of such approaches lies beyond the scope of this paper.) We thus consider the eigenvalue problem for an L-shaped membrane with a smoothed corner, as discussed in what follows.To decompose the membrane, we smooth the corner by locally approximating the boundary curve withC
∞ curve that is sampled at a fixed set of points. Evaluation at an arbitrary location is produced by means of Neville’s algorithm for cubic interpolation on a sufficiently fine mesh to ensure that approximation errors are below the tolerances otherwise required. The smooth curve that replaces the corner is given (up to rigid transformations) by(9)
g
x
=
∫
-
∞
x
∫
-
∞
y
ϕ
z
d
z
d
y
,where the function ϕ
(
z
) is a smooth bump function on R with support on [
-
1,1
]:(10)
ϕ
z
=
0
,
if
z
∉
-
1,1
,
e
-
1
/
1
-
z
e
-
1
/
1
-
z
+
e
-
1
/
z
,
if
z
∈
-
1,1
.Thus, outside the interval [
-
1,1
], g
(
x
)
=
|
x
| and matches the rest of the boundary of the L-shaped membrane near the corner. Within it, we have a smooth curve that deviates from the corner, as shown in Figure 4. (This can be thought of as a way of approximating the Dirac delta function, which is the second “derivative” of the absolute value function, with a smooth function that has a peak of finite height and nonzero width.)Figure 4
The smoothed corner coincides exactly with the absolute value function outside of[
-
1,1
].Integral (9) is obtained numerically using Fourier techniques that take advantage of the fact that ϕ
(
z
) is C
∞ periodic function on [
-
1,1
]. We sample ϕ at 1024 equally spaced points on [
-
1,1
] and compute the FFT of the discrete data; at higher frequencies, the Fourier coefficients are effectively zero. Integration in the frequency domain amounts to multiplying each coefficient by a constant. Thus, taking the inverse FFT after multiplication yields the integral of ϕ on [
-
1,1
] up to a constant of integration. The second integral is computed the same way, with the constant term integrated analytically. This yields a highly accurate representation of g
(
x
) and its derivatives, which are necessary for the calculation of the Laplacian. Note that this method can be generalized for smoothing reentrant corners of any angle by choosing the first integration constant or adjusting the height of the peak in ϕ accordingly.Using this description of the approximate boundary curve we parametrize the region of the L-shaped membrane near the corner by(11)
x
~
s
,
t
=
2
+
2
l
s
-
l
+
1
,
y
~
s
,
t
=
t
b
g
2
+
2
l
s
-
l
+
1
,
x
y
=
d
cos
⁡
θ
-
sin
⁡
θ
sin
⁡
θ
cos
⁡
θ
x
~
y
~
+
x
0
y
0
+
h
sin
⁡
θ
cos
⁡
θ
,for s
∈
[
0,1
] and t
∈
[
0
,
b
]. Here θ
=
-
π
/
4, (
x
0
,
y
0
)
=
(
1,1
), and l and h are chosen to be of the same order as d. The decomposition of the smoothed L-shaped membrane is shown in Figure 5. This particular parametrization is chosen so that as d
→
0 the curve approximates the corner of the L-shaped membrane arbitrarily well.Figure 5
Domain decomposition for the smoothed L-shaped membrane. The corner patch (d
=
0.05) is drawn using a dashed line for emphasis.Figure7 indicates that 4th-order convergence is achieved, just as in the tests for the simpler unit square and unit disk geometries; run time information is given in Table 3. Of equal importance is the fact that the computed eigenfunctions (Figure 6) closely resemble those of the true L-shaped membrane [2] and the corresponding eigenvalues agree to anywhere from two to four digits, as one might expect. In practice, we have found that the eigenvalues converge linearly as the rounding radius of curvature tends to zero.Table 3
Computational time required for evaluation of the Laplacian for the smoothed L-shaped membrane. All computations were performed on a 2.67 GHz Intel Xeon processor.
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
9216
0.012912
6820
88.0598
18225
0.018062
14736
266.1676
36864
0.031643
21469
679.3593
72990
0.049035
44212
2167.9487
147456
0.149842
88070
13196.5794Figure 6
The first 5 eigenfunctions (modes of vibration) of the smoothed L-shaped membrane (cf. Figure  1 in [2]): (a) λ
1
=
-
9.6013; (b) λ
2
=
-
15.1969; (c) λ
3
=
-
19.7392; (d) λ
4
=
-
29.5214; (e) λ
5
=
-
31.8179.
(a)
(b)
(c)
(d)
(e)Figure 7
Convergence in the first 5 eigenvalues for the smoothed L-shaped membrane (d
=
0.05) is roughly 4th order when compared to the finest discretization tested (N
=
147456). Note that the eigenvalues do not necessarily converge to those of the true L-shaped membrane with a sharp corner due to the approximation of the corner geometry. For example, here λ
1
=
9.6013 when N
=
147456, even though λ
1
=
9.63972 (e.g., [2, 4]) for the membrane with the sharp corner.


### 4.4. Thinly Tethered Membrane

In this section we explore a more physically relevant membrane geometry inspired by quantum optics experiments. The particular geometry explored here can be described as a roughly square-shaped membrane (with side length on the order of 100μm) attached by thin tethers coming off of the four corners to a square frame (which has side length on the order of 600 μm). The edges of the membrane that are not attached to the frame are left free. In addition, the tethers smoothly continue into the boundaries of the inner square and also taper smoothly at the corners of the frame. (The purpose of this geometry is to decrease stress on the main part of the membrane, which would theoretically drive down the fundamental frequency.) For this experiment, we take the stress across the membrane to be constant in order to obtain a first approximation to the mechanics of such a complex membrane. Under this approximation, the operator in (2) is once again the Laplace operator (up to a constant factor).The decomposition developed for this membrane is shown in Figure8. Nearly all of the patches used are rectangles, but the smooth boundary curves at the corner of the membrane and the frame provide a challenge. To address this issue, we use the methodology previously developed for smoothing corners to parametrize these regions; we treat the segments of the boundary shown in Figures 8(b) and 8(c) as having smoothed reentrant corners with angles of 90° and 135°, respectively. The dimensions of the membrane are scaled so that the frame is the square of side length 2 centered at the origin.Figure 8
Domain decomposition for the thinly tethered membrane: (a) full decomposition (37 patches); (b) curving boundary near point where tether attaches to center square; (c) curving boundary near corner of frame.
(a)
(b)
(c)As in all the previous cases, we consider several successively finer discretizations for the geometry. Because this problem is significantly larger than those considered previously in this paper, we only increase the number of points in each direction by a factor of 1.2 at each level. Convergence rates for the first 5 eigenvalues are shown in Figure9 and timings are given in Table 4. In addition, plots of the eigenfunctions corresponding to these computed eigenvalues are shown in Figure 10, which show interesting oscillation modes.Table 4
Computational times and number of iterations required for evaluation of the Laplacian on the thinly tethered membrane. All computations were performed on a 2.67 GHz Intel Xeon processor to a tolerance of1
0
-
4.
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
56320
0.061338
464727
28505.4983
79420
0.072279
1758160
127079.4402
116380
0.115435
2679346
309291.7615
166045
0.167378
2793136
467510.6446Figure 9
Convergence in the first 5 eigenvalues for the thinly tethered membrane is roughly 4th order (when compared to the finest discretization tested;N
=
239580).Figure 10
The first 6 eigenfunctions (modes of vibration) of the thinly tethered membrane: (a)u
1; (b) u
2; (c) u
3; (d) u
4; (e) u
5; (f) u
6.
(a)
(b)
(c)
(d)
(e)
(f)


### 4.5. Variable Coefficients

In this section we consider a more general eigenvalue problem with nonconstant coefficients, which generalizes the eigenvalue problem for the Laplacian:(12)
∇
·
β
x
,
y
∇
u
=
λ
u
,
u
∂
Ω
1
=
0
,
∂
u
∂
n
^
∂
Ω
2
=
0
.We take β to be a scalar function of position. While the Laplacian eigenvalue problem allows us to study the fundamental modes of vibration for a homogeneous membrane, these more general operators allow for the simulation of membranes with spatially varying tension, for example, which are of great physical relevance.For our tests, we solve (12) for the unit disk with fixed boundary and β
(
x
,
y
)
=
x
2
+
y
2
+
1, using the same decomposition and set of discretizations as used previously. The absolute error is computed with respect to finest discretization. Results here, shown in Figure 11 and Table 5, are similar to those seen in the case of the Laplacian eigenvalue problem on the unit disk. We observe, most notably, that the convergence rate for the computed eigenvalues (with respect to the discretization) is roughly fourth order again. Similarly, the time spent computing a single realization of the operator and the number of calls to the operator scale sublinearly with respect to N, just as in the constant coefficient case. However, the time per L
u calculation is slightly larger here since the gradient and divergence operators are computed in succession.Table 5
Computational times and number of iterations required for evaluation of the variable coefficient elliptic operator (12) with β
(
x
,
y
)
=
x
2
+
y
2
+
1. All computations were performed on a 2.67 GHz Intel Xeon processor.
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
3950
0.009051
1242
11.2415
15800
0.023501
3072
72.1960
63200
0.109634
7610
834.3157
252800
0.370401
20630
7641.3635Figure 11
Convergence in the first 10 eigenvalues for the unit disk with the variable coefficientβ
(
x
,
y
)
=
x
2
+
y
2
+
1 is shown with respect to the finest discretization; N
=
252800.


### 4.6. Parallelization

The differential operatorL can be evaluated independently on each patch P
i, so the algorithm presented here can be parallelized by distributing this step to multiple processors. This is because the only step that requires communication between different patches is the setting of continuity conditions (via interpolation from the interiors of overlapping patches). To implement parallelization across p processors, we use a simple load balancing algorithm for n
≥
p (cf. Algorithm  2a in [21]) that assigns a given task to the processor with the least amount of work. The task in this case is computing the restriction of L
u to P
i, and the amount of work for this task is given by the total number of points in P
i.To test the efficiency of the parallelization, we consider again the thinly tethered membrane geometry shown in Figure8 (n
=
37), using a discretization with N
=
3604480. The parallel efficiency, defined by(13)
E
p
=
T
1
p
T
p
,where T
1 is the time taken to apply L
u using a single processor and T
p is the time taken to apply L
u after distributing the tasks to p processors, is shown in Figure 12 with respect to increasing p. Sharp jumps in the efficiency are most likely because of the prime number of patches, which makes distribution of tasks in as even a manner as possible somewhat difficult under this framework. Our parallel efficiency tests included runs up to twelve processors; the parallel efficiency decreased slowly as the number of processors increased to a level around 80%.Figure 12
Parallel efficiency tests performed on 12 Intel Xeon 2.66 GHz processor computer: (a) parallel efficiency (E
p) decreases with the number of processors to approximately 80%; (b) time per call to L
u is shown versus the increasing number of processors. The dashed line denotes perfect efficiency (i.e., E
p
=
1, or p
T
p
=
T
1).
(a)
(b)


## 4.1. Simple Geometries

To test the efficiency and accuracy of the eigenvalue solver as well as the domain decomposition approach we first present numerical results for two simple geometries, the unit square and the unit disk, usingL
=
∇
2. For each of these geometries, we consider both Dirichlet and Neumann homogeneous boundary conditions. The exact eigenvalues and eigenfunctions are well known and easily derived. We report them here for reference:(1)
Unit square (Ω
=
[
0,1
]
×
[
0,1
]) with all sides fixed: λ
m
n
=
π
2
(
m
2
+
n
2
) for m
,
n
∈
Z
+.
(2)
Unit square with one side (x
=
0) free: λ
m
n
=
π
2
(
(
m
+
1
/
2
)
2
+
n
2
) for m
∈
Z
+
∪
{
0
} and n
∈
Z
+.
(3)
Unit disk (Ω
=
{
(
x
,
y
)
∈
R
2
∣
x
2
+
y
2
≤
1
}) with fixed boundary: λ
m
n
=
j
n
m
2, where j
n
m is the nth root of the mth Bessel function.
(4)
Unit disk with free boundary:λ
m
n
=
j
′
2
n
m, where j
′
n
m is the nth root of the derivative of the mth Bessel function.For these cases, the domain decompositions we use are not optimal (in the sense of minimizing the number of patches), as they are meant to test how well continuity conditions are enforced. In addition, we intentionally discretize neighboring patches with slightly different numbers of grid points (to prevent grid lines from overlapping and artificially increasing the accuracy of interpolation). The decompositions used for the square and disk cases are shown in Figure2.Figure 2
Decomposition of simple geometries: (a) the unit square, decomposed into two overlapping rectangles; (b) the unit disk, decomposed into eight overlapping sections of an annulus and a square centered at the origin. These examples demonstrate the usefulness of domain decomposition to avoid difficulties in mapping entire domains.
(a)
(b)The error in the first 10 computed eigenvalues for all the cases and for increasingly finer discretizations is shown in Figure3. In each case, we observe fourth-order convergence (e.g., doubling the number of points along one direction of each patch decreases the error in the eigenvalues by no less than 2
-
4). In addition, for the finest discretization in each case, we obtain at least 8 digits of accuracy in the principal eigenvalue. Finally, the number of calls to the numerical operator (denoted by L
u) and the time taken per case are reported in Table 1.Table 1
Computational time required for evaluation of the differential operatorL
u for four test cases, with breakdown of computing time per L
u calculation, and the number of L
u calculations needed for the eigenvalue solver to converge to the required tolerance. All computations were performed on a 2.67 GHz Intel Xeon processor. (a) Square with fixed edges; (b) square with one free edge; (c) disk with fixed boundary; (d) disk with free boundary.
(a)
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
1525
0.001691
728
1.2313
6100
0.006096
2179
13.2835
24400
0.020961
6003
125.8296
97600
0.103420
22861
2364.2878
(b)
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
1525
0.001891
923
1.7452
6100
0.005459
2702
14.7507
24400
0.019813
8328
165.0070
97600
0.094591
23440
2217.2179
(c)
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
3950
0.006348
1065
6.7605
15800
0.017500
2037
35.6459
63200
0.072966
4981
363.4445
252800
0.257892
12338
3181.8739
(d)
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
3950
0.006638
1584
10.5144
15800
0.012105
3336
40.3828
63200
0.061148
7027
429.6900
252800
0.222666
14820
3299.9156Figure 3
Convergence in the first ten eigenvalues for each of the four cases in Section4.1: (a) unit square and fixed edges; (b) unit square and one free edge; (c) unit disk and fixed boundary; (d) unit disk and free boundary (with zero eigenvalues omitted). In all cases, the dashed line for fourth-order convergence is shown for comparison. N denotes the number of points used to discretize the overall two-dimensional geometry in each case.
(a)
(b)
(c)
(d)


## 4.2. Comparison to Finite Differences

Because the Arnoldi iteration is a black box algorithm that depends only on the definition of the numerical operator, one might consider a very similar algorithm that replaces the FC(Gram) algorithm with another method for computing derivatives. In this section, we use the examples of the previous section to show that the choice of FC(Gram) is justified when compared to a standard finite difference scheme and that the advantages of FC(Gram) are quite apparent. (It is worth noting the recent contribution in [20], which relies on finite differences and Richardson extrapolation, can provide highly accurate approximations of Laplace eigenvalues provided that the domain can be discretized by means of grids which exactly sample the boundary of the domain.)In these tests, we compare the number of iterations, time taken, and number of grid points needed by the two methods to compute the principal eigenvalue of the Laplacian to a fixed amount of error. We consider both the unit disk and the unit square with Dirichlet boundary conditions, and a second-order accurate centered finite difference scheme is used in the comparison. As the results in Table2 indicate, the FC(Gram) algorithm requires significantly fewer iterations and grid points than the finite difference algorithm to reach the desired accuracy. This also results in a much lower overall run time, even though the time per operation is smaller for the finite difference scheme. In addition, the FC(Gram) algorithm scales well to higher accuracies, while the cost of the finite difference algorithm grows prohibitively large for the finer discretizations. These results are even more pronounced in the case of the unit disk. Finally, simple tests using higher order Gram bases for the FC(Gram) algorithm demonstrate that the order of accuracy can be increased without much change in the computational cost. At the same time, there is no significant alteration to the FC(Gram) algorithm at high orders, while high-order finite difference schemes must be specially reformulated at boundaries to take into account the extended stencils. Thus, FC(Gram) seems to be well suited to the elliptic eigenvalue problems we consider here and competitive with other methods.Table 2
Number of iterations and computational time required by the FC solver and the finite difference solver to computeλ
1 to a fixed accuracy for the unit square example (fixed boundaries). (The first 10 eigenvalues are calculated to the tolerance of 1
0
-
7, but we only show the error in the first eigenvalue for conciseness.) For each given accuracy, the FC(Gram) algorithm outperforms the finite difference algorithm by a large margin. All computations were performed on a 2.67 GHz Intel Xeon processor. (a) Results using centered finite difference; (b) results using FC(Gram) with 5th-order basis.
(a)
Absolute error
N
Time perL
u operation (s)
Number ofL
u operations
3.00 × 10−4
74725
0.009501
9693
7.57 × 10−5
298800
0.05130
37760
(b)
Absolute error
N
Time perL
u operation (s)
Number ofL
u operations
3.53 × 10−4
1525
0.001235
728
1.41 × 10−5
6100
0.004211
2245


## 4.3. Smoothed L-Shaped Membrane

Collocation methods were used to calculate the eigenvalues of the L-shaped membrane with Dirichlet boundary conditions in order to show how polygonal boundaries, including those with geometric singularities, could be handled. Improvements to Fox et al.’s work [4] on the method of particular solutions have yielded the first several eigenvalues to at least 10 digits of accuracy [5]. In its present form the FC method does not handle geometries containing corner singularities, but it can be applied to approximations thereof containing sharp but smooth corners. (In the presence of corners in the boundary of the domain, PDE solutions and eigenfunctions are singular and therefore cannot be approximated with high accuracy by means of the present version of the FC method. An extension of this method could be envisioned which takes into account such singularities explicitly, but the development of such approaches lies beyond the scope of this paper.) We thus consider the eigenvalue problem for an L-shaped membrane with a smoothed corner, as discussed in what follows.To decompose the membrane, we smooth the corner by locally approximating the boundary curve withC
∞ curve that is sampled at a fixed set of points. Evaluation at an arbitrary location is produced by means of Neville’s algorithm for cubic interpolation on a sufficiently fine mesh to ensure that approximation errors are below the tolerances otherwise required. The smooth curve that replaces the corner is given (up to rigid transformations) by(9)
g
x
=
∫
-
∞
x
∫
-
∞
y
ϕ
z
d
z
d
y
,where the function ϕ
(
z
) is a smooth bump function on R with support on [
-
1,1
]:(10)
ϕ
z
=
0
,
if
z
∉
-
1,1
,
e
-
1
/
1
-
z
e
-
1
/
1
-
z
+
e
-
1
/
z
,
if
z
∈
-
1,1
.Thus, outside the interval [
-
1,1
], g
(
x
)
=
|
x
| and matches the rest of the boundary of the L-shaped membrane near the corner. Within it, we have a smooth curve that deviates from the corner, as shown in Figure 4. (This can be thought of as a way of approximating the Dirac delta function, which is the second “derivative” of the absolute value function, with a smooth function that has a peak of finite height and nonzero width.)Figure 4
The smoothed corner coincides exactly with the absolute value function outside of[
-
1,1
].Integral (9) is obtained numerically using Fourier techniques that take advantage of the fact that ϕ
(
z
) is C
∞ periodic function on [
-
1,1
]. We sample ϕ at 1024 equally spaced points on [
-
1,1
] and compute the FFT of the discrete data; at higher frequencies, the Fourier coefficients are effectively zero. Integration in the frequency domain amounts to multiplying each coefficient by a constant. Thus, taking the inverse FFT after multiplication yields the integral of ϕ on [
-
1,1
] up to a constant of integration. The second integral is computed the same way, with the constant term integrated analytically. This yields a highly accurate representation of g
(
x
) and its derivatives, which are necessary for the calculation of the Laplacian. Note that this method can be generalized for smoothing reentrant corners of any angle by choosing the first integration constant or adjusting the height of the peak in ϕ accordingly.Using this description of the approximate boundary curve we parametrize the region of the L-shaped membrane near the corner by(11)
x
~
s
,
t
=
2
+
2
l
s
-
l
+
1
,
y
~
s
,
t
=
t
b
g
2
+
2
l
s
-
l
+
1
,
x
y
=
d
cos
⁡
θ
-
sin
⁡
θ
sin
⁡
θ
cos
⁡
θ
x
~
y
~
+
x
0
y
0
+
h
sin
⁡
θ
cos
⁡
θ
,for s
∈
[
0,1
] and t
∈
[
0
,
b
]. Here θ
=
-
π
/
4, (
x
0
,
y
0
)
=
(
1,1
), and l and h are chosen to be of the same order as d. The decomposition of the smoothed L-shaped membrane is shown in Figure 5. This particular parametrization is chosen so that as d
→
0 the curve approximates the corner of the L-shaped membrane arbitrarily well.Figure 5
Domain decomposition for the smoothed L-shaped membrane. The corner patch (d
=
0.05) is drawn using a dashed line for emphasis.Figure7 indicates that 4th-order convergence is achieved, just as in the tests for the simpler unit square and unit disk geometries; run time information is given in Table 3. Of equal importance is the fact that the computed eigenfunctions (Figure 6) closely resemble those of the true L-shaped membrane [2] and the corresponding eigenvalues agree to anywhere from two to four digits, as one might expect. In practice, we have found that the eigenvalues converge linearly as the rounding radius of curvature tends to zero.Table 3
Computational time required for evaluation of the Laplacian for the smoothed L-shaped membrane. All computations were performed on a 2.67 GHz Intel Xeon processor.
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
9216
0.012912
6820
88.0598
18225
0.018062
14736
266.1676
36864
0.031643
21469
679.3593
72990
0.049035
44212
2167.9487
147456
0.149842
88070
13196.5794Figure 6
The first 5 eigenfunctions (modes of vibration) of the smoothed L-shaped membrane (cf. Figure  1 in [2]): (a) λ
1
=
-
9.6013; (b) λ
2
=
-
15.1969; (c) λ
3
=
-
19.7392; (d) λ
4
=
-
29.5214; (e) λ
5
=
-
31.8179.
(a)
(b)
(c)
(d)
(e)Figure 7
Convergence in the first 5 eigenvalues for the smoothed L-shaped membrane (d
=
0.05) is roughly 4th order when compared to the finest discretization tested (N
=
147456). Note that the eigenvalues do not necessarily converge to those of the true L-shaped membrane with a sharp corner due to the approximation of the corner geometry. For example, here λ
1
=
9.6013 when N
=
147456, even though λ
1
=
9.63972 (e.g., [2, 4]) for the membrane with the sharp corner.


## 4.4. Thinly Tethered Membrane

In this section we explore a more physically relevant membrane geometry inspired by quantum optics experiments. The particular geometry explored here can be described as a roughly square-shaped membrane (with side length on the order of 100μm) attached by thin tethers coming off of the four corners to a square frame (which has side length on the order of 600 μm). The edges of the membrane that are not attached to the frame are left free. In addition, the tethers smoothly continue into the boundaries of the inner square and also taper smoothly at the corners of the frame. (The purpose of this geometry is to decrease stress on the main part of the membrane, which would theoretically drive down the fundamental frequency.) For this experiment, we take the stress across the membrane to be constant in order to obtain a first approximation to the mechanics of such a complex membrane. Under this approximation, the operator in (2) is once again the Laplace operator (up to a constant factor).The decomposition developed for this membrane is shown in Figure8. Nearly all of the patches used are rectangles, but the smooth boundary curves at the corner of the membrane and the frame provide a challenge. To address this issue, we use the methodology previously developed for smoothing corners to parametrize these regions; we treat the segments of the boundary shown in Figures 8(b) and 8(c) as having smoothed reentrant corners with angles of 90° and 135°, respectively. The dimensions of the membrane are scaled so that the frame is the square of side length 2 centered at the origin.Figure 8
Domain decomposition for the thinly tethered membrane: (a) full decomposition (37 patches); (b) curving boundary near point where tether attaches to center square; (c) curving boundary near corner of frame.
(a)
(b)
(c)As in all the previous cases, we consider several successively finer discretizations for the geometry. Because this problem is significantly larger than those considered previously in this paper, we only increase the number of points in each direction by a factor of 1.2 at each level. Convergence rates for the first 5 eigenvalues are shown in Figure9 and timings are given in Table 4. In addition, plots of the eigenfunctions corresponding to these computed eigenvalues are shown in Figure 10, which show interesting oscillation modes.Table 4
Computational times and number of iterations required for evaluation of the Laplacian on the thinly tethered membrane. All computations were performed on a 2.67 GHz Intel Xeon processor to a tolerance of1
0
-
4.
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
56320
0.061338
464727
28505.4983
79420
0.072279
1758160
127079.4402
116380
0.115435
2679346
309291.7615
166045
0.167378
2793136
467510.6446Figure 9
Convergence in the first 5 eigenvalues for the thinly tethered membrane is roughly 4th order (when compared to the finest discretization tested;N
=
239580).Figure 10
The first 6 eigenfunctions (modes of vibration) of the thinly tethered membrane: (a)u
1; (b) u
2; (c) u
3; (d) u
4; (e) u
5; (f) u
6.
(a)
(b)
(c)
(d)
(e)
(f)


## 4.5. Variable Coefficients

In this section we consider a more general eigenvalue problem with nonconstant coefficients, which generalizes the eigenvalue problem for the Laplacian:(12)
∇
·
β
x
,
y
∇
u
=
λ
u
,
u
∂
Ω
1
=
0
,
∂
u
∂
n
^
∂
Ω
2
=
0
.We take β to be a scalar function of position. While the Laplacian eigenvalue problem allows us to study the fundamental modes of vibration for a homogeneous membrane, these more general operators allow for the simulation of membranes with spatially varying tension, for example, which are of great physical relevance.For our tests, we solve (12) for the unit disk with fixed boundary and β
(
x
,
y
)
=
x
2
+
y
2
+
1, using the same decomposition and set of discretizations as used previously. The absolute error is computed with respect to finest discretization. Results here, shown in Figure 11 and Table 5, are similar to those seen in the case of the Laplacian eigenvalue problem on the unit disk. We observe, most notably, that the convergence rate for the computed eigenvalues (with respect to the discretization) is roughly fourth order again. Similarly, the time spent computing a single realization of the operator and the number of calls to the operator scale sublinearly with respect to N, just as in the constant coefficient case. However, the time per L
u calculation is slightly larger here since the gradient and divergence operators are computed in succession.Table 5
Computational times and number of iterations required for evaluation of the variable coefficient elliptic operator (12) with β
(
x
,
y
)
=
x
2
+
y
2
+
1. All computations were performed on a 2.67 GHz Intel Xeon processor.
N
Time perL
u operation (s)
Number ofL
u operations
Total time inL
u (s)
3950
0.009051
1242
11.2415
15800
0.023501
3072
72.1960
63200
0.109634
7610
834.3157
252800
0.370401
20630
7641.3635Figure 11
Convergence in the first 10 eigenvalues for the unit disk with the variable coefficientβ
(
x
,
y
)
=
x
2
+
y
2
+
1 is shown with respect to the finest discretization; N
=
252800.


## 4.6. Parallelization

The differential operatorL can be evaluated independently on each patch P
i, so the algorithm presented here can be parallelized by distributing this step to multiple processors. This is because the only step that requires communication between different patches is the setting of continuity conditions (via interpolation from the interiors of overlapping patches). To implement parallelization across p processors, we use a simple load balancing algorithm for n
≥
p (cf. Algorithm  2a in [21]) that assigns a given task to the processor with the least amount of work. The task in this case is computing the restriction of L
u to P
i, and the amount of work for this task is given by the total number of points in P
i.To test the efficiency of the parallelization, we consider again the thinly tethered membrane geometry shown in Figure8 (n
=
37), using a discretization with N
=
3604480. The parallel efficiency, defined by(13)
E
p
=
T
1
p
T
p
,where T
1 is the time taken to apply L
u using a single processor and T
p is the time taken to apply L
u after distributing the tasks to p processors, is shown in Figure 12 with respect to increasing p. Sharp jumps in the efficiency are most likely because of the prime number of patches, which makes distribution of tasks in as even a manner as possible somewhat difficult under this framework. Our parallel efficiency tests included runs up to twelve processors; the parallel efficiency decreased slowly as the number of processors increased to a level around 80%.Figure 12
Parallel efficiency tests performed on 12 Intel Xeon 2.66 GHz processor computer: (a) parallel efficiency (E
p) decreases with the number of processors to approximately 80%; (b) time per call to L
u is shown versus the increasing number of processors. The dashed line denotes perfect efficiency (i.e., E
p
=
1, or p
T
p
=
T
1).
(a)
(b)


## 5. Conclusions

We have presented a high-order method for evaluation of the eigenvalues of elliptic operators. Our approach combines the FC method for evaluation of accurate Fourier approximations of nonperiodic functions, domain decomposition, and Arnoldi iteration. Our domain decomposition strategy enables consideration of general geometries. Additionally, the decomposition strategy leads easily to efficient parallel computations. The resulting eigenvalue algorithm converges at a roughly fourth-order rate in all cases considered. Furthermore, a methodology to deal with corner singularities is suggested by the results from the smoothed L-shaped membrane. Comparisons to solutions provided by a finite difference scheme demonstrate the advantages provided by the proposed algorithm as well as the viability of our Fourier continuation algorithm for evaluation of eigenvalues for complex geometries and general variable-coefficient operators.



---

*Source: 184786-2015-11-24.xml*