Search

Numerical methods of approximate solution of the Cauchy problem for coupled systems of evolution equations are considered. Separating simpler subproblems for individual components of the solution achieves simplification of the problem at a new level in time. The decoupling method, a significant approach to simplifying the problem, is based on the decomposition of the problem's operator matrix. The approximate solution is constructed based on the linear composition of solutions to auxiliary problems. The paper investigates decoupling variants based on extracting the diagonal part of the operator matrix and the lower and upper triangular submatrices. The study introduces a new decomposition approach, which involves splitting the operator matrix into rows and columns. The composition stage utilizes various variants of splitting schemes, showcasing the versatility of the approach. In additive operator-difference schemes, we can distinguish explicit-implicit schemes, factorized schemes for two-component splitting, and regularized schemes for general multicomponent splitting. The study of stability of two- and three-level decoupling composition schemes is carried out using the theory of stability (correctness) of operator-difference schemes for finite-dimensional Hilbert spaces. The theoretical results of the decoupling technique for systems of evolution equations are illustrated on a test two-dimensional problem for a coupled system of two diffusion equations with inhomogeneous self- and cross-diffusion coefficients., Comment: 28 pages, 10 figures

In this paper, we introduce a novel category of central compact schemes inspired by existing cell-node and cell-centered compact finite difference schemes, that offer a superior spectral resolution for solving the dispersive wave equation. In our approach, we leverage both the function values at the cell nodes and cell centers to calculate third-order spatial derivatives at the cell nodes. To compute spatial derivatives at the cell centers, we employ a technique that involves half-shifting the indices within the formula initially designed for the cell-nodes. In contrast to the conventional compact interpolation scheme, our proposed method effectively sidesteps the introduction of transfer errors. We employ the Taylor-series expansion-based method to calculate the finite difference coefficients. By conducting systematic Fourier analysis and numerical tests, we note that the methods exhibit exceptional characteristics such as high order, superior resolution, and low dissipation. Computational findings further illustrate the effectiveness of high-order compact schemes, particularly in addressing problems with a third derivative term., Comment: 28 pages, 15 figures, 10 Tables

Stable computational algorithms for the approximate solution of the Cauchy problem for nonstationary problems are based on implicit time approximations. Computational costs for boundary value problems for systems of coupled multidimensional equations can be reduced by additive decomposition of the problem operator(s) and composition of the approximate solution using particular explicit-implicit time approximations. Such a technique is currently applied in conditions where the decomposition step is uncomplicated. A general approach is proposed to construct decomposition-composition algorithms for evolution equations in finite-dimensional Hilbert spaces. It is based on two main variants of the decomposition of the unit operator in the corresponding spaces at the decomposition stage and the application of additive operator-difference schemes at the composition stage. The general results are illustrated on the boundary value problem for a second-order parabolic equation by constructing standard splitting schemes on spatial variables and region-additive schemes (domain decomposition schemes)., Comment: 22 pages, 3 figures

We provide a new theoretical framework for the variable-step deferred correction (DC) methods based on the well-known BDF2 formula. By using the discrete orthogonal convolution kernels, some high-order BDF2-DC methods are proven to be stable on arbitrary time grids according to the recent definition of stability (SINUM, 60: 2253-2272). It significantly relaxes the existing step-ratio restrictions for the BDF2-DC methods (BIT, 62: 1789-1822). The associated sharp error estimates are established by taking the numerical effects of the starting approximations into account, and they suggest that the BDF2-DC methods have no aftereffect, that is, the lower-order starting scheme for the BDF2 scheme will not cause a loss in the accuracy of the high-order BDF2-DC methods. Extensive tests on the graded and random time meshes are presented to support the new theory., Comment: 27 pages, 12 tables, 8 figures

In the present work, we demonstrate how the pseudoinverse concept from linear algebra can be used to represent and analyze the boundary conditions of linear systems of partial differential equations. This approach has theoretical and practical implications; the theory applies even if the boundary operator is rank deficient, or near rank deficient. If desired, the pseudoinverse can be implemented directly using standard tools like Matlab. We also introduce a new and simplified version of the semidiscrete approximation of the linear PDE system, which completely avoids taking the time derivative of the boundary data. The stability results are valid for general, nondiagonal summation-by-parts norms. Another key result is the extension of summation-by-parts operators to multi-domains by means of carefully crafted embedding operators. No extra numerical boundary conditions are required at the grid interfaces. The aforementioned pseudoinverse allows for a compact representation of these multi-block operators, which preserves all relevant properties of the single-block operators. The embedding operators can be constructed for multiple space dimensions. Numerical results for the two-dimensional Maxwell's equations are presented, and they show very good agreement with theory., Comment: 109 pages, 7 figures

For solving two-dimensional incompressible flow in the vorticity form by the fourth-order compact finite difference scheme and explicit strong stability preserving (SSP) temporal discretizations, we show that the simple bound-preserving limiter in [Li H., Xie S., Zhang X., SIAM J. Numer. Anal., 56 (2018)]. can enforce the strict bounds of the vorticity, if the velocity field satisfies a discrete divergence free constraint. For reducing oscillations, a modified TVB limiter adapted from [Cockburn B., Shu CW., SIAM J. Numer. Anal., 31 (1994)] is constructed without affecting the bound-preserving property. This bound-preserving finite difference method can be used for any passive convection equation with a divergence free velocity field.

We show that the classical fourth order accurate compact finite difference scheme with high order strong stability preserving time discretizations for convection diffusion problems satisfies a weak monotonicity property, which implies that a simple limiter can enforce the bound-preserving property without losing conservation and high order accuracy. Higher order accurate compact finite difference schemes satisfying the weak monotonicity will also be discussed.

Due to the lack of corresponding analysis on appropriate mapping operator between two grids, high-order two-grid difference algorithms are rarely studied. In this paper, we firstly discuss the boundedness of a local bi-cubic Lagrange interpolation operator. And then, taking the semilinear parabolic equation as an example, we first construct a variable-step high-order nonlinear difference algorithm using compact difference technique in space and the second-order backward differentiation formula (BDF2) with variable temporal stepsize in time. With the help of discrete orthogonal convolution (DOC) kernels and a cut-off numerical technique, the unique solvability and corresponding error estimates of the high-order nonlinear difference scheme are established under assumptions that the temporal stepsize ratio satisfies rk < 4.8645 and the maximum temporal stepsize satisfies tau = o(h^1/2 ). Then, an efficient two-grid high-order difference algorithm is developed by combining a small-scale variable-step high-order nonlinear difference algorithm on the coarse grid and a large-scale variable-step high-order linearized difference algorithm on the fine grid, in which the constructed piecewise bi-cubic Lagrange interpolation mapping operator is adopted to project the coarse-grid solution to the fine grid. Under the same temporal stepsize ratio restriction rk < 4.8645 and a weaker maximum temporal stepsize condition tau = o(H^1.2 ), optimal fourth-order in space and second-order in time error estimates of the two-grid difference scheme is established if the coarse-fine grid stepsizes satisfy H = O(h^4/7). Finally, several numerical experiments are carried out to demonstrate the effectiveness and efficiency of the proposed scheme., Comment: 34 pages, 6 figures

In this paper, we consider the variable-order time fractional mobile/immobile diffusion (TF-MID) equation in two-dimensional spatial domain, where the fractional order $\alpha(t)$ satisfies $0<\alpha_{*}\leq \alpha(t)\leq \alpha^{*}<1$. We combine the quadratic spline collocation (QSC) method and the $L1^+$ formula to propose a QSC-$L1^+$ scheme. It can be proved that, the QSC-$L1^+$ scheme is unconditionally stable and convergent with $\mathcal{O}(\tau^{\min{\{3-\alpha^*-\alpha(0),2\}}} + \Delta x^{2}+\Delta y^{2})$, where $\tau$, $\Delta x$ and $\Delta y$ are the temporal and spatial step sizes, respectively. With some proper assumptions on $\alpha(t)$, the QSC-$L1^+$ scheme has second temporal convergence order even on the uniform mesh, without any restrictions on the solution of the equation. We further construct a novel alternating direction implicit (ADI) framework to develop an ADI-QSC-$L1^+$ scheme, which has the same unconditionally stability and convergence orders. In addition, a fast implementation for the ADI-QSC-$L1^+$ scheme based on the exponential-sum-approximation (ESA) technique is proposed. Moreover, we also introduce the optimal QSC method to improve the spatial convergence to fourth-order. Numerical experiments are attached to support the theoretical analysis, and to demonstrate the effectiveness of the proposed schemes., Comment: 37 pages

A stencil-adaptive SBP-SAT finite difference scheme is shown to display superconvergent behavior. Applied to the linear advection equation, it has a convergence rate $\mathcal{O}(\Delta x^4)$ in contrast to a conventional scheme, which converges at a rate $\mathcal{O}(\Delta x^3)$., Comment: 9 pages, 1 figure

We present a novel numerical method for solving the anisotropic diffusion equation in magnetic fields confined to a periodic box which is accurate and provably stable. We derive energy estimates of the solution of the continuous initial boundary value problem. A discrete formulation is presented using operator splitting in time with the summation by parts finite difference approximation of spatial derivatives for the perpendicular diffusion operator. Weak penalty procedures are derived for implementing both boundary conditions and parallel diffusion operator obtained by field line tracing. We prove that the fully-discrete approximation is unconditionally stable. Discrete energy estimates are shown to match the continuous energy estimate given the correct choice of penalty parameters. A nonlinear penalty parameter is shown to provide an effective method for tuning the parallel diffusion penalty and significantly minimises rounding errors. Several numerical experiments, using manufactured solutions, the ``NIMROD benchmark'' problem and a single island problem, are presented to verify numerical accuracy, convergence, and asymptotic preserving properties of the method. Finally, we present a magnetic field with chaotic regions and islands and show the contours of the anisotropic diffusion equation reproduce key features in the field., Comment: 28 pages, 11 figures

We solve the anisotropic diffusion equation in 2D, where the dominant direction of diffusion is defined by a vector field which does not conform to a Cartesian grid. Our method uses operator splitting to separate the diffusion perpendicular and parallel to the vector field. The slow time scale is solved using a provably stable finite difference formulation in the perpendicular to the vector field, and an integral operator for the diffusion parallel to it. Energy estimates are shown to for the continuous and semi-discrete cases. Numerical experiments are performed showing convergence of the method, and examples is given to demonstrate the capabilities of the method.

We develop a hybrid spatial discretization for the wave equation in second order form, based on high-order accurate finite difference methods and discontinuous Galerkin methods. The hybridization combines computational efficiency of finite difference methods on Cartesian grids and geometrical flexibility of discontinuous Galerkin methods on unstructured meshes. The two spatial discretizations are coupled by a penalty technique at the interface such that the overall semidiscretization satisfies a discrete energy estimate to ensure stability. In addition, optimal convergence is obtained in the sense that when combining a fourth order finite difference method with a discontinuous Galerkin method using third order local polynomials, the overall convergence rate is fourth order. Furthermore, we use a novel approach to derive an error estimate for the semidiscretization by combining the energy method and the normal mode analysis for a corresponding one dimensional model problem. The stability and accuracy analysis are verified in numerical experiments.

In this paper, we present the stability analysis of the perfectly matched layer (PML) in two-space dimensional layered elastic media. Using normal mode analysis we prove that all interface wave modes present at a planar interface of bi-material elastic solids are dissipated by the PML. Our analysis builds upon the ideas presented in [SIAM Journal on Numerical Analysis 52 (2014) 2883-2904] and extends the stability results of boundary waves (such as Rayleigh waves) on a half-plane elastic solid to interface wave modes (such as Stoneley waves) transmitted into the PML at a planar interface separating two half-plane elastic solids. Numerical experiments in two-layer and multi-layer elastic solids corroborate the theoretical analysis, and generalise the results to complex elastic media. Numerical examples using the Marmousi model demonstrates the utility of the PML and our numerical method for seismological applications.

The reduction of computational costs in the numerical solution of nonstationary problems is achieved through splitting schemes. In this case, solving a set of less computationally complex problems provides the transition to a new level in time. The traditional construction approach of splitting schemes is based on an additive representation of the problem operator(s) and uses explicit-implicit approximations for individual terms. Recently (Y. Efendiev, P.N. Vabishchevich. Splitting methods for solution decomposition in nonstationary problems. \textit{Applied Mathematics and Computation}. \textbf{397}, 125785, 2021), a new class of methods of approximate solution of nonstationary problems has been introduced based on decomposition not of operators but of the solution itself. This new approach with subdomain solution selection is used in this paper to construct domain decomposition schemes. The boundary value problem for a second-order parabolic equation in a rectangle with a difference approximation in space is typical. Two and three-level schemes for decomposition of the domain with and without overlapping subdomains are investigated. Our numerical experiments complement the theoretical results., Comment: 22 pages, 10 figures

This is one of our series works on discrete energy analysis of the variable-step BDF schemes. In this part, we present stability and convergence analysis of the third-order BDF (BDF3) schemes with variable steps for linear diffusion equations, see e.g. [SIAM J. Numer. Anal., 58:2294-2314] and [Math. Comp., 90: 1207-1226] for our previous works on the BDF2 scheme. To this aim, we first build up a discrete gradient structure of the variable-step BDF3 formula under the condition that the adjacent step ratios are less than 1.4877, by which we can establish a discrete energy dissipation law. Mesh-robust stability and convergence analysis in the $L^2$ norm are then obtained. Here the mesh robustness means that the solution errors are well controlled by the maximum time-step size but independent of the adjacent time-step ratios. We also present numerical tests to support our theoretical results., Comment: 21 pages, 5 tables

The recently developed technique of DOC kernels has been a great success in the stability and convergence analysis for BDF2 scheme with variable time steps. However, such an analysis technique seems not directly applicable to problems with initial singularity. In the numerical simulations of solutions with initial singularity, variable time-steps schemes like the graded mesh are always adopted to achieve the optimal convergence, whose first adjacent time-step ratio may become pretty large so that the acquired restriction is not satisfied. In this paper, we revisit the variable time-step implicit-explicit two-step backward differentiation formula (IMEX BDF2) scheme presented in [W. Wang, Y. Chen and H. Fang, \emph{SIAM J. Numer. Anal.}, 57 (2019), pp. 1289-1317] to compute the partial integro-differential equations (PIDEs) with initial singularity. We obtain the sharp error estimate under a mild restriction condition of adjacent time-step ratios $r_{k}: =\tau_{k}/\tau_{k-1} \; (k\geq 3) < r_{\max} = 4.8645 $ and a much mild requirement on the first ratio, i.e., $r_2>0$. This leads to the validation of our analysis of the variable time-step IMEX BDF2 scheme when the initial singularity is dealt by a simple strategy, i.e., the graded mesh $t_k=T(k/N)^{\gamma}$. In this situation, the convergence of order $\mathcal{O}(N^{-\min\{2,\gamma \alpha\}})$ is achieved with $N$ and $\alpha$ respectively representing the total mesh points and indicating the regularity of the exact solution. This is, the optical convergence will be achieved by taking $\gamma_{\text{opt}}=2/\alpha$. Numerical examples are provided to demonstrate our theoretical analysis.

In this paper we consider a linearized variable-time-step two-step backward differentiation formula (BDF2) scheme for solving nonlinear parabolic equations. The scheme is constructed by using the variable time-step BDF2 for the linear term and a Newton linearized method for the nonlinear term in time combining with a Galerkin finite element method (FEM) in space. We prove the unconditionally optimal error estimate of the proposed scheme under mild restrictions on the ratio of adjacent time-steps, i.e. $0

We prove that the two-step backward differentiation formula (BDF2) method is stable on arbitrary time grids; while the variable-step BDF3 scheme is stable if almost all adjacent step ratios are less than 2.553. These results relax the severe mesh restrictions in the literature and provide a new understanding of variable-step BDF methods. Our main tools include the discrete orthogonal convolution kernels and an elliptic-type matrix norm., Comment: 20 pages, 2 figures, 3 tables

We develop a stable finite difference method for the elastic wave equation in bounded media, where the material properties can be discontinuous at curved interfaces. The governing equation is discretized in second order form by a fourth or sixth order accurate summation-by-parts operator. The mesh size is determined by the velocity structure of the material, resulting in nonconforming grid interfaces with hanging nodes. We use order-preserving interpolation and the ghost point technique to couple adjacent mesh blocks in an energy-conserving manner, which is supported by a fully discrete stability analysis. In our previous work for the wave equation, two pairs of order-preserving interpolation operators are needed when imposing the interface conditions weakly by a penalty technique. Here, we only use one pair in the ghost point method. In numerical experiments, we demonstrate that the convergence rate is optimal, and is the same as when a globally uniform mesh is used in a single domain. In addition, with a predictor-corrector time integration method, we obtain time stepping stability with stepsize almost the same as given by the usual Courant-Friedrichs-Lewy condition., Comment: 34 pages

We study necessary conditions for stability of a Numerov-type compact higher-order finite-difference scheme for the 1D homogeneous wave equation in the case of non-uniform spatial meshes. We first show that the uniform in time stability cannot be valid in any spatial norm provided that the complex eigenvalues appear in the associated mesh eigenvalue problem. Moreover, we prove that then the solution norm grows exponentially in time making the scheme strongly non-dissipative and therefore impractical. Numerical results confirm this conclusion. In addition, for some sequences of refining spatial meshes, an excessively strong condition between steps in time and space is necessary (even for the non-uniform in time stability) which is familiar for explicit schemes in the parabolic case., Comment: 6 pages

We consider finite difference approximations of the second derivative, exemplified in Poisson's equation, the heat equation and the wave equation. The finite difference operators satisfy a summation-by-parts property, which mimics the integration-by-parts. Since the operators approximate the second derivative, they are singular by construction. To impose boundary conditions, these operators are modified using Simultaneous Approximation Terms. This makes the modified matrices non-singular, for most choices of boundary conditions. Recently, inverses of such matrices were derived. However, when considering Neumann boundary conditions on both boundaries, the modified matrix is still singular. For such matrices, we have derived an explicit expression for the Moore-Penrose pseudoinverse, which can be used for solving elliptic problems and some time-dependent problems. The condition for this new pseudoinverse to be valid, is that the modified matrix does not have more than one zero eigenvalue. We have reconstructed the sixth order accurate narrow-stencil operator with a free parameter and show that more than one zero eigenvalue can occur. We have performed a detailed analysis on the free parameter to improve the properties of the second derivative operator. We complement the derivations by numerical experiments to demonstrate the improvements of the new second derivative operator.

In this work, we are concerned with the stability and convergence analysis of the second order BDF (BDF2) scheme with variable steps for the molecular beam epitaxial model without slope selection. We first show that the variable-step BDF2 scheme is convex and uniquely solvable under a weak time-step constraint. Then we show that it preserves an energy dissipation law if the adjacent time-step ratios $r_k:=\tau_k/\tau_{k-1}<3.561.$ Moreover, with a novel discrete orthogonal convolution kernels argument and some new discrete convolutional inequalities, the $L^2$ norm stability and rigorous error estimates are established, under the same step-ratios constraint that ensuring the energy stability., i.e., $0

We propose a structure-preserving finite difference scheme for the Cahn-Hilliard equation with a dynamic boundary condition using the discrete variational derivative method (DVDM). In this approach, it is important and essential how to discretize the energy which characterizes the equation. By modifying the conventional manner and using an appropriate summation-by-parts formula, we can use a standard central difference operator as an approximation of an outward normal derivative on the discrete boundary condition of the scheme. We show that our proposed scheme is second-order accurate in space, although the previous structure-preserving scheme by Fukao-Yoshikawa-Wada (Commun. Pure Appl. Anal. 16 (2017), 1915-1938) is first-order accurate in space. Also, we show the stability, the existence, and the uniqueness of the solution for the proposed scheme. Computation examples demonstrate the effectiveness of the proposed scheme. Especially through computation examples, we confirm that numerical solutions can be stably obtained by our proposed scheme., Comment: 32 pages, 18 figures

In this paper we study the stability of explicit finite difference discretizations of linear advection-diffusion equations (ADE) with arbitrary order of accuracy in the context of method of lines. The analysis first focuses on the stability of the system of ordinary differential equations (ODE) that is obtained by discretizing the ADE in space and then extends to fully discretized methods where explicit Runge-Kutta methods are used for integrating the ODE system. In particular, it is proved that all stable semi-discretization of the ADE gives rise to a conditionally stable fully discretized method if the time-integrator is at least first-order accurate, whereas high-order spatial discretization of the advection equation cannot yield a stable method if the temporal order is too low. In the second half of this paper, we extend the analysis to a partially dissipative wave system and obtain the stability results for both semi-discretized and fully-discretized methods. Finally, the major theoretical predictions are verified numerically., Comment: 27 pages, 11 figures

We develop a fourth order accurate finite difference method for the three dimensional elastic wave equation in isotropic media with the piecewise smooth material property. In our model, the material property can be discontinuous at curved interfaces. The governing equations are discretized in second order form on curvilinear meshes by using a fourth order finite difference operator satisfying a summation-by-parts property. The method is energy stable and high order accurate. The highlight is that mesh sizes can be chosen according to the velocity structure of the material so that computational efficiency is improved. At the mesh refinement interfaces with hanging nodes, physical interface conditions are imposed by using ghost points and interpolation. With a fourth order predictor-corrector time integrator, the fully discrete scheme is energy conserving. Numerical experiments are presented to verify the fourth order convergence rate and the energy conserving property.

The variable two-step backward differentiation formula (BDF2) is revisited via a new theoretical framework using the positive semi-definiteness of BDF2 convolution kernels and a class of orthogonal convolution kernels. We prove that, if the adjacent time-step ratios $r_k:=\tau_k/\tau_{k-1}\le(3+\sqrt{17})/2\approx3.561$, the adaptive BDF2 time-stepping scheme for linear reaction-diffusion equations is unconditionally stable and (maybe, first-order) convergent in the $L^2$ norm. The second-order temporal convergence can be recovered if almost all of time-step ratios $r_k\le 1+\sqrt{2}$ or some high-order starting scheme is used. Specially, for linear dissipative diffusion problems, the stable BDF2 method preserves both the energy dissipation law (in the $H^1$ seminorm) and the $L^2$ norm monotonicity at the discrete levels. An example is included to support our analysis., Comment: 20 pages

The threshold dynamics algorithm of Merriman, Bence, and Osher is only first order accurate in the two-phase setting. Its accuracy degrades further to half order in the multi-phase setting, a shortcoming it has in common with other related, more recent algorithms such as the equal surface tension version of the Voronoi implicit interface method. As a first, rigorous step in addressing this shortcoming, we present two different second order accurate versions of two-phase threshold dynamics. Unlike in previous efforts in this direction, we present careful consistency calculations for both of our algorithms. The first algorithm is consistent with its limit (motion by mean curvature) up to second order in any space dimension. The second achieves second order accuracy only in dimension two, but comes with a rigorous stability guarantee (unconditional energy stability) in any dimension -- a first for high order schemes of its type.

Efficient and accurate numerical simulation of 3D acoustic wave propagation in heterogeneous media plays an important role in the success of seismic full waveform inversion (FWI) problem. In this work, we employed the combined scheme and developed a new explicit compact high-order finite difference scheme to solve the 3D acoustic wave equation with spatially variable acoustic velocity. The boundary conditions for the second derivatives of spatial variables have been derived by using the equation itself and the boundary condition for $u$. Theoretical analysis shows that the new scheme has an accuracy order of $O(\tau^2) + O(h^4)$, where $\tau$ is the time step and $h$ is the grid size. Combined with Richardson extrapolation or Runge-Kutta method, the new method can be improved to 4th-order in time. Three numerical experiments are conducted to validate the efficiency and accuracy of the new scheme. The stability of the new scheme has been proved by an energy method, which shows that the new scheme is conditionally stable with a Courant - Friedrichs - Lewy (CFL) number which is slightly lower than that of the Pad\'{e} approximation based method., Comment: Submitted to Journal of Computational Science

Due to the intrinsically initial singularity of solution and the discrete convolution form in numerical Caputo derivatives, the traditional $H^1$-norm analysis (corresponding to the case for a classical diffusion equation) to the time approximations of a fractional subdiffusion problem always leads to suboptimal error estimates (a loss of time accuracy). To recover the theoretical accuracy in time, we propose an improved discrete Gr\"{o}nwall inequality and apply it to the well-known L1 formula and a fractional Crank-Nicolson scheme. With the help of a time-space error-splitting technique and the global consistency analysis, sharp $H^1$-norm error estimates of the two nonuniform approaches are established for a reaction-subdiffusion problems. Numerical experiments are included to confirm the sharpness of our analysis., Comment: 22 pages, 8 tables

We analyze two types of summation-by-parts finite difference operators for approximating the second derivative with variable coefficient. The first type uses ghost points, while the second type does not use any ghost points. A previously unexplored relation between the two types of summation-by-parts operators is investigated. By combining them we develop a new fourth order accurate finite difference discretization with hanging nodes on the mesh refinement interface. We take the model problem as the two-dimensional acoustic wave equation in second order form in terms of acoustic pressure, and prove energy stability for the proposed method. Compared to previous approaches using ghost points, the proposed method leads to a smaller system of linear equations that needs to be solved for the ghost point values. Another attractive feature of the proposed method is that the explicit time step does not need to be reduced relative to the corresponding periodic problem. Numerical experiments, both for smoothly varying and discontinuous material properties, demonstrate that the proposed method converges to fourth order accuracy. A detailed comparison of the accuracy and the time-step restriction with the simultaneous-approximation-term penalty method is also presented.

We consider wave propagation in a coupled fluid-solid region, separated by a static but possibly curved interface. The wave propagation is modeled by the acoustic wave equation in terms of a velocity potential in the fluid, and the elastic wave equation for the displacement in the solid. At the fluid solid interface, we impose suitable interface conditions to couple the two equations. We use a recently developed, energy based discontinuous Galerkin method to discretize the governing equations in space. Both energy conserving and upwind numerical fluxes are derived to impose the interface conditions. The highlights of the developed scheme include provable energy stability and high order accuracy. We present numerical experiments to illustrate the accuracy property and robustness of the developed scheme.

We study the Back and Forth Error Compensation and Correction (BFECC) method for linear hyperbolic PDE systems. The BFECC method has been applied to schemes for advection equations to improve their stability and order of accuracy. Similar results are established in this paper for schemes for linear hyperbolic PDE systems with constant coefficients. We apply the BFECC method to central difference scheme and Lax-Friedrichs scheme for the Maxwell's equations and obtain second order accurate schemes with larger CFL number than the classical Yee scheme. The method is further applied to schemes on non-orthogonal unstructured grids. The new BFECC schemes for the Maxwell's equations operate on a single non-staggered grid and are simple to implement on unstructured grids. Numerical examples are given to demonstrate the effectiveness of the new schemes.

In this paper, by using Strang's second-order splitting method, the numerical procedure for the three-dimensional (3D) space fractional Allen-Cahn equation can be divided into three steps. The first and third steps involve an ordinary differential equation, which can be solved analytically. The intermediate step involves a 3D linear fractional diffusion equation, which is solved by the Crank-Nicolson alternating directional implicit (ADI) method. The ADI technique can convert the multidimensional problem into a series of one-dimensional problems, which greatly reduces the computational cost. A fourth-order difference scheme is adopted for discretization of the space fractional derivatives. Finally, Richardson extrapolation is exploited to increase the temporal accuracy. The proposed method is shown to be unconditionally stable by Fourier analysis. Another contribution of this paper is to show that the numerical solutions satisfy the discrete maximum principle under reasonable time step constraint. For fabricated smooth solutions, numerical results show that the proposed method is unconditionally stable and fourth-order accurate in both time and space variables. In addition, the discrete maximum principle is also numerically verified., Comment: 24 pages, 7 figures, 10 tables

A new parametric class of semi-implicit numerical schemes for a level set advection equation on Cartesian grids is derived and analyzed. An accuracy and a stability study is provided for a linear advection equation with a variable velocity using partial Lax-Wendroff procedure and numerical von Neumann stability analysis. The obtained semi-implicit kappa-scheme is 2nd order accurate in space and time in any dimensional case when using a dimension by dimension extension of the one-dimensional scheme that is not the case for analogous fully explicit or fully implicit kappa-schemes. A further improvement is obtained by using so-called Corner Transport Upwind extension in two-dimensional case. The extended semi-implicit kappa-scheme with a specific (velocity dependent) value of kappa is 3rd order accurate in space and time for a constant advection velocity, and it is unconditional stable according to the numerical von Neumann stability analysis for the linear advection equation in general., Comment: arXiv admin note: substantial text overlap with arXiv:1611.04153 Comment from the authors - this is a corrected paper where typesetting error in formula (37) is removed

The modified Hunter--Saxton equation models the propagation of short capillary-gravity waves. As it involves a mixed derivative, its initial value problem on the periodic domain is much more complicated than the standard evolutionary equations. Although its local well-posedness has recently been proved, the behavior of its solution is yet to be investigated. In this paper, to develop a reliable numerical method for this problem, we derive a conservative finite difference scheme. Then, we rigorously prove not only its stability in the sense of the uniform norm but also its uniform convergence to sufficiently smooth exact solutions. Discrete conservation laws are used to overcome the difficulty due to the mixed derivative.

The stability bounds and error estimates for a compact higher order Numerov-Crank-Nicolson scheme on non-uniform space meshes for the 1D time-dependent Schr\"odinger equation have been recently derived. This analysis has been done in $L^2$ and $H^1$ mesh norms and used the non-standard "converse" condition $h_\omega\leq c_0\tau$, where $h_\omega$ is the mean space step, $\tau$ is the time step and $c_0>0$. Now we prove that such condition is necessary for some families of non-uniform meshes and any space norm. Also numerical results show unacceptably wrong behavior of numerical solutions (their dramatic mass non-conservation) when this condition is violated., Comment: 6 pages, 2 figures, 1 table

This paper presents an extension of a recently developed high order finite difference method for the wave equation on a grid with non-conforming interfaces. The stability proof of the existing methods relies on the interpolation operators being norm-contracting, which is satisfied by the second and fourth order operators, but not by the sixth order operator. We construct new penalty terms to impose interface conditions such that the stability proof does not require the norm-contracting condition. As a consequence, the sixth order accurate scheme is also provably stable. Numerical experiments demonstrate the improved stability and accuracy property.

This paper investigates oscillation-free stability conditions of numerical methods for linear parabolic partial differential equations with some example extrapolations to nonlinear equations. Not clearly understood, numerical oscillations can create infeasible results. Since oscillation-free behavior is not ensured by stability conditions, a more precise condition would be useful for accurate solutions. Using Von Neumann and spectral analyses, we find and explore oscillation-free conditions for several finite difference schemes. Further relationships between oscillatory behavior and eigenvalues is supported with numerical evidence and proof. Also, evidence suggests that the oscillation-free stability condition for a consistent linearization may be sufficient to provide oscillation-free stability of the nonlinear solution. These conditions are verified numerically for several example problems by visually comparing the analytical conditions to the behavior of the numerical solution for a wide range of mesh sizes., Comment: 12 pages, 6 figures, 1 table

An initial-boundary value problem for the 1D self-adjoint parabolic equation on the half-axis is solved. We study a broad family of two-level finite-difference schemes with two parameters related to averagings both in time and space. Stability in two norms is proved by the energy method. Also discrete transparent boundary conditions are rigorously derived for schemes by applying the method of reproducing functions. Results of numerical experiments are included as well., Comment: 28 pages, 4 figures

To solve numerically boundary value problems for parabolic equations with mixed derivatives, the construction of difference schemes with prescribed quality faces essential difficulties. In parabolic problems, some possibilities are associated with the transition to a new formulation of the problem, where the fluxes (derivatives with respect to a spatial direction) are treated as unknown quantities. In this case, the original problem is rewritten in the form of a boundary value problem for the system of equations in the fluxes. This work deals with studying schemes with weights for parabolic equations written in the flux coordinates. Unconditionally stable flux locally one-dimensional schemes of the first and second order of approximation in time are constructed for parabolic equations without mixed derivatives. A peculiarity of the system of equations written in flux variables for equations with mixed derivatives is that there do exist coupled terms with time derivatives.

We propose a fully discretised numerical scheme for the hyperelastic rod wave equation on the line. The convergence of the method is established. Moreover, the scheme can handle the blow-up of the derivative which naturally occurs for this equation. By using a time splitting integrator which preserves the invariants of the problem, we can also show that the scheme preserves the positivity of the energy density.

Several relaxation approximations to partial differential equations have been recently proposed. Examples include conservation laws, Hamilton-Jacobi equations, convection-diffusion problems, gas dynamics problems. The present paper focuses onto diffusive relaxation schemes for the numerical approximation of nonlinear parabolic equations. These schemes are based on suitable semilinear hyperbolic system with relaxation terms. High order methods are obtained by coupling ENO and WENO schemes for space discretization with IMEX schemes for time integration. Error estimates and convergence analysis are developed for semidiscrete schemes with numerical analysis for fully discrete relaxed schemes. Various numerical results in one and two dimension illustrate the high accuracy and good properties of the proposed numerical schemes. These schemes can be easily implemented for parallel computer and applied to more general system of nonlinear parabolic equations in two- and three-dimensional cases., Comment: 27 pages, 3 tables, 2 figures Submitted to: SINUM