We develop a new moving-water equilibria preserving numerical scheme for the Saint-Venant system. The new scheme is designed in two major steps. First, the geometric source term is incorporated into the discharge flux, which results in a hyperbolic system with a global flux. Second, the discharge equation is relaxed so that the nonlinearity is moved into the stiff right-hand side of the added auxiliary equation. The main advantages of the new scheme are that (i) no special treatment of the geometric source term is required, and (ii) no nonlinear (cubic) equations should be solved to obtain the point values of the water depth out of the reconstructed equilibrium variables, as it must be done in the existing alternative methods. We also develop a hybrid numerical flux, which helps to handle various flow regimes in a stable manner. Several numerical experiments are performed to verify that the proposed scheme is capable of exactly preserving general moving-water steady states and accurately capturing their small perturbations.

In this paper, we design and analyze third order positivity-preserving discontinuous Galerkin (DG) schemes for solving the time-dependent system of Poisson–Nernst–Planck (PNP) equations, which have found much use in diverse applications. Our DG method with Euler forward time discretization is shown to preserve the positivity of cell averages at all time steps. The positivity of numerical solutions is then restored by a scaling limiter in reference to positive weighted cell averages. The method is also shown to preserve steady states. Numerical examples are presented to demonstrate the third order accuracy and illustrate the positivity-preserving property in both one and two dimensions.

We present a regularized finite difference method for the logarithmic Schr ̈odingerequation (LogSE) and establish its error bound. Due to the blowup of the logarithmic nonlinearity, i.e., ln ρ→ −\infty when ρ→0+v with ρ = |u|^2 being the density and u being the complex-valuedwave function or order parameter, there are significant difficulties in designing numerical methodsand establishing their error bounds for the LogSE. In order to suppress the roundoff error and toavoid blowup, a regularized LogSE (RLogSE) is proposed with a small regularization parameter 0 < ε << 1 and linear convergence is established between the solutions of RLogSE and LogSE interm of ε. Then a semi-implicit finite difference method is presented for discretizing the RLogSEand error estimates are established in terms of the mesh sizehand time stepτas well as the smallregularization parameterε. Finally numerical results are reported to illustrate our error bounds.

We present two uniformly accurate numerical methods for discretizing the Zakharovsystem (ZS) with a dimensionless parameter 0< ε ≤ 1, which is inversely proportional to theacoustic speed. In the subsonic limit regime, i.e., 0< ε << 1, the solution of ZS propagates waves with O(ε)- andO(1)-wavelengths in time and space, respectively, and/or rapid outgoing initial layerswith speed O(1/ε) in space due to the singular perturbation of the wave operator in ZS and/or theincompatibility of the initial data. By adopting an asymptotic consistent formulation of ZS, wepresent a time-splitting exponential wave integrator (TS-EWI) method via applying a time-splittingtechnique and an exponential wave integrator for temporal derivatives in the nonlinear Schr ̈odingerequation and wave-type equation, respectively. By introducing a multiscale decomposition of ZS, wepropose a time-splitting multiscale time integrator (TS-MTI) method. Both methods are explicitand convergent exponentially in space for all kinds of initial data, which is uniformly for ε ∈ (0,1].The TS-EWI method is simpler to be implemented and it is only uniformly and optimally accuratein time for well-prepared initial data, while the TS-MTI method is uniformly and optimally accuratein time for any kind of initial data. Extensive numerical results are reported to show their efficiencyand accuracy, especially in the subsonic limit regime. Finally, the TS-MTI method is applied tostudy numerically convergence rates of ZS to its limiting models when ε→0+.

We propose an explicit numerical method for the periodic Korteweg–de Vries equation. Our method is based on a Lawson-type exponential integrator for time integration and the Rusanov scheme for Burgers’ nonlinearity. We prove first-order convergence in both space and time under a mild Courant–Friedrichs–Lewy condition τ=O(h)⁠, where τ and h represent the time step and mesh size for solutions in the Sobolev space H^3((−π,π))⁠, respectively. Numerical examples illustrating our convergence result are given.

We establish uniform error bounds of a finite difference method for the Klein-Gordon-Zakharov system (KGZ) with a dimensionless parameter ε∈(0,1], which is inversely proportional to the acoustic speed. In the subsonic limit regime, i.e. 0<ε≪1, the solution propagates highly oscillatory waves in time and/or rapid outgoing initial layers in space due to the singular perturbation in the Zakharov equation and/or the incompatibility of the initial data. Specifically, the solution propagates waves with O(ε)-wavelength in time and O(1)-wavelength in space as well as outgoing initial layers in space at speed O(1/ε). This high oscillation in time and rapid outgoing waves in space of the solution cause significant burdens in designing numerical methods and establishing error estimates for KGZ. By adapting an asymptotic consistent formulation, we propose a uniformly accurate finite difference method and rigorously establish two independent error bounds at O(h^2+τ^2/ε) and O(h^2+τ+ε) with h mesh size and τ time step. Thus we obtain a uniform error bound at O(h^2+τ) for 0<ε≤1. The main techniques in the analysis include the energy method, cut-off of the nonlinearity to bound the numerical solution, the integral approximation of the oscillatory term, and ε-dependent error bounds between the solutions of KGZ and its limiting model when ε→0+. Finally, numerical results are reported to confirm our error bounds.

We present a uniformly accurate finite difference method and establish rigorously its uniform error bounds for the Zakharov system (ZS) with a dimensionless parameter 0 < ε ≤ 1, which is inversely proportional to the speed of sound. In the subsonic limit regime, i.e., 0 < ε << 1, the solution propagates highly oscillatory waves and/or rapid outgoing initial layers due to the perturbation of the wave operator in the ZS and/or the incompatibility of the initial data, which is characterized by two parameters α ≥ 0 and β ≥ −1. Specifically, the solution propagates waves with wavelength of O(ε) and O(1) in time and space, respectively, and amplitude at O(ε^{min{2,α,1+β}}). This high oscillation of the solution in time brings significant difficulties in designing numerical methods and establishing their error bounds, especially in the subsonic limit regime. A uniformly accurate finite difference method is proposed by reformulating the ZS into an asymptotic consistent formulation and adopting an integral approximation of the oscillatory term. By applying the energy method and using the limiting equation via a nonlinear Schr ̈odinger equation with an oscillatory potential, we rigorously establish two independent error bounds at O(h^2 + τ^2/ε) and O(h^2 + τ^2 + τε^{α^∗} + ε^{1+α^∗}), respectively, with h the mesh size, τ the time step and α^∗ = min{1, α, 1 + β}. Thus we obtain error bounds at O(h^2 + τ^{4/3}) and O(h^2 + τ^{1+α^∗/(2+α^∗)}) for well-prepared (α^∗ = 1) and ill-prepared (0 ≤ α^∗ < 1) initial data, respectively, which are uniform in both space and time for 0 < ε ≤ 1 and optimal at the second order in space. Other techniques in the analysis include the cut-off technique for treating the nonlinearity and inverse estimates to bound the numerical solution. Numerical results are reported to demonstrate our error bounds.

We consider the Schrodinger equation with a logarithmic nonlinearity and a repulsive harmonic potential. Depending on the parameters of the equation, the solution may or may not be dispersive. When dispersion occurs, it does with an exponential rate in time. To control this, we change the unknown function through a generalized lens transform. This approach neutralizes the possible boundary effects, and could be used in the case of the Schrodinger equation without potential. We then employ standard splitting methods on the new equation via a nonuniform grid, after the logarithmic nonlinearity has been regularized. We also discuss the case of a power nonlinearity and give some results concerning the error estimates of the first-order Lie-Trotter splitting method for both cases of nonlinearities. Finally extensive numerical experiments are reported to investigate the dynamics of the equations.

We establish uniform error bounds of time-splitting Fourier pseudospectral (TSFP) methods for the nonlinear Klein–Gordon equation (NKGE) with weak power-type nonlinearity and O(1) initial data, while the nonlinearity strength is characterized by εp with a constant p ∈ N+ and a dimensionless parameter ε ∈ (0, 1], for the long-time dynamics up to the time at O(ε^{−β}) with 0 ≤ β ≤ p. In fact, when 0 < ε ≪ 1, the problem is equivalent to the long-time dynamics of NKGE with small initial data and O(1) nonlinearity strength, while the amplitude of the initial data (and the solution) is at O(ε). By reformulating the NKGE into a relativistic nonlinear Schr ̈odinger equation, we adapt the TSFP method to discretize it numerically. By using the method of mathematical induction to bound the numerical solution, we prove uniform error bounds at O(h^m + ε^{p−β} τ^2) of the TSFP method with h mesh size, τ time step and m ≥ 2 depending on the regularity of the solution. The error bounds are uniformly accurate for the long-time simulation up to the time at O(ε^{−β}) and uniformly valid for ε ∈ (0, 1]. Especially, the error bounds are uniformly at the second order rate for the large time step τ = O(ε^{−(p−β)/2}) in the parameter regime 0 ≤ β < p. Numerical results are reported to confirm our error bounds in the long-time regime. Finally, the TSFP method and its error bounds are extended to a highly oscillatory complex NKGE which propagates waves with wavelength at O(1) in space and O(ε^β) in time and wave velocity at O(ε^{−β}).

The logarithmic nonlinearity has been used in many partial differential equations (PDEs) for modeling problems in various applications. Due to the singularity of the logarithmic function, it introduces tremendous difficulties in establishing mathematical theories, as well as in designing and analyzing numerical methods for PDEs with such nonlinearity. Here, we take the logarithmic Schrödinger equation (LogSE) as a prototype model. Instead of regularizing f(ρ)=lnρ in the LogSE directly and globally as being done in the literature, we propose a local energy regularization (LER) for the LogSE by first regularizing F(ρ)=ρlnρ−ρ locally near ρ=0+ with a polynomial approximation in the energy functional of the LogSE and then obtaining an energy regularized logarithmic Schrödinger equation (ERLogSE) via energy variation. Linear convergence is established between the solutions of ERLogSE and LogSE in terms of a small regularization parameter 0<ε≪1. Moreover, the conserved energy of the ERLogSE converges to that of LogSE quadratically, which significantly improves the linear convergence rate of the regularization method in the literature. Error estimates are also presented for solving the ERLogSE by using Lie–Trotter splitting integrators. Numerical results are reported to confirm our error estimates of the LER and of the time-splitting integrators for the ERLogSE. Finally, our results suggest that the LER performs better than regularizing the logarithmic nonlinearity in the LogSE directly.

This paper is devoted to developing and analysing a highly accurate conservative method for solving the quantum Zakharov system. The scheme is based on a linearly-implicit compact finite difference discretization and conserve the mass as well as energy in discrete level. Detailed numerical analysis is presented which shows the method is fourth-order accurate in space and second-order accurate in time. Several numerical examples are reported to confirm the conservation properties and high accuracy of the proposed scheme. Finally the compact scheme is applied to study the convergence rate of the quantum Zakharov system to its limiting model in the semi-classical limit.

A Deuflhard-type exponential integrator sine pseudospectral (DEI-SP) method is proposed and analyzed for solving the generalized improved Boussinesq (GIBq) equation. The numerical scheme is based on a second-order exponential integrator for time integration and a sine pseudospectral discretization in space. Rigorous analysis and abundant experiments show that the method converges quadratically and spectrally in time and space, respectively. Finally the DEI-SP method is applied to investigate the complicated and interesting long-time dynamics of the GIBq equation.

We present a uniformly first order accurate numerical method for solving the Klein-Gordon-Zakharov (KGZ) system with two dimensionless parameters 0 < ε ≤ 1 and 0 <γ ≤ 1, which are inversely proportional to the plasma frequency and the acoustic speed, respectively. In the simultaneous high-plasma-frequency and subsonic limit regime, i.e. ε < γ → 0+, the KGZ system collapses to a cubic Schrödinger equation, and the solution propagates waves with O ( ε2)-wavelength in time and meanwhile contains rapid outgoing initial layers with speed O (1/ γ ) in space due to the incompatibility of the initial data. By presenting a multiscale decomposition of the KGZ system, we propose a multiscale time integrator Fourier pseudospectral method which is explicit, efficient and uniformly accurate for solving the KGZ system for all 0 < ε < γ ≤ 1. Numerical results are reported to show the efficiency and accuracy of scheme. Finally, the method is applied to investigate the convergence rates of the KGZ system to its limiting models when ε < γ → 0+.

Two fundamental difficulties are encountered in the numerical evaluation of time- dependent layer potentials. One is the quadratic cost of history dependence, which has been successfully addressed by splitting the potentials into two parts—a local part that contains the most recent contributions and a history part that contains the con- tributions from all earlier times. The history part is smooth, easily discretized using high-order quadratures, and straightforward to compute using a variety of fast algo- rithms. The local part, however, involves complicated singularities in the underlying Green’s function. Existing methods, based on exchanging the order of integration in space and time, are able to achieve high-order accuracy, but are limited to the case of stationary boundaries. Here, we present a new quadrature method that leaves the order of integration unchanged, making use of a change of variables that converts the singular integrals with respect to time into smooth ones. We have also derived asymptotic formulas for the local part that lead to fast and accurate hybrid schemes, extending earlier work for scalar heat potentials and applicable to moving boundaries. The performance of the overall scheme is demonstrated via numerical examples.

With rheology applications in mind, we present a fast solver for the time-dependent effective viscosity of an infinite lattice containing one or more neutrally buoyant smooth rigid particles per unit cell, in a two-dimensional Stokes fluid with given shear rate. At each time, the mobility problem is reformulated as a 2nd-kind boundary integral equation, then discretized to spectral accuracy by the Nyström method and solved iteratively, giving typically 10 digits of accuracy. Its solution controls the evolution of particle locations and angles in a first-order system of ordinary differential equations. The formulation is placed on a rigorous footing by defining a generalized periodic Green’s function for the skew lattice. Numerically, the periodized integral operator is split into a near image sum— applied in linear time via the fast multipole method—plus a correction field solved cheaply via proxy Stokeslets. We use barycentric quadratures to evaluate particle interactions and velocity fields accurately, even at distances much closer than the node spacing. Using first- order time-stepping we simulate, for example, 25 ellipses per unit cell to 3-digit accuracy on a desktop in 1 hour per shear time. Our examples show equilibration at long times, force chains, and two types of blow-ups (jamming) whose power laws match lubrication theory asymptotics.

We study the stability properties of explicit marching schemes for second-kind Volterra integral equations that arise when solving boundary value problems for the heat equation by means of potential theory. It is well known that explicit finite-difference or finite-element schemes for the heat equation are stable only if the time step 􏰃t is of the order O(􏰃x2), where 􏰃x is the finest spatial grid spacing. In contrast, for the Dirichlet and Neumann problems on the unit ball in all dimensions d ≥ 1, we show that the simplest Volterra marching scheme, i.e., the forward Euler scheme, is unconditionally stable. Our proof is based on an explicit spectral radius bound of the marching matrix, leading to an estimate that an L2-norm of the solution to the integral equation is bounded by cd T d /2 times the norm of the right-hand side. For the Robin problem on the half-space in any dimension, with constant Robin (heat transfer) coefficient κ, we exhibit a constant C such that the forward Euler scheme is stable if 􏰃t < C /κ 2 , independent of any spatial discretization. This relies on new lower bounds on the spectrum of real symmetric Toeplitz matrices defined by convex sequences. Finally, we show that the forward Euler scheme is unconditionally stable for the Dirichlet problem on any smooth convex domain in any dimension, in the L∞-norm.

We present a hybrid asymptotic/numerical method for the accurate computation of single- and double-layer heat potentials in two dimensions. It has been shown in previous work that simple quadrature schemes suffer from a phenomenon called “geometrically induced stiffness,” meaning that formally high-order accurate methods require excessively small time steps before the rapid convergence rate is observed. This can be overcome by analytic integration in time, requiring the evaluation of a collection of spatial boundary integral operators with non-physical, weakly singular kernels. In our hybrid scheme, we combine a local asymptotic approximation with the evaluation of a few boundary integral operators involving only Gaussian kernels, which are easily accelerated by a new version of the fast Gauss transform. This new scheme is robust, avoids geometrically induced stiffness, and is easy to use in the presence of moving geometries. Its extension to three dimensions is natural and straightforward, and should permit layer heat potentials to become flexible and powerful tools for modeling diffusion processes.

A variety of problems in computational physics and engineering require the convolu- tion of the heat kernel (a Gaussian) with either discrete sources, densities supported on boundaries, or continuous volume distributions. We present a unified fast Gauss transform for this purpose in two dimensions, making use of an adaptive quad-tree discretization on a unit square which is assumed to contain all sources. Our implementation permits either free-space or periodic boundary conditions to be imposed, and is efficient for any choice of variance in the Gaussian.

In this paper, we propose an unfitted Nitsche's method for computing wave modes in topological materials. The proposed method is based on Nitsche's technique to study the performance-enhanced topological materials which have strongly heterogeneous structures (e.g., the refractive index is piecewise constant with high contrasts). For periodic bulk materials, we use Floquet-Bloch theory and solve an eigenvalue problem on a torus with unfitted meshes. For the materials with a line defect, a sufficiently large domain with zero boundary conditions is used to compute the localized eigenfunctions corresponding to the edge modes. The interfaces are handled by Nitsche's method on an unfitted uniform mesh. We prove the proposed methods converge optimally, and present numerical examples to validate the theoretical results and demonstrate the capability of simulating topological materials.

In this paper, we propose an unfitted Nitsche's method to compute the band structures of phononic crystal with impurities of general geometry. The proposed method does not require the background mesh to fit the interfaces of impurities, and thus avoids the expensive cost of generating body-fitted meshes and simplifies the inclusion of interface conditions in the formulation. The quasi-periodic boundary conditions are handled by the Floquet-Bloch transform, which converts the computation of band structures into an eigenvalue problem with periodic boundary conditions. More importantly, we show the well-posedness of the proposed method using a delicate argument based on the trace inequality, and further prove the convergence by the Babuška-Osborn theory. We achieve the optimal convergence rate at the presence of the impurities of general geometry. We confirm the theoretical results by two numerical examples, and show the capability of the proposed methods for computing the band structures without fitting the interfaces of impurities.

Photonic graphene, a photonic crystal with honeycomb structures, has been intensively studied in both theoretical and applied fields. Similar to graphene which admits Dirac Fermions and topological edge states, photonic graphene supports novel and subtle propagating modes (edge modes) of electromagnetic waves. These modes have wide applications in many optical systems. In this paper, we propose a novel gradient recovery method based on Bloch theory for the computation of topological edge modes in the honeycomb structure. Compared to standard finite element methods, this method provides higher order accuracy with the help of gradient recovery technique. This high order accuracy is highly desired for constructing the propagating electromagnetic modes in applications. We analyze the accuracy and prove the superconvergence of this method. Numerical examples are presented to show the efficiency by computing the edge mode for the P-symmetry and C-symmetry breaking cases in honeycomb structures.

No abstract uploaded!

A key difficulty in the analysis and numerical approximation of the shallow water equations is the nonconservative product of measures due to the gravitational force acting on a sloped bottom. Solutions may be nonunique, and numerical schemes are not only consistent discretizations of the shallow water equations, but they also determine how to model the physics. Our derivation is based on a subcell reconstruction using infinitesimal singular layers at the cell boundaries, as inspired by S. Noelle, Y. Xing, and C.-W. Shu [J. Comput. Phys., 226 (2007), pp. 29-58]. One key step is to separate the singular measures. Another aspect is the reconstruction of the solution variables in the singular layers. We study three reconstructions. The first leads to the well-known scheme of Audusse et al., [SIAM J. Sci. Comput., 25 (2004), pp. 2050-2065], which introduces the hydrostatic reconstruction. The second is a modification proposed in [T. Morales de Luna, M. J. Castro D´ıaz, and C. Par´es, Appl. Math. Comput., 219 (2013), pp. 9012-9032], which analyzes whether a wave has enough energy to overcome a step. The third is our new scheme, which borrows its structure from the wet-dry front. For a number of cases discussed in recent years, where water runs down a hill, Audusse’s scheme converges slowly or fails. Morales’ scheme gives a visible improvement. Both schemes are clearly outperformed by our new scheme.

A novel discontinuous Galerkin (DG) method is developed to solve timedependent bi-harmonic type equations involving fourth derivatives in one and multiple space dimensions. We present the spatial DG discretization based on a mixed formulation and central interface numerical fluxes so that the resulting semi-discrete schemes are L^2 stable even without interior penalty. For time discretization, we use Crank–Nicolson so that the resulting scheme is unconditionally stable and second order in time. We present the optimal L^2 error estimate of O(h^{k+1}) for polynomials of degree k for semi-discrete DG schemes, and the L2 error of O(h^{k+1} + (\Delta t)^2) for fully discrete DG schemes. Extensions to more general fourth order partial differential equations and cases with non-homogeneous boundary conditions are provided. Numerical results are presented to verify the stability and accuracy of the schemes. Finally, an application to the one-dimensional Swift–Hohenberg equation endowed with a decay free energy is presented.

No abstract uploaded!