Search

In this paper we prove that for stable semi-discretizations of the wave equation for the WaveHoltz iteration is guaranteed to converge to an approximate solution of the corresponding frequency domain problem, if it exists. We show that for certain classes of frequency domain problems, the WaveHoltz iteration without acceleration converges in $O({\omega})$ iterations with the constant factor depending logarithmically on the desired tolerance. We conjecture that the Helmholtz problem in open domains with no trapping waves is one such class of problems and we provide numerical examples in one and two dimensions using finite differences and discontinuous Galerkin discretizations which demonstrate these converge results.

The implicit boundary integral method (IBIM) provides a framework to construct quadrature rules on regular lattices for integrals over irregular domain boundaries. This work provides a systematic error analysis for IBIMs on uniform Cartesian grids for boundaries with different degree of regularities. We first show that the quadrature error gains an addition order of $\frac{d-1}{2}$ from the curvature for a strongly convex smooth boundary due to the ``randomness'' in the signed distances. This gain is discounted for degenerated convex surfaces. We then extend the error estimate to general boundaries under some special circumstances, including how quadrature error depends on the boundary's local geometry relative to the underlying grid. Bounds on the variance of the quadrature error under random shifts and rotations of the lattices are also derived.

In this work we show an error estimate for a first order Gaussian beam at a fold caustic, approximating time-harmonic waves governed by the Helmholtz equation. For the caustic that we study the exact solution can be constructed using Airy functions and there are explicit formulae for the Gaussian beam parameters. Via precise comparisons we show that the pointwise error on the caustic is of the order $O(k^{-5/6})$ where $k$ is the wave number in Helmholtz.

In this paper, we solve the linearized Poisson-Boltzmann equation, used to model the electric potential of macromolecules in a solvent. We derive a corrected trapezoidal rule with improved accuracy for a boundary integral formulation of the linearized Poisson-Boltzmann equation. More specifically, in contrast to the typical boundary integral formulations, the corrected trapezoidal rule is applied to integrate a system of compacted supported singular integrals using uniform Cartesian grids in $\mathbb{R}^3$, without explicit surface parameterization. A Krylov method, accelerated by a fast multipole method, is used to invert the resulting linear system. We study the efficacy of the proposed method, and compare it to an existing, lower order method. We then apply the method to the computation of electrostatic potential of macromolecules immersed in solvent. The solvent excluded surfaces, defined by a common approach, are merely piecewise smooth, and we study the effectiveness of the method for such surfaces., Comment: 22 pages, 6 figures

We present convergence theory for corrected quadrature rules on uniform Cartesian grids for functions with a point singularity. We begin by deriving an error estimate for the punctured trapezoidal rule, and then derive error expansions. We define the corrected trapezoidal rules, based on the punctured trapezoidal rule, where the weights for the nodes close to the singularity are judiciously corrected based on these expansions. Then we define the composite corrected trapezoidal rules for a larger family of functions using series expansions around the point singularity and applying corrected trapezoidal rules appropriately. We prove that we can achieve high order accuracy by using a sufficient number of correction nodes around the point singularity and of expansion terms., Comment: 23 pages, 1 figure

In this paper we extend analysis of the WaveHoltz iteration -- a time-domain iterative method for the solution of the Helmholtz equation. We expand the previous analysis of energy conserving problems and prove convergence of the WaveHoltz iteration for problems with impedance boundary conditions in a single spatial dimension. We then consider interior Dirichlet/Neumann problems with damping in any spatial dimension, and show that for a sufficient level of damping the WaveHoltz iteration converges in a number of iteration independent of the frequency. Finally, we present a discrete analysis of the WaveHoltz iteration for a family of higher order time-stepping schemes. We show that the fixed-point of the discrete WaveHoltz iteration converges to the discrete Helmholtz solution with the order of the time-stepper chosen. We present numerical examples and demonstrate that it is possible to completely remove time discretization error from the WaveHoltz solution through careful analysis of the discrete iteration together with updated quadrature formulas.

We consider the application of the WaveHoltz iteration to time-harmonic elastic wave equations with energy conserving boundary conditions. The original WaveHoltz iteration for acoustic Helmholtz problems is a fixed-point iteration that filters the solution of the wave equation with time-harmonic forcing and boundary data. As in the original WaveHoltz method, we reformulate the fixed point iteration as a positive definite linear system of equations that is iteratively solved by a Krylov method. We present two time-stepping schemes, one explicit and one (novel) implicit, which completely remove time discretization error from the WaveHoltz solution by performing a simple modification of the initial data and time-stepping scheme. Numerical experiments indicate an iteration scaling similar to that of the original WaveHoltz method, and that the convergence rate is dictated by the shortest (shear) wave speed of the problem. We additionally show that the implicit scheme can be advantageous in practice for meshes with disparate element sizes.

We present a family of high order trapezoidal rule-based quadratures for a class of singular integrals, where the integrand has a point singularity. The singular part of the integrand is expanded in a Taylor series involving terms of increasing smoothness. The quadratures are based on the trapezoidal rule, with the quadrature weights for Cartesian nodes close to the singularity judiciously corrected based on the expansion. High order accuracy can be achieved by utilizing a sufficient number of correction nodes around the singularity to approximate the terms in the series expansion. The derived quadratures are applied to the Implicit Boundary Integral formulation of surface integrals involving the Laplace layer kernels., Comment: 41 pages, 9 figures, 2 tables, 2 appendices

In this paper, we present a finite difference heterogeneous multiscale method for the Landau-Lifshitz equation with a highly oscillatory diffusion coefficient. The approach combines a higher order discretization and artificial damping in the so-called micro problem to obtain an efficient implementation. The influence of different parameters on the resulting approximation error is discussed. Numerical examples for both periodic as well as more general coefficients are given to demonstrate the functionality of the approach.

We present new higher-order quadratures for a family of boundary integral operators re-derived using the approach introduced in [Kublik, Tanushev, and Tsai - J. Comp. Phys. 247: 279-311, 2013]. In this formulation, a boundary integral over a smooth, closed hypersurface is transformed into an equivalent volume integral defined in a sufficiently thin tubular neighborhood of the surface. The volumetric formulation makes it possible to use the simple trapezoidal rule on uniform Cartesian grids and relieves the need to use parameterization for developing quadrature. Consequently, typical point singularities in a layer potential extend along the surface's normal lines. We propose new higher-order corrections to the trapezoidal rule on the grid nodes around the singularities. This correction is based on local decompositions of the singularity and is dependent on the angle of approach to the singularity relative to the surface's principal curvature directions. The proposed decomposition, combined with the volumetric formulation, leads to a special quadrature error cancellation., Comment: 44 pages, 13 figures, 2 tables

In this paper, we consider several possible ways to set up Heterogeneous Multiscale Methods for the Landau-Lifshitz equation with a highly oscillatory diffusion coefficient, which can be seen as a means to modeling rapidly varying ferromagnetic materials. We then prove estimates for the errors introduced when approximating the relevant quantity in each of the models given a periodic problem, using averaging in time and space of the solution to a corresponding micro problem. In our setup, the Landau-Lifshitz equation with highly oscillatory coefficient is chosen as the micro problem for all models. We then show that the averaging errors only depend on $\varepsilon$, the size of the microscopic oscillations, as well as the size of the averaging domain in time and space and the choice of averaging kernels.

In this paper, we consider homogenization of the Landau-Lifshitz equation with a highly oscillatory material coefficient with period $\varepsilon$ modeling a ferromagnetic composite. We derive equations for the homogenized solution to the problem and the corresponding correctors and obtain estimates for the difference between the exact and homogenized solution as well as corrected approximations to the solution. Convergence rates in $\varepsilon$ over times $O(\varepsilon^\sigma)$ with $0\leq \sigma\leq 2$ are given in the Sobolev norm $H^q$, where $q$ is limited by the regularity of the solution to the detailed Landau-Lifshitz equation and the homogenized equation. The rates depend on $q$, $\sigma$ and the the number of correctors.

A new idea for iterative solution of the Helmholtz equation is presented. We show that the iteration which we denote WaveHoltz and which filters the solution to the wave equation with harmonic data evolved over one period, corresponds to a coercive operator or a positive definite matrix in the discretized case.

We consider the wave equation with uncertain initial data and medium, when the wavelength $\varepsilon$ of the solution is short compared to the distance traveled by the wave. We are interested in the statistics for quantities of interest (QoI), defined as functionals of the wave solution, given the probability distributions of the uncertain parameters in the wave equation. Fast methods to compute this statistics require considerable smoothness in the mapping from parameters to the QoI, which is typically not present in the high frequency case, as the oscillations on the $\varepsilon$ scale in the wave field is inherited by the QoIs. The main contribution of this work is to identify certain non-oscillatory quadratic QoIs and show $\varepsilon$-independent estimates for the derivatives of the QoI with respect to the parameters, when the wave solution is replaced by a Gaussian beam approximation.

This paper concerns the analysis of a multiscale method for wave propagation problems in microscopically nonhomogeneous media. A direct numerical approximation of such problems is prohibitively expensive as it requires resolving the microscopic variations over a much larger physical domain of interest. The heterogeneous multiscale method (HMM) is an efficient framework to approximate the solutions of multiscale problems. In HMM, one assumes an incomplete macroscopic model which is coupled to a known but expensive microscopic model. The micromodel is solved only locally to upscale the parameter values which are missing in the macromodel. The resulting macroscopic model can then be solved at a cost independent of the small scales in the problem. In general, the accuracy of the HMM is related to how good the upscaling step approximates the right macroscopic quantities. The analysis of the method, that we consider here, was previously addressed only in purely periodic media although the method itself is numerically shown to be applicable to more general settings. In the present study, we consider a more realistic setting by assuming a locally-periodic medium where slow and fast variations are allowed at the same time. We then prove that HMM captures the right macroscopic effects. The generality of the tools and ideas in the analysis allows us to establish convergence rates in a multi-dimensional setting. The theoretical findings here imply an improved convergence rate in one-dimension, which also justifies the numerical observations from our earlier study., Comment: 29 pages

This work is concerned with the accuracy of Gaussian beam superpositions, which are asymptotically valid high frequency solutions to linear hyperbolic partial differential equations and the Schr\"odinger equation. We derive Sobolev and max norms estimates for the difference between an exact solution and the corresponding Gaussian beam approximation, in terms of the short wavelength $\varepsilon$. The estimates are performed for the scalar wave equation and the Schr\"odinger equation. Our result demonstrates that a Gaussian beam superposition with $k$-th order beams converges to the exact solution as $O(\varepsilon^{k/2-s})$ in order $s$ Sobolev norms. This result is valid in any number of spatial dimensions and it is unaffected by the presence of caustics in the solution. In max norm, we show that away from caustics the convergence rate is $O(\varepsilon^{\lceil k/2\rceil})$ and away from the essential support of the solution, the convergence is spectral in $\varepsilon$. However, in the neighborhood of a caustic point we are only able to show the slower, and dimensional dependent, rate $O(\varepsilon^{(k-n)/2})$ in $n$ spatial dimensions.

We consider the wave equation with highly oscillatory initial data, where there is uncertainty in the wave speed, initial phase and/or initial amplitude. To estimate quantities of interest related to the solution and their statistics, we combine a high-frequency method based on Gaussian beams with sparse stochastic collocation. Although the wave solution, $u^\varepsilon$, is highly oscillatory in both physical and stochastic spaces, we provide theoretical arguments and numerical evidence that quantities of interest based on local averages of $|u^\varepsilon|^2$ are smooth, with derivatives in the stochastic space uniformly bounded in $\varepsilon$, where $\varepsilon$ denotes the short wavelength. This observable related regularity makes the sparse stochastic collocation approach more efficient than Monte Carlo methods. We present numerical tests that demonstrate this advantage.

We present a new computational approach for a class of large-scale nonlinear eigenvalue problems (NEPs) that are nonlinear in the eigenvalue. The contribution of this paper is two-fold. We derive a new iterative algorithm for NEPs, the tensor infinite Arnoldi method (TIAR), which is applicable to a general class of NEPs, and we show how to specialize the algorithm to a specific NEP: the waveguide eigenvalue problem. The waveguide eigenvalue problem arises from a finite-element discretization of a partial differential equation (PDE) used in the study waves propagating in a periodic medium. The algorithm is successfully applied to accurately solve benchmark problems as well as complicated waveguides. We study the complexity of the specialized algorithm with respect to the number of iterations m and the size of the problem n, both from a theoretical perspective and in practice. For the waveguide eigenvalue problem, we establish that the computationally dominating part of the algorithm has complexity O(nm^2 + sqrt(n) m^3). Hence, the asymptotic complexity of TIAR applied to the waveguide eigenvalue problem, for n that goes to infinity, is the same as for Arnoldi's method for standard eigenvalue problems.

In this work we construct Gaussian beam approximations to solutions of the high frequency Helmholtz equation with a localized source. Under the assumption of non-trapping rays we show error estimates between the exact outgoing solution and Gaussian beams in terms of the wave number $k$, both for single beams and superposition of beams. The main result is that the relative local $L^2$ error in the beam approximations decay as {$k^{-N/2}$ independent of dimension and presence of caustics, for $N$-th order beams.

Multiscale problems are computationally costly to solve by direct simulation because the smallest scales must be represented over a domain determined by the largest scales of the problem. We have developed and analyzed new numerical methods for multiscale wave propagation following the framework of the heterogeneous multiscale method. The numerical methods couple simulations on macro- and microscales for problems with rapidly fluctuating material coefficients. The computational complexity of the new method is significantly lower than that of traditional techniques. We focus on HMM approximation applied to long time integration of one-dimensional wave propagation problems in both periodic and non-periodic medium and show that the dispersive effect that appear after long time is fully captured.

Gaussian beams are asymptotically valid high frequency solutions to hyperbolic partial differential equations, concentrated on a single curve through the physical domain. They can also be extended to some dispersive wave equations, such as the Schr\"odinger equation. Superpositions of Gaussian beams provide a powerful tool to generate more general high frequency solutions that are not necessarily concentrated on a single curve. This work is concerned with the accuracy of Gaussian beam superpositions in terms of the wavelength $\epsilon$. We present a systematic construction of Gaussian beam superpositions for all strictly hyperbolic and Schr\"odinger equations subject to highly oscillatory initial data of the form $Ae^{i\Phi/\epsilon}$. Through a careful estimate of an oscillatory integral operator, we prove that the $k$-th order Gaussian beam superposition converges to the original wave field at a rate proportional to $\epsilon^{k/2}$ in the appropriate norm dictated by the well-posedness estimate. In particular, we prove that the Gaussian beam superposition converges at this rate for the acoustic wave equation in the standard, $\epsilon$-scaled, energy norm and for the Schr\"odinger equation in the $L^2$ norm. The obtained results are valid for any number of spatial dimensions and are unaffected by the presence of caustics. We present a numerical study of convergence for the constant coefficient acoustic wave equation in $\Real^2$ to analyze the sharpness of the theoretical results.

Multi-scale wave propagation problems are computationally costly to solve by traditional techniques because the smallest scales must be represented over a domain determined by the largest scales of the problem. We have developed and analyzed new numerical methods for multi-scale wave propagation in the framework of heterogeneous multi-scale method. The numerical methods couples simulations on macro- and micro-scales for problems with rapidly oscillating coefficients. We show that the complexity of the new method is significantly lower than that of traditional techniques with a computational cost that is essentially independent of the micro-scale. A convergence proof is given and numerical results are presented for periodic problems in one, two and three dimensions. The method is also successfully applied to non-periodic problems and for long time integration where dispersive effects occur.

We expose in full detail a constructive procedure to invert the so--called "finite Markov moment problem". The proofs rely on the general theory of Toeplitz matrices together with the classical Newton's relations.

The Gaussian beam superposition method is an asymptotic method for computing high frequency wave fields in smoothly varying inhomogeneous media. In this paper we study the accuracy of the Gaussian beam superposition method and derive error estimates related to the discretization of the superposition integral and the Taylor expansion of the phase and amplitude off the center of the beam. We show that in the case of odd order beams, the error is smaller than a simple analysis would indicate because of error cancellation effects between the beams. Since the cancellation happens only when odd order beams are used, there is no remarkable gain in using even order beams. Moreover, applying the error estimate to the problem with constant speed of propagation, we show that in this case the local beam width is not a good indicator of accuracy, and there is no direct relation between the error and the beam width. We present numerical examples to verify the error estimates.

We consider a class of finite Markov moment problems with arbitrary number of positive and negative branches. We show criteria for the existence and uniqueness of solutions, and we characterize in detail the non-unique solution families. Moreover, we present a constructive algorithm to solve the moment problems numerically and prove that the algorithm computes the right solution.

We present and discuss a framework for computer-aided multiscale analysis, which enables models at a "fine" (microscopic/stochastic) level of description to perform modeling tasks at a "coarse" (macroscopic, systems) level. These macroscopic modeling tasks, yielding information over long time and large space scales, are accomplished through appropriately initialized calls to the microscopic simulator for only short times and small spatial domains. Our equation-free (EF) approach, when successful, can bypass the derivation of the macroscopic evolution equations when these equations conceptually exist but are not available in closed form. We discuss how the mathematics-assisted development of a computational superstructure may enable alternative descriptions of the problem physics (e.g. Lattice Boltzmann (LB), kinetic Monte Carlo (KMC) or Molecular Dynamics (MD) microscopic simulators, executed over relatively short time and space scales) to perform systems level tasks (integration over relatively large time and space scales,"coarse" bifurcation analysis, optimization, and control) directly. In effect, the procedure constitutes a systems identification based, "closure on demand" computational toolkit, bridging microscopic/stochastic simulation with traditional continuum scientific computation and numerical analysis. We illustrate these ideas through examples from chemical kinetics (LB, KMC), rheology (Brownian Dynamics), homogenization and the computation of "coarsely self-similar" solutions, and discuss various features, limitations and potential extensions of the approach.