Generalized impedance boundary conditions for the Helmholtz equation arise in a number of acoustic and electromagnetic models that approximate boundary effects, including thin coatings and the viscous and thermal losses that can occur in boundary layers. Here we present a well-conditioned numerical method for this class of boundary conditions based on a boundary integral re-formulation of the equations. The formulation applies a combined layer representation together with a surface preconditioner to obtain an integral equation of the second kind. Additionally, by using appropriate image sources, we construct well-conditioned representations of impedance-to-impedance maps for subdomains carrying these boundary conditions on part of their boundary, so that the method can be used within a domain-decomposition framework. Several numerical examples in geometries inspired by acoustic applications demonstrate the efficacy of the method.

The dynamics of surface waves traveling along the boundary of a liquid medium are changed by the presence of floating plates and membranes, contributing to a number of important phenomena in a wide range of applications. Mathematically, if the fluid is only partly covered by a plate or membrane, the order of derivatives of the surface-boundary conditions jump between regions of the surface. In this work, we consider a general class of problems for infinite depth linearized surface waves in which the plate or membrane has a compact hole or multiple holes. For this class of problems, we describe a general integral equation approach, and for two important examples, the partial membrane and the polynya, we analyze the resulting boundary integral equations. In particular, we show that they are Fredholm second kind and discuss key properties of their solutions. We develop flexible and fast algorithms for discretizing and solving these equations, and demonstrate their robustness and scalability in resolving surface wave phenomena through several numerical examples.

In this paper, we develop second kind integral formulations for flexural wave scattering problems involving the clamped, supported, and free plate boundary conditions. While the clamped plate problem can be solved with layer potentials developed for the biharmonic equation, the free plate problem is more difficult due to the order and complexity of the boundary conditions. In this work, we describe a representation for the free plate problem that uses the Hilbert transform to cancel singularities of certain layer potentials, ultimately leading to a Fredholm integral equation of the second kind. Additionally, for the supported plate problem, we improve on an existing representation to obtain a second kind integral equation formulation. With these representations it is possible to solve flexural wave scattering problems with high-order-accurate methods, examine the far field patterns of scattering objects, and solve large problems involving multiple scatterers.

In this work, we develop a fast and accurate method for the scattering of flexural-gravity waves by a thin plate of varying thickness overlying a fluid of infinite depth. This problem commonly arises in the study of sea ice and ice shelves, which can have complicated heterogeneities that include ridges and rolls. With certain natural assumptions on the thickness, we present an integral equation formulation for solving this class of problems and analyze its mathematical properties. The integral equation is then discretized and solved using a high-order-accurate, FFT-accelerated algorithm. The speed, accuracy, and scalability of this approach are demonstrated through a variety of illustrative examples.

In inverse scattering problems, a model that allows for the simultaneous recovery of both the domain shape and an impedance boundary condition covers a wide range of problems with impenetrable domains, including recovering the shape of sound-hard and sound-soft obstacles and obstacles with thin coatings. This work develops an optimization framework for recovering the shape and material parameters of a penetrable, dissipative obstacle in the multifrequency setting, using a constrained class of curvature-dependent impedance function models proposed by Antoine, Barucq, and Vernhet. We find that this constrained model improves the robustness of the recovery problem, compared to more general models, and provides meaningfully better obstacle recovery than simpler models. We explore the effectiveness of the model for varying levels of dissipation, for noise-corrupted data, and for limited aperture data in the numerical examples.

A new scheme is presented for imposing periodic boundary conditions on unit cells with arbitrary source distributions. We restrict our attention here to the Poisson, modified Helmholtz, Stokes and modified Stokes equations. The approach extends to the oscillatory equations of mathematical physics, including the Helmholtz and Maxwell equations, but we will address these in a companion paper, since the nature of the problem is somewhat different and includes the consideration of quasiperiodic boundary conditions and resonances. Unlike lattice sum-based methods, the scheme is insensitive to the unit cell's aspect ratio and is easily coupled to adaptive fast multipole methods (FMMs). Our analysis relies on classical “plane-wave” representations of the fundamental solution, and yields an explicit low-rank representation of the field due to all image sources beyond the first layer of neighboring unit cells. When the aspect ratio of the unit cell is large, our scheme can be coupled with the nonuniform fast Fourier transform (NUFFT) to accelerate the evaluation of the induced field. Its performance is illustrated with several numerical examples.

Inverse obstacle scattering is the recovery of an obstacle boundary from the scattering data produced by incident waves. This shape recovery can be done by iteratively solving a PDE-constrained optimization problem for the obstacle boundary. While it is well known that this problem is typically non-convex and ill-posed, previous investigations have shown that in many settings these issues can be alleviated by using a continuation-in-frequency method and introducing a regularization that limits the frequency content of the obstacle boundary. It has been recently observed that these techniques can fail for obstacles with pronounced cavities, even in the case of penetrable obstacles where similar optimization and regularization methods work for the equivalent problem of recovering a piecewise constant wave speed. The present work investigates the recovery of obstacle boundaries for impenetrable, sound-soft media with pronounced cavities, given multi-frequency scattering data. Numerical examples demonstrate that the problem is sensitive to the choice of iterative solver used at each frequency and the initial guess at the lowest frequency. We propose a modified continuation-in-frequency method which follows a random walk in frequency, as opposed to the standard monotonically increasing path. This method shows some increased robustness in recovering cavities, but can also fail for more extreme examples. An interesting phenomenon is observed that while the obstacle reconstructions obtained over several random trials can vary significantly near the cavity, the results are consistent for non-cavity parts of the boundary.

With the funding provided by this award, we developed numerical codes for the study of magnetically confined plasmas for fusion applications. Accordingly, our work can be divided into two separate categories: 1) the design and analysis of novel numerical methods providing high accuracy and high efficiency; 2) the study of the equilibrium and stability of magnetically confined plasmas with some of these numerical codes, as well as the study of the nature of the turbulent behavior which may arise in the presence of instabilities. We first developed new numerical schemes based on integral equation methods for the computation of steady-state magnetic configurations in fusion experiments, providing high accuracy for the magnetic field and its derivatives, which are required for stability and turbulence calculations. We employed different integral formulations depending on the application of interest: axisymmetric or non-axisymmetric equilibria, force-free or magnetohydrodynamic equilibria, fixed-boundary equilibria or free-boundary equilibria. While efficient, these methods do not yet apply to plasma boundaries which are not smooth, a situation which is fairly common in magnetic confinement experiments. To address this temporary weakness, we also constructed a new steady-state solver based on the Hybridizable Discontinuous Galerkin (HDG) method, which provides full geometric flexibility. In addition to these numerical tools focused on steady-states, we also contributed to the improvement of the speed and accuracy of codes simulating the plasma dynamics of fusion plasmas, by developing a novel velocity space representation for the efficient solution of kinetic equations, which most accurately describe the time evolution of hot plasmas in fusion experiments. Using the tools discussed above, we studied several questions pertaining to the equilibrium and stability of magnetically confined plasmas. In particular, we derived a new simple model for axisymmetric devices called tokamaks, to predict how elongated a fusion plasma can be before it becomes unstable and collapses. We also looked at the effect of the shape of the outer plasma surface on key properties of the steady-state magnetic configurations, and how these properties impact turbulence in fusion plasmas, and the corresponding transport of momentum. Likewise, we studied the role of large localized flows on the steady-state magnetic configurations, and how they may influence plasma stability and turbulence. Non-axisymmetric steady-state magnetic configurations are inherently more complex than axisymmetric steady-state configurations, and the subject of ongoing controversies regarding the regularity of the equations determining such steady-states, and their solutions. Implementing an existing NYU code in a new geometry, we studied the nature of the singularity of the solutions observed in the code, and methods to eliminate them. Our main conclusion is that by appropriately tailoring the plasma boundary, it is possible to eliminate the singularities otherwise appearing in our simulations, and to obtain steady-states which appear to be smooth. To gain further insights on incompletely understood turbulence phenomena, we proposed a new reduced model capturing most of these phenomena, which is simple enough to not require expensive numerical simulations on massive supercomputers to investigate them. We demonstrated the strong similarity between our simulations and published results obtained from computationally expensive simulations, and plan to rely on our reduced model to identify the key mechanisms determining the evolution and strength turbulent driven transport in fusion plasmas. Finally, we proposed a new framework for tokamak reactor design studies, enabling us to consider the relative merits of steady-state versus pulsed fusion reactors. We found that pulsed fusion reactors may benefit most from recent advances in magnet technology, and the availability of very high field magnets. As such, they may become more desirable than steady-state tokamak reactors for cost efficient electricity generation.

The dynamic mode decomposition (DMD) is a broadly applicable dimensionality reduction algorithm that approximates a matrix containing time-series data by the outer product of a matrix of exponentials, representing Fourier-like time dynamics, and a matrix of coefficients, representing spatial structures. This interpretable spatio-temporal decomposition is commonly computed using linear algebraic techniques in its simplest formulation or a nonlinear optimization procedure within the variable projection framework. For data with sparse outliers or data which are not well-represented by exponentials in time, the standard Frobenius norm fit of the data creates significant biases in the recovered time dynamics. As a result, practitioners are left to clean such defects from the data manually or to use a black-box cleaning approach like robust PCA. As an alternative, we propose a framework and a set of algorithms for incorporating robust features into the nonlinear optimization used to compute the DMD itself. The algorithms presented are flexible, allowing for regu- larizers and constraints on the optimization, and scalable, using a stochastic approach to decrease the computational cost for data in high dimensional space. Both synthetic and real data examples are provided.

A new scheme is presented for imposing periodic boundary conditions on unit cells with arbitrary source distributions. We restrict our attention here to the Poisson, modified Helmholtz, Stokes and modified Stokes equations. The approach extends to the oscillatory equations of mathematical physics, including the Helmholtz and Maxwell equations, but we will address these in a companion paper, since the nature of the problem is somewhat different and includes the consideration of quasiperiodic boundary conditions and resonances. Unlike lattice sum-based methods, the scheme is insensitive to the unit cell's aspect ratio and is easily coupled to adaptive fast multipole methods (FMMs). Our analysis relies on classical "plane-wave" representations of the fundamental solution, and yields an explicit low-rank representation of the field due to all image sources beyond the first layer of neighboring unit cells. When the aspect ratio of the unit cell is large, our scheme can be coupled with the nonuniform fast Fourier transform (NUFFT) to accelerate the evaluation of the induced field. Its performance is illustrated with several numerial examples.

The eigenvalues and eigenfunctions of the Stokes operator have been the subject of intense analytical investigation and have applications in the study and simulation of the Navier–Stokes equations. As the Stokes operator is second order and has the divergence-free constraint, computing these eigenvalues and the corresponding eigenfunctions is a challenging task, particularly in complex geometries and at high frequencies. The boundary integral equation (BIE) framework provides robust and scalable eigenvalue computations due to (a) the reduction in the dimension of the problem to be discretized and (b) the absence of high-frequency “pollution” when using Green’s function to represent propagating waves. In this paper, we detail the theoretical justification for a BIE approach to the Stokes eigenvalue problem on simply- and multiply-connected planar domains, which entails a treatment of the uniqueness theory for oscillatory Stokes equations on exterior domains. Then, using well-established techniques for discretizing BIEs, we present numerical results which confirm the analytical claims of the paper and demonstrate the efficiency of the overall approach.

The integral equation approach to partial differential equations (PDEs) provides significant advantages in the numerical solution of the incompressible Navier-Stokes equations. In particular, the divergence-free condition and boundary conditions are handled naturally, and the ill-conditioning caused by high order terms in the PDE is preconditioned analytically. Despite these advantages, the adoption of integral equation methods has been slow due to a number of difficulties in their implementation. This work describes a complete integral equation-based flow solver that builds on recently developed methods for singular quadrature and the solution of PDEs on complex domains, in combination with several more well-established numerical methods. We apply this solver to flow problems on a number of geometries, both simple and challenging, studying its convergence properties and computational performance. This serves as a demonstration that it is now relatively straightforward to develop a robust, efficient, and flexible Navier-Stokes solver, using integral equation methods.

The problem of optimally placing sensors under a cost constraint arises naturally in the design of industrial and commercial products, as well as in scientific experiments. We consider a relaxation of the full optimization formulation of this problem and then extend a well-established greedy algorithm for the optimal sensor placement problem without cost constraints. We demonstrate the effectiveness of this algorithm on the datasets related to facial recognition, climate science, and fluid mechanics. This algorithm is scalable and often identifies sparse sensors with near-optimal reconstruction performance, while dramatically reducing the overall cost of the sensors. We find that the cost-error landscape varies by application, with intuitive connections to the underlying physics. In addition, we include experiments for various pre-processing techniques and find that a popular technique based on the singular value decomposition is often suboptimal.

Regularized regression problems are ubiquitous in statistical modeling, signal processing, and machine learning. Sparse regression, in particular, has been instrumental in scientific model discovery, including compressed sensing applications, variable selection, and high-dimensional analysis. We propose a broad framework for sparse relaxed regularized regression, called SR3. The key idea is to solve a relaxation of the regularized problem, which has three advantages over the state-of-the-art: 1) solutions of the relaxed problem are superior with respect to errors, false positives, and conditioning; 2) relaxation allows extremely fast algorithms for both convex and nonconvex formulations; and 3) the methods apply to composite regularizers, essential for total variation (TV) as well as sparsity-promoting formulations using tight frames. We demonstrate the advantages of SR3 (computational efficiency, higher accuracy, faster convergence rates, and greater flexibility) across a range of regularized regression problems with synthetic and real data, including applications in compressed sensing, LASSO, matrix completion, TV regularization, and group sparsity. Following standards of reproducible research, we also provide a companion MATLAB package that implements these examples.

Hybrid systems are traditionally difficult to identify and analyse using classical dynamical systems theory. Moreover, recently developed model identification methodologies largely focus on identifying a single set of governing equations solely from measurement data. In this article, we develop a new methodology, Hybrid-Sparse Identification of Nonlinear Dynamics, which identifies separate nonlinear dynamical regimes, employs information theory to manage uncertainty and characterizes switching behaviour. Specifically, we use the nonlinear geometry of data collected from a complex system to construct a set of coordinates based on measurement data and augmented variables. Clustering the data in these measurement-based coordinates enables the identification of nonlinear hybrid systems. This methodology broadly empowers nonlinear system identification without constraining the data locally in time and has direct connections to hybrid systems theory. We demonstrate the success of this method on numerical examples including a mass-spring hopping model and an infectious disease model. Characterizing complex systems that switch between dynamic behaviours is integral to overcoming modern challenges such as eradication of infectious diseases, the design of efficient legged robots and the protection of cyber infrastructures.

Regularized regression problems are ubiquitous in statistical modeling, signal processing, and machine learning. Sparse regression in particular has been instrumental in scientific model discovery, including compressed sensing applications, variable selection, and high-dimensional analysis. We propose a new and highly effective approach for regularized regression, called SR3. The key idea is to solve a relaxation of the regularized problem, which has three advantages over the state-of-the-art: (1) solutions of the relaxed problem are superior with respect to errors, false positives, and conditioning, (2) relaxation allows extremely fast algorithms for both convex and nonconvex formulations, and (3) the methods apply to composite regularizers such as total variation (TV) and its nonconvex variants. We demonstrate the improved performance of SR3 across a range of regularized regression problems with synthetic and real data, including compressed sensing, LASSO, matrix completion and TV regularization. To promote reproducible research, we include a companion MATLAB package that implements these popular applications.

The modified biharmonic equation is encountered in a variety of application areas, including streamfunction formulations of the Navier-Stokes equations. We develop a separation of variables representation for this equation in polar coordinates, for either the interior or exterior of a disk, and derive a new class of special functions which makes the approach stable. We discuss how these functions can be used in conjunction with fast algorithms to accelerate the solution of the modified biharmonic equation or the bi-Helmholtz equation in more complex geometries.

We present a novel integral representation for the biharmonic Dirichlet problem. To obtain the representation, the Dirichlet problem is first converted into a related Stokes problem for which the Sherman–Lauricella integral representation can be used. Not all potentials for the Dirichlet problem correspond to a potential for Stokes flow, and vice-versa, but we show that the integral representation can be augmented and modified to handle either simply or multiply connected domains. The resulting integral representation has a kernel which behaves better on domains with high curvature than existing representations. Thus, this representation results in more robust computational methods for the solution of the Dirichlet problem of the biharmonic equation and we demonstrate this with several numerical examples.

We consider a new class of periodic solutions to the Lugiato-Lefever equations (LLE) that govern the electromagnetic field in a microresonator cavity. Specifically, we rigorously characterize the stability and dynamics of the Jacobi elliptic function solutions of LLE and show that the dn solution is stabilized by the pumping of the microresonator. In analogy with soliton perturbation theory, we also derive a microcomb perturbation theory that allows one to consider the effects of physically realizable perturbations on the comb line stability, including effects of Raman scattering and stimulated emission. Our results are verified through full numerical simulations of the LLE cavity dynamics. The perturbation theory gives a simple analytic platform for potentially engineering new resonator designs.

The dynamic mode decomposition (DMD) has become a leading tool for data-driven modeling of dynamical systems, providing a regression framework for fitting linear dynamical models to time-series measurement data. We present a simple algorithm for computing an optimized version of the DMD for data which may be collected at unevenly spaced sample times. By making use of the variable projection method for nonlinear least squares problems, the algorithm is capable of solving the underlying nonlinear optimization problem efficiently. We explore the performance of the algorithm with some numerical examples for synthetic and real data from dynamical systems and find that the resulting decomposition displays less bias in the presence of noise than standard DMD algorithms. Because of the flexibility of the algorithm, we also present some interesting new options for DMD-based analysis.