Michael O'Neil

Andreas Klöckner, Alexander Barnett, Leslie Greengard et al.

Integral equation methods for the solution of partial differential equations, when coupled with suitable fast algorithms, yield geometrically flexible, asymptotically optimal and well-conditioned schemes in either interior or exterior domains. The practical application of these methods, however, requires the accurate evaluation of boundary integrals with singular, weakly singular or nearly singular kernels. Historically, these issues have been handled either by low-order product integration rules (computed semi-analytically), by singularity subtraction/cancellation, by kernel regularization and asymptotic analysis, or by the construction of special purpose "generalized Gaussian quadrature" rules. In this paper, we present a systematic, high-order approach that works for any singularity (including hypersingular kernels), based only on the assumption that the field induced by the integral operator is locally smooth when restricted to either the interior or the exterior. Discontinuities in the field across the boundary are permitted. The scheme, denoted QBX (quadrature by expansion), is easy to implement and compatible with fast hierarchical algorithms such as the fast multipole method. We include accuracy tests for a variety of integral operators in two dimensions on smooth and corner domains.

Charles L. Epstein, Leslie Greengard, Michael O'Neil

In this paper, we develop a new integral representation for the solution of the time harmonic Maxwell equations in media with piecewise constant dielectric permittivity and magnetic permeability in R^3. This representation leads to a coupled system of Fredholm integral equations of the second kind for four scalar densities supported on the material interface. Like the classical Muller equation, it has no spurious resonances. Unlike the classical approach, however, the representation does not suffer from low frequency breakdown. We illustrate the performance of the method with numerical examples.

Michael O'Neil, Leslie Greengard, Andras Pataki

A classical problem in acoustic (and electromagnetic) scattering concerns the evaluation of the Green's function for the Helmholtz equation subject to impedance boundary conditions on a half-space. The two principal approaches used for representing this Green's function are the Sommerfeld integral and the (closely related) method of complex images. The former is extremely efficient when the source is at some distance from the half-space boundary, but involves an unwieldy range of integration as the source gets closer and closer. Complex image-based methods, on the other hand, can be quite efficient when the source is close to the boundary, but they do not easily permit the use of the superposition principle since the selection of complex image locations depends on both the source and the target. We have developed a new, hybrid representation which uses a finite number of real images (dependent only on the source location) coupled with a rapidly converging Sommerfeld-like integral. While our method applies in both two and three dimensions, we restrict the detailed analysis and numerical experiments here to the two-dimensional case.

Andras Pataki, Antoine J. Cerfon, Jeffrey P. Freidberg et al.

We present a new fast solver to calculate fixed-boundary plasma equilibria in toroidally axisymmetric geometries. By combining conformal mapping with Fourier and integral equation methods on the unit disk, we show that high-order accuracy can be achieved for the solution of the equilibrium equation and its first and second derivatives. Smooth arbitrary plasma cross-sections as well as arbitrary pressure and poloidal current profiles are used as initial data for the solver. Equilibria with large Shafranov shifts can be computed without difficulty. Spectral convergence is demonstrated by comparing the numerical solution with a known exact analytic solution. A fusion-relevant example of an equilibrium with a pressure pedestal is also presented.

Charles L. Epstein, Leslie Greengard, Michael O'Neil

The generalized Debye source representation of time-harmonic electromagnetic fields yields well-conditioned second-kind integral equations for a variety of boundary value problems, including the problems of scattering from perfect electric conductors and dielectric bodies. Furthermore, these representations, and resulting integral equations, are fully stable in the static limit as $ω\to 0$ in multiply connected geometries. In this paper, we present the first high-order accurate solver based on this representation for bodies of revolution. The resulting solver uses a Nyström discretization of a one-dimensional generating curve and high-order integral equation methods for applying and inverting surface differentials. The accuracy and speed of the solvers are demonstrated in several numerical examples.

Charles L. Epstein, Zydrunas Gimbutas, Leslie Greengard et al.

A classical problem in electromagnetics concerns the representation of the electric and magnetic fields in the low-frequency or static regime, where topology plays a fundamental role. For multiply connected conductors, at zero frequency the standard boundary conditions on the tangential components of the magnetic field do not uniquely determine the vector potential. We describe a (gauge-invariant) consistency condition that overcomes this non-uniqueness and resolves a longstanding difficulty in inverting the magnetic field integral equation.

Jun Lai, Leslie Greengard, Michael O'Neil

Scattering from large, open cavity structures is of importance in a variety of electromagnetic applications. In this paper, we propose a new well conditioned integral equation for scattering from general open cavities embedded in an infinite, perfectly conducting half-space. The integral representation permits the stable evaluation of both the electric and magnetic field, even in the low-frequency regime, using the continuity equation in a post-processing step. We establish existence and uniqueness results, and demonstrate the performance of the scheme in the cavity-of-revolution case. High-order accuracy is obtained using a Nystrom discretization with generalized Gaussian quadratures.

Manas Rachh, Andreas Klöckner, Michael O'Neil

The use of integral equation methods for the efficient numerical solution of PDE boundary value problems requires two main tools: quadrature rules for the evaluation of layer potential integral operators with singular kernels, and fast algorithms for solving the resulting dense linear systems. Classically, these tools were developed separately. In this work, we present a unified numerical scheme based on coupling Quadrature by Expansion, a recent quadrature method, to a customized Fast Multipole Method (FMM) for the Helmholtz equation in two dimensions. The method allows the evaluation of layer potentials in linear-time complexity, anywhere in space, with a uniform, user-chosen level of accuracy as a black-box computational method. Providing this capability requires geometric and algorithmic considerations beyond the needs of standard FMMs as well as careful consideration of the accuracy of multipole translations. We illustrate the speed and accuracy of our method with various numerical examples. Keywords: Layer Potentials; Singular Integrals; Quadrature; High-order accuracy; Integral equations; Helmholtz equation; Fast multipole method.

Charles L. Epstein, Michael O'Neil

We introduce an arbitrary order, computationally efficient method to smooth corners on curves in the plane, as well as edges and vertices on surfaces in $\mathbb R^3$. The method is local, only modifying the original surface in a neighborhood of the geometric singularity, and preserves desirable features like convexity and symmetry. The smoothness of the final surface is an explicit parameter in the method, and the bandlimit of the smoothed surface is proportional to its smoothness. Several numerical examples are provided in the context of acoustic scattering. In particular, we compare scattered fields from smoothed geometries in two dimensions with those from polygonal domains. We observe that significant reductions in computational cost can be obtained if merely approximate solutions are desired in the near- or far-field. Provided that it is sub-wavelength, the error of the scattered field is proportional to the size of the geometry that is modified.

Dhairya Malhotra, Antoine Cerfon, Lise-Marie Imbert-Gérard et al.

We present a boundary integral equation solver for computing Taylor relaxed states in non-axisymmetric solid and shell-like toroidal geometries. The computation of Taylor states in these geometries is a key element for the calculation of stepped pressure stellarator equilibria. The integral representation of the magnetic field in this work is based on the generalized Debye source formulation, and results in a well-conditioned second-kind boundary integral equation. The integral equation solver is based on a spectral discretization of the geometry and unknowns, and the computation of the associated weakly-singular integrals is performed with high-order quadrature based on a partition of unity. The resulting scheme for applying the integral operator is then coupled with an iterative solver and suitable preconditioners. Several numerical examples are provided to demonstrate the accuracy and efficiency of our method, and a direct comparison with the leading code in the field is reported.

Michael O'Neil

The Laplace-Beltrami problem $Δ_Γψ= f$ has several applications in mathematical physics, differential geometry, machine learning, and topology. In this work, we present novel second-kind integral equations for its solution which obviate the need for constructing a suitable parametrix to approximate the in-surface Green's function. The resulting integral equations are well-conditioned and compatible with standard fast multipole methods and iterative linear algebraic solvers, as well as more modern fast direct solvers. Using layer-potential identities known as Calderón projectors, the Laplace-Beltrami operator can be pre-conditioned from the left and/or right to obtain second-kind integral equations. We demonstrate the accuracy and stability of the scheme in several numerical examples along surfaces described by curvilinear triangles.

Michael O'Neil, Antoine J. Cerfon

We develop an algorithm for the numerical calculation of Taylor states in toroidal and toroidal shell geometries using an analytical framework developed for the solution to the time-harmonic Maxwell equations. Taylor states are a special case of what are known as Beltrami fields, or linear force-free fields. The scheme of this work relies on the generalized Debye source representation of Maxwell fields and an integral representation of Beltrami fields which immediately yields a well-conditioned second-kind integral equation. This integral equation has a unique solution whenever the Beltrami parameter $λ$ is not a member of a discrete, countable set of resonances which physically correspond to spontaneous symmetry breaking. Several numerical examples relevant to magnetohydrodynamic equilibria calculations are provided. Lastly, our approach easily generalizes to arbitrary geometries, both bounded and unbounded, and of varying genus.

Paul G. Beckman, Samuel F. Potter, Michael O'Neil

The eigenfunctions of the Laplacian are a natural basis of functions for many tasks in computational mathematics. On the circle and sphere, the eigenfunctions are given by complex periodic exponentials and spherical harmonics, respectively, and much work has been done to develop fast algorithms for analyzing and synthesizing data in these bases. In this work, we generalize these special-case transforms to Laplace-Beltrami eigenfunctions of arbitrary surfaces, referred to as manifold harmonics. The resulting fast algorithm for computing linear combinations of the manifold harmonics is based on a butterfly factorization, which hierarchically compresses the transform matrix by constructing nested low-rank approximations of carefully selected submatrices. Several numerical examples are provided which demonstrate the speedups and reduction in memory requirements achieved by our algorithm for a variety of geometries, discretizations, and applications. In addition, a detailed analysis of the algorithm is given in the case that the underlying manifold is the flat periodic square.

Jun Lai, Michael O'Neil

Fast, high-order accurate algorithms for electromagnetic scattering from axisymmetric objects are of great importance when modeling physical phenomena in optics, materials science (e.g. meta-materials), and many other fields of applied science. In this paper, we develop an FFT-accelerated separation of variables solver that can be used to efficiently invert integral equation formulations of Maxwell's equations for scattering from axisymmetric penetrable (dielectric) bodies. Using a standard variant of Müller's integral representation of the fields, our numerical solver rapidly and directly inverts the resulting second-kind integral equation. In particular, the algorithm of this work (1) rapidly evaluates the modal Green's functions, and their derivatives, via kernel splitting and the use of novel recursion formulas, (2) discretizes the underlying integral equation using generalized Gaussian quadratures on adaptive meshes, and (3) is applicable to geometries containing edges. Several numerical examples are provided to demonstrate the efficiency and accuracy of the aforementioned algorithm in various geometries.

Jun Lai, Leslie Greengard, Michael O'Neil

A variety of problems in acoustic and electromagnetic scattering require the evaluation of impedance or layered media Green's functions. Given a point source located in an unbounded half-space or an infinitely extended layer, Sommerfeld and others showed that Fourier analysis combined with contour integration provides a systematic and broadly effective approach, leading to what is generally referred to as the Sommerfeld integral representation. When either the source or target is at some distance from an infinite boundary, the number of degrees of freedom needed to resolve the scattering response is very modest. When both are near an interface, however, the Sommerfeld integral involves a very large range of integration and its direct application becomes unwieldy. Historically, three schemes have been employed to overcome this difficulty: the method of images, contour deformation, and asymptotic methods of various kinds. None of these methods make use of classical layer potentials in physical space, despite their advantages in terms of adaptive resolution and high-order accuracy. The reason for this is simple: layer potentials are impractical in layered media or half-space geometries since they require the discretization of an infinite boundary. In this paper, we propose a hybrid method which combines layer potentials (physical-space) on a finite portion of the interface together with a Sommerfeld-type (Fourier) correction. We prove that our method is efficient and rapidly convergent for arbitrarily located sources and targets, and show that the scheme is particularly effective when solving scattering problems for objects which are close to the half-space boundary or even embedded across a layered media interface.

Sivaram Ambikasaran, Daniel Foreman-Mackey, Leslie Greengard et al.

A number of problems in probability and statistics can be addressed using the multivariate normal (Gaussian) distribution. In the one-dimensional case, computing the probability for a given mean and variance simply requires the evaluation of the corresponding Gaussian density. In the $n$-dimensional setting, however, it requires the inversion of an $n \times n$ covariance matrix, $C$, as well as the evaluation of its determinant, $\det(C)$. In many cases, such as regression using Gaussian processes, the covariance matrix is of the form $C = σ^2 I + K$, where $K$ is computed using a specified covariance kernel which depends on the data and additional parameters (hyperparameters). The matrix $C$ is typically dense, causing standard direct methods for inversion and determinant evaluation to require $\mathcal O(n^3)$ work. This cost is prohibitive for large-scale modeling. Here, we show that for the most commonly used covariance functions, the matrix $C$ can be hierarchically factored into a product of block low-rank updates of the identity matrix, yielding an $\mathcal O (n\log^2 n) $ algorithm for inversion. More importantly, we show that this factorization enables the evaluation of the determinant $\det(C)$, permitting the direct calculation of probabilities in high dimensions under fairly broad assumptions on the kernel defining $K$. Our fast algorithm brings many problems in marginalization and the adaptation of hyperparameters within practical reach using a single CPU core. The combination of nearly optimal scaling in terms of problem size with high-performance computing resources will permit the modeling of previously intractable problems. We illustrate the performance of the scheme on standard covariance kernels.