In this paper, a parallelizable high-order split-exponential integrator is presented for solving the scalar wave equation, which is matrix-free with respect to the Laplacian operator. The spatial operator is split into a nilpotent component, which holds the Laplacian operator, and a scalar component, which is a 2 × 2 scalar matrix. Each component’s exponential can be computed exactly, allowing for favorable stability and performance when high-order optimized splitting strategies are used. This is paired with Gauss-Lobatto quadrature for integration, allowing for fewer evaluations of the exponential and simple parallelization. Both numerical experiments and analysis are provided to verify the efficacy of the method.
This work addresses a weakly coupled system of two singularly perturbed linear convection-diffusion equations subject to Robin boundary conditions. The solution of the system exhibits overlapping boundary layers. The problem is solved on the Vulanović-Bakhvalov mesh using the standard upwind scheme, and a rigorous analysis is presented to establish first-order convergence of the proposed numerical method. To the best of our knowledge, parameter-uniform convergence results for weakly coupled systems with Robin boundary conditions on the Vulanović-Bakhvalov mesh have not been previously reported. Numerical experiments are provided to validate the theoretical results, and a comparison with solutions obtained on the standard Shishkin mesh is included. To demonstrate the efficiency of the mesh for more general systems, we also consider a problem involving three coupled equations.
The numerical approximation of some Boussinesq systems in two spatial dimensions is here considered. The differential systems under study are proposed as asymptotic models for the propagation of waves along the interface of two layers of fluids with different densities and subjected to a Boussinesq physical regime in each layer. Well-posedness of the periodic initial-value problem (ivp) of the systems is first analyzed. Then, a discretization in space based on the spectral Fourier-Galerkin method is introduced and error estimates for the semidiscrete approximation are derived. Using an efficient time integrator, some numerical experiments to illustrate the performance of the discretization are presented.
Singular behavior near the initial time, induced by weakly singular kernels in Volterra integro-differential equations, often degrades the accuracy of classical spectral collocation methods. To mitigate this effect, considerable attention has been devoted to the development of piecewise collocation methods employing fractional-order polynomials as basis functions, which effectively capture weakly singular or irregular solution behavior. In this paper, we consider the double weakly singular Volterra integro-differential equations and apply a piecewise collocation method based on a fractional-order polynomial basis. In addition, we establish a rigorous regularity analysis of the solution to this problem. We employ a graded mesh to overcome the low convergence orders usually given in collocation methods based on uniform meshes. Based on the presented analysis, we derive error estimates for the proposed method and investigate the selection of the fractional-order and the mesh grading parameter to achieve an optimal convergence rate of O(M−N), where M denotes the number of subintervals and N denotes the number of collocation points in each subinterval. Numerical results demonstrate the validity and accuracy of the theoretical results and confirm the effectiveness of the numerical approach.
A compact finite difference hybrid scheme is proposed for the Savage–Hutter equations. The method combines compact average–point-value conversion, adaptive central–Rusanov fluxes based on a Weighted Essentially Non-Oscillatory Z (WENO-Z) weight-deviation sensor, third-order Strong Stability Preserving Runge–Kutta (SSP-RK3) time integration, and wet–dry treatment. The compact spatial discretization provides nearly fourth-order accuracy for smooth solutions, while the sensor-based hybrid flux preserves the central compact flux in smooth regions and introduces Rusanov dissipation only near non-smooth regions and wet–dry fronts. Numerical tests demonstrate stable front capturing and excellent global mass conservation in moving-front granular flows.
We propose a novel dual-inertial subgradient extragradient method with extrapolation from the past for solving variational inequalities in Hilbert spaces. The algorithm uniquely combines Nesterov-type momentum, single operator evaluation per iteration, half-space projections, and a fully adaptive step-size rule requiring no Lipschitz constant. Under pseudo-monotonicity and Lipschitz continuity we establish weak convergence, the first such result for a subgradient extragradient scheme with dual Nesterov inertia and past extrapolation, while under strong pseudo-monotonicity we prove linear convergence at an explicit rate 1−ηαL with α≥3+2. Properties of the adaptive step-size sequence are analyzed, and practical parameter guidelines are provided. Numerical experiments on finite- and infinite-dimensional problems as well as compressed sensing show that the proposed algorithm consistently outperforms state-of-the-art methods in iteration count and computational time while maintaining accuracy and robustness.
In this paper, we develop a novel bound-preserving flux limiting scheme within the discontinuous Galerkin framework for one-dimensional scalar nonlinear conservation laws and convection-diffusion equations. The proposed scheme integrates the global monolithic convex (GMC) limiter with the generalized local Lax–Friedrichs flux, which effectively eliminates maximum principle violations around shocks. Compared with the classical flux-corrected transport technique, the GMC limiter introduces fewer constraints and enjoys natural compatibility with spatial semi-discretizations. In addition, the generalized local Lax–Friedrichs flux contains two tunable parameters, which facilitate accurate shock capturing and high-resolution simulation of smooth solutions. The established framework is further extended to handle one-dimensional convection-diffusion problems. In addition, we couple the GMC strategy with several typical high-order Runge–Kutta methods and analyze the bound-preserving property of the corresponding fully discrete scheme. A series of numerical examples involving one- and two-dimensional problems are presented. Numerical results demonstrate that the developed scheme preserves physical bounds near shocks while maintaining high-order accuracy in smooth regions.
In this paper, we propose and analyze a second-order finite difference method for a one-dimensional variable-coefficient damped nonlinear Klein–Gordon equation subject to periodic boundary conditions. The heterogeneous spatial operator is discretized by a conservative centered difference approximation, while the nonlinear force is treated using a discrete-gradient approximation associated with the original potential. The resulting fully discrete scheme satisfies an exact discrete energy-dissipation law, which reduces to exact energy conservation in the absence of damping. Unique solvability of the nonlinear algebraic system is established under an explicit time-step restriction that is independent of the spatial mesh size. We further derive a defect-corrected conformal symplectic balance, thereby clarifying the structural effect of the discrete-gradient treatment of the nonlinearity. By combining the discrete energy estimate, a one-dimensional discrete Sobolev inequality, and consistency estimates, we prove second-order convergence in the coefficient-weighted discrete H1-norm for the solution and in the discrete L2-norm for its velocity. Numerical experiments confirm the predicted convergence rate, the exact discrete energy identity, the long-time dissipative behavior, and the applicability of the method to wave propagation in heterogeneous media.
We study the conforming discontinuous Galerkin finite element method for solving the displacement obstacle problem of Kirchhoff plates on bounded polygonal domains in two dimensions. Under the weak complementarity condition, the error estimate in a discrete H2 norm for the quadratic method is O(hα), where α∈(12,1] is determined by the geometry of the polygonal domain. By imposing additional constraints on the contact set to ensure enhanced solution regularity, we establish the optimal error estimate with α∈(1,32) for the cubic method. The theoretical results are verified by numerical experiments.
In this paper, we present a fast numerical method for solving the parabolic Monge-Ampère (PMA) equation derived in the work of Sulman et al. [J. Comput. Phys., 230 (2011), pp. 3302–3330] for computing adaptive moving mesh in higher dimensions. The method is based on the exploitation of the Anderson acceleration algorithm developed by H. F. Walker and P. NI [SIAM J. Numer. Anal., Vol. 49, No. 4, pp. 1715–1735] for the fixed-point iteration. Anderson acceleration method has been successfully employed in electronic structure computations, and it has been used for various other engineering applications. Here, a standard finite difference scheme is used for the temporal discretization of the PMA equation which we then write as a fixed point iterations. We employ the Anderson acceleration (AA) method to the fixed point iterations to expedited the process of computing the steady state solution of the PMA equation. The proposed method is applied for computing adaptive meshes in two and three spatial dimensions. Several numerical experiments are presented to demonstrate the performance of the proposed accelerated AAPMA adaptive moving mesh method.
The power method is a simple and effective iteration method for computing the dominant eigenvalue and its corresponding eigenvector of a given real matrix. To further improve its convergence property, we propose randomized power methods in which only a subset of columns of the matrix and the corresponding subset of elements of a certain vector are sampled randomly so that the matrix-vector multiplications involved are computed approximately. In this way, the power method is technically randomized and can be effectively executed to solve the eigenvalue problems with respect to extremely large matrices especially when they can not be stored entirely in a computing resource. We adopt two sample strategies: sampling with replacement and sampling without replacement. By minimizing the expected Euclidean norm of error at each iteration step we obtain the optimal sample probabilities that are utilized in the randomized power methods. In addition, we discuss the basic properties associated with the randomized power methods, and show their numerical advantages over the classical power method.
We present an original preconditioner to solve the discretized Poisson equation with the conjugate gradient method. The preconditioner is based on the approximation of the inverse Poisson operator. In particular, we construct the inverse operator corresponding to a layered model which is close to the original one. The action of the preconditioner at each iteration requires the system of equations to be solved corresponding to the layered media. It is done by applying the spectral decomposition of the matrix corresponding to the 1D problem (spectral decomposition along one spatial direction) and further direct solution of a series of 1D problems along the other direction. We provide the numerical and theoretical study of the convergence rate depending on the way the layered model is constructed to approximate the original. We consider four cases: maximal value over the layer, minimal value, arithmetic averaging, and harmonic averaging. For the first two cases, we prove analytically that the convergence rate of the preconditioned conjugate gradient method is independent of the size of the problem but depends on the coefficient contrast. For the other two cases, the weak dependence on the size of the problem is illustrated numerically.
In this paper, we propose a structure-preserving Galerkin method for the Gross-Pitaevskii equation with angular momentum rotation based on the scalar auxiliary variable (SAV) formulation, which consists of a Crank-Nicolson temporal discretization and the finite element spatial discretization. The proposed scheme is proved to conserve both mass and modified energy exactly at the discrete level. By employing Schaefer’s fixed point theorem together with an H−1 duality argument, we establish the well-posedness and optimal H1-norm error estimates for the numerical solution without any grid-ratio condition. Optimal-order L2-norm convergence is further recovered by means of a discrete Abel summation formula and a temporal regularity estimate that compensate for the order reduction introduced by the SAV reformulation. Numerical examples are provided to validate the conservation properties and convergence rates of the proposed scheme.
In this work, we investigate a Thermo-Hydro-Mechanical (THM) model describing compressible flow in deformable porous media. Such models play a significant role in various areas of geomechanics, with applications ranging from underground energy storage to oil and gas reservoir engineering. The simulations of THM models are essential for the analysis and design of safe and efficient energy storage cavities. We begin by presenting the mathematical formulation of the model, which consists of a system of parabolic partial differential equations governing the conservation of fluid mass, entropy, and skeleton momentum. First, we derive energy estimates for the compressible model. Next, we focus in particular on the development of numerical discretizations that preserve the energy estimates of the continuous system at the discrete level. To this end, we employ the backward Euler scheme for the time discretization together with the finite volume two-point flux approximation (TPFA) method for the spatial discretization. Numerical experiments are presented to validate the proposed approach and to illustrate its accuracy, efficiency, and robustness.
This paper proposes a direct discontinuous Galerkin scheme to numerically solve stochastic convection-diffusion problems. First, the coercivity property of the constructed bilinear form corresponding to our numerical method is established by introducing interface correction terms. Following this, we rigorously prove the L2-stability of the proposed method. Then, the optimal convergence rate can be proved by using a special global projection operator to eliminate or control difficult terms at cell interfaces. Finally, we implement a set of numerical tests to confirm the rationality of our theoretical results.
In this paper, we construct and analyze two efficient numerical methods for a special class of highly oscillatory integrals arising from solving the Helmholtz equation and oscillatory integral equations. Such integrals have an integrand that is the product of the Hankel function of the first kind and a complex exponential function, and contains algebraic singularities at two endpoints. By employing the Meijer-G function, we build both crucial closed-form expressions and asymptotic estimates for the useful generalized moment integral(1)(0) x(alpha) (1 - x)(beta) e(t omega xr) IIv(1)(omega x(p))dx, where IIv(1) is the Hankel function of the first kind of order v. Moreover, we rigorously derive a fundamental recursive relation governing the required modified moment in question. These key results establish the theoretical foundation for implementing and analyzing the two proposed methods. Notably, we conduct rigorous error analysis through theoretical verification. The validity of our theoretical framework and the efficacy of the proposed methods are subsequently confirmed via numerical experiments. Finally, we provide a comprehensive numerical comparison between the two approaches, concluding with a systematic summary of their respective merits and limitations.
