Mathematical Problems in Engineering

Volume 2015, Article ID 184786, 15 pages

http://dx.doi.org/10.1155/2015/184786

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

Applied and Computational Mathematics, California Institute of Technology, Pasadena, CA 91125, USA

Received 2 July 2015; Accepted 25 October 2015

Academic Editor: Salvatore Alfonzetti

Copyright © 2015 Oscar P. Bruno et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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: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 operations, where 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 by and . 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 be a bounded domain and let be a linear elliptic operator of the form , where is a two-dimensional multi-index. Throughout this paper, we search for a subset of the scalar values and corresponding functions such thatwhere are the eigenvalues and are the eigenfunctions. The boundary conditions are defined so that is a disjoint union that partitions and 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

Let be open sets relative to such that there exist invertible smooth functions , where . We let denote coordinates on , which correspond to on under the mapping . Using the chain rule to translate derivatives in domain coordinates to derivatives in reference coordinates (e.g., ) and inverting the resulting system yields the desired expressions for derivatives in and in terms of and . The derivatives in and 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 on modes is . 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 function , produces Fourier series for smooth periodic extensions of to a suitably larger interval. Thus, given a function defined on (say) which is sampled on the discrete grid given by , where and , the algorithm constructs a new function on for some such that , , and , . In view of these conditions the function is smooth and periodic, and therefore its -periodic trigonometric polynomial,(where is the set of modes distributed symmetrically about 0) is highly accurate. To enforce the approximation of the function , one must compute the coefficients in (3) such that and agree on the original sample points. This could, for example, be accomplished by solving a (possibly oversampled) minimization problem:Written as a matrix equation, this amounts to solvingin the least-squares sense, where is a column vector whose entries are the coefficients , is a column vector whose entries are , and is a matrix with .

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 of . 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 , where and are fixed integers chosen independently of . By computing the trigonometric polynomial fit via (4) on this set of points, we can construct a function such thatis -periodic. A Fourier approximation of the smooth and periodic function in the interval (which also provides an approximation for in the interval ) is highly accurate, and it provides, in particular, the needed Fourier expansion of the function in the interval .

It is possible to select and 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 and . 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 and a Gram polynomial basis up to degree 5.