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.
2.5. Numerical Differentiation
If is a 1-dimensional function whose continuation (computed via FC(Gram)) is , then the derivative of the Fourier series of is computed as follows:The derivative of is recovered by restricting the domain of . 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.
| ||||||||||||||||||||
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 mapping , the boundaries of correspond exactly to the boundaries of , 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 . 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 ):The values and depend only on the unit outward normal and the partial derivatives of components of (for a given patch ).
By treating the points on the boundary, where a Neumann condition is prescribed as unknowns, we will already have data from the candidate function along the boundary. As a result, either or 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. Suppose overlaps with near one of its boundaries. A layer of fringe points is added to the patch on the overlap region. Then the function values for on the fringe points are prescribed by interpolating from the function values in the interior of (avoiding values in 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 (when is not explicitly known) for a given point at the boundary of .
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 as an upper Hessenberg matrix . The eigenvalues of (“Ritz values”) are, in practice, good estimates of the eigenvalues of , as are the eigenvectors (up to an orthogonal transformation). The special structure of 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 rather than the cost of computing the eigenvalues of a large 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 using our FC algorithm. All eigenvalues in the following numerical experiments are computed to a tolerance of 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, using . 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 () with all sides fixed: for .(2)Unit square with one side () free: for and .(3)Unit disk () with fixed boundary: , where is the th root of the th Bessel function.(4)Unit disk with free boundary: , where is the th root of the derivative of the th 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 Figure 2.
(a)
(b)
The error in the first 10 computed eigenvalues for all the cases and for increasingly finer discretizations is shown in Figure 3. 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 ). 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 ) and the time taken per case are reported in Table 1.
(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 Table 2 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.
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 with 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) bywhere the function is a smooth bump function on with support on :Thus, outside the interval , 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.)
Integral (9) is obtained numerically using Fourier techniques that take advantage of the fact that is periodic function on . We sample at 1024 equally spaced points on 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 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 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 byfor and . Here , , and and are chosen to be of the same order as . The decomposition of the smoothed L-shaped membrane is shown in Figure 5. This particular parametrization is chosen so that as the curve approximates the corner of the L-shaped membrane arbitrarily well.
Figure 7 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.
(a)
(b)
(c)
(d)
(e)
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 Figure 8. 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.
(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 Figure 9 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.
(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: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 , 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 , just as in the constant coefficient case. However, the time per calculation is slightly larger here since the gradient and divergence operators are computed in succession.
4.6. Parallelization
The differential operator can be evaluated independently on each patch , 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 processors, we use a simple load balancing algorithm for (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 to , and the amount of work for this task is given by the total number of points in .
To test the efficiency of the parallelization, we consider again the thinly tethered membrane geometry shown in Figure 8 (), using a discretization with . The parallel efficiency, defined bywhere is the time taken to apply using a single processor and is the time taken to apply after distributing the tasks to processors, is shown in Figure 12 with respect to increasing . 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%.
(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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge support from NSF and AFOSR.