Stability Preserving Research Articles

Strong Stability Preserving IMEX Methods for Partitioned Systems of Differential Equations

We investigate strong stability preserving (SSP) implicit-explicit (IMEX) methods for partitioned systems of differential equations with stiff and nonstiff subsystems. Conditions for order p and stage order q=p are derived, and characterization of SSP IMEX methods is provided following the recent work by Spijker. Stability properties of these methods with respect to the decoupled linear system with a complex parameter, and a coupled linear system with real parameters are also investigated. Examples of methods up to the order p=4 and stage order q=p are provided. Numerical examples on six partitioned test systems confirm that the derived methods achieve the expected order of convergence for large range of stepsizes of integration, and they are also suitable for preserving the accuracy in the stiff limit or preserving the positivity of the numerical solution for large stepsizes.

Spectral discretizations analysis with time strong stability preserving properties for pseudo-parabolic models

In this work, we study the numerical approximation of the initial-boundary-value problem of nonlinear pseudo-parabolic equations with Dirichlet boundary conditions. We propose a discretization in space with spectral schemes based on Jacobi polynomials and in time with robust schemes attending to qualitative features such as stiffness and preservation of strong stability for a more correct simulation of non-regular data. Error estimates for the corresponding semidiscrete Galerkin and collocation schemes are derived. The performance of the fully discrete methods is analyzed in a computational study.

The meshless local Petrov–Galerkin method based on moving Taylor polynomial approximation to investigate unsteady diffusion–convection problems of anisotropic functionally graded materials related to incompressible flow

This paper concerns to a meshless local Petrov–Galerkin (MLPG) method for studying the unsteady diffusion–convection problems of anisotropic functionally graded materials. A new version of MLPG method based on the moving Taylor polynomial approximation is developed to discrete the spatial variable. Then, we obtain a system of ODEs which depends to the time variable. A strong stability preserving (SSP) Runge–Kutta idea is provided to solve the final ODEs with enough accuracy and stability. Also, the grading function which defines the variable elastic coefficient can be any types of continuous functions. The developed numerical formulation is applied for different examples of non-rectangular domains to check its accuracy.

High Order Semi-implicit Multistep Methods for Time-Dependent Partial Differential Equations

We consider the construction of semi-implicit linear multistep methods that can be applied to time-dependent PDEs where the separation of scales in additive form, typically used in implicit-explicit (IMEX) methods, is not possible. As shown in Boscarino et al. (J. Sci. Comput. 68: 975–1001, 2016) for Runge-Kutta methods, these semi-implicit techniques give a great flexibility, and allow, in many cases, the construction of simple linearly implicit schemes with no need of iterative solvers. In this work, we develop a general setting for the construction of high order semi-implicit linear multistep methods and analyze their stability properties for a prototype linear advection-diffusion equation and in the setting of strong stability preserving (SSP) methods. Our findings are demonstrated on several examples, including nonlinear reaction-diffusion and convection-diffusion problems.

Computation of Optimal Linear Strong Stability Preserving Methods Via Adaptive Spectral Transformations of Poisson–Charlier Measures

Strong stability preserving (SSP) coefficients govern the maximally allowable step-size at which positivity or contractivity preservation of integration methods for initial value problems is guaranteed. In this paper, we show that the task of computing optimal linear SSP coefficients of explicit one-step methods is, to a certain extent, equivalent to the problem of characterizing positive quadratures with integer nodes with respect to Poisson–Charlier measures. Using this equivalence, we provide sharp upper and lower bounds for the optimal linear SSP coefficients in terms of the zeros of generalized Laguerre orthogonal polynomials. This in particular provides us with a sharp upper bound for the optimal SSP coefficients of explicit Runge–Kutta methods. Also based on this equivalence, we propose a highly efficient and stable algorithm for computing these coefficients, and their associated optimal linear SSP methods, based on adaptive spectral transformations of Poisson–Charlier measures. The algorithm possesses the remarkable property that its complexity depends only on the order of the method and thus is independent of the number of stages. Our results are achieved by adapting and extending an ingenious technique by Bernstein (Acta Math 52:1–66, 1928) in his seminal work on absolutely monotonic functions. Moreover, the techniques introduced in this work can be adapted to solve the integer quadrature problem for any positive discrete multi-parametric measure supported on \({\mathbb {N}}\) under some mild conditions on the zeros of the associated orthogonal polynomials.

Krylov SSP Integrating Factor Runge–Kutta WENO Methods

Weighted essentially non-oscillatory (WENO) methods are especially efficient for numerically solving nonlinear hyperbolic equations. In order to achieve strong stability and large time-steps, strong stability preserving (SSP) integrating factor (IF) methods were designed in the literature, but the methods there were only for one-dimensional (1D) problems that have a stiff linear component and a non-stiff nonlinear component. In this paper, we extend WENO methods with large time-stepping SSP integrating factor Runge–Kutta time discretization to solve general nonlinear two-dimensional (2D) problems by a splitting method. How to evaluate the matrix exponential operator efficiently is a tremendous challenge when we apply IF temporal discretization for PDEs on high spatial dimensions. In this work, the matrix exponential computation is approximated through the Krylov subspace projection method. Numerical examples are shown to demonstrate the accuracy and large time-step size of the present method.

Maximum bound principle preserving integrating factor Runge–Kutta methods for semilinear parabolic equations

A large class of semilinear parabolic equations satisfy the maximum bound principle (MBP) in the sense that the time-dependent solution preserves for any time a uniform pointwise bound imposed by its initial and boundary conditions. The MBP plays a crucial role in understanding the physical meaning and the wellposedness of the mathematical model. Investigation on numerical algorithms with preservation of the MBP has attracted increasingly attentions in recent years, especially for the temporal discretizations, since the violation of MBP may lead to nonphysical solutions or even blow-ups of the algorithms. In this paper, we study high-order MBP-preserving time integration schemes by means of the integrating factor Runge–Kutta (IFRK) method. Beginning with the space-discrete system of semilinear parabolic equations, we present the IFRK method in general form and derive the sufficient conditions for the method to preserve the MBP. In particular, we show that the classic four-stage, fourth-order IFRK scheme is MBP preserving for some typical semilinear systems although not strong stability preserving, which can be instantly applied to the Allen–Cahn type of equations. To our best knowledge, this is the first time to present a fourth-order linear numerical method preserving the MBP. In addition, convergence of these numerical schemes is proved theoretically and verified numerically, as well as their efficiency by simulations of 2D and 3D long-time evolutional behaviors. Numerical experiments are also carried out for a model which is not a typical gradient flow as the Allen–Cahn type of equations.

A High-Order Maximum-Principle-Satisfying Discontinuous Galerkin Method for the Level Set Problem

Level set (LS) method is a widely used interface capturing method. In the simulations of incompressible two-phase flows, in order to avoid discontinuities at interfaces, the LS function is usually taken as a smeared-out Heaviside function bounded on [0, 1] and advected by a given velocity field $$\mathbf {u}$$ obtained from the solution of the incompressible Navier-Stokes equations. In the incompressible limit $$\nabla \cdot \mathbf {u}=0$$ , the advection equation for the LS function can be written and discretized in conservative form. However, due to numerical errors, the resulting velocity field is in general not divergence free which leads to the solution of the advection equation in conservative form does not satisfy the maximum principle. To overcome this issue, in this work, we develop a high-order discontinuous Galerkin (DG) method to directly solve the advection equation for the LS function in non-conservative form. Moreover, we prove that by applying a linear scaling limiter, the proposed method together with a strong stability preserving (SSP) time discretization scheme can satisfy the strict maximum principle under a suitable CFL condition. Numerical simulations of several well-known benchmark problems, including the application to incompressible two-phase flows, are presented to demonstrate the high-order accuracy and maximum-principle-satisfying property of the proposed method.

A comparative analysis of explicit, IMEX and implicit strong stability preserving Runge-Kutta schemes

Selecting between explicit and implicit time marching schemes, as well as non-iterative explicit-implicit ones, is well understood on a linear stability basis. Implicit schemes are favored when either searching for steady-states or performing time accurate simulations of stiff equations due to their larger time steps. When nonlinear stability properties must be preserved, however, solution monotonicity is guaranteed by a time step restriction that puts into question the usefulness of implicit schemes. Five test cases are used in this paper for a comparative efficiency analysis of explicit, implicit and implicit-explicit Strong Stability Preserving (SSP) Runge-Kutta (RK) schemes. They indicate that explicit SSP-RK schemes are, in general, the most efficient ones. It is important to emphasize that great care must be exercised when generalizing these conclusions, since such an analysis is inherently software, hardware, numerics and model dependent. Nevertheless, this study provides the first systematic evidence in the literature confirming these trends, to the best of our knowledge. Finally, implicit-explicit SSP-RK schemes do perform well, even though they are not optimal. Hence, there is a need for the development of their optimal versions before any conclusions can be drawn about their efficiency.

High order conservative Lagrangian schemes for one-dimensional radiation hydrodynamics equations in the equilibrium-diffusion limit

Radiation hydrodynamics (RH) describes the interaction between matter and radiation which affects the thermodynamic states and the dynamic flow characteristics of the matter-radiation system. Its application areas are mainly in high-temperature hydrodynamics, including gaseous stars in astrophysics, combustion phenomena, reentry vehicles fusion physics and inertial confinement fusion (ICF). Solving the radiation hydrodynamics equations (RHE), even in the equilibrium-diffusion limit, is a difficult task. In this paper, we will discuss the methodology to construct fully explicit and implicit-explicit (IMEX) high order Lagrangian schemes solving one dimensional RHE in the equilibrium-diffusion limit respectively, which can be used to simulate multi-material problems with the coupling of radiation and hydrodynamics. The schemes are based on the HLLC numerical flux, the essentially non-oscillatory (ENO) reconstruction for the advection term, ENO reconstruction or high order central reconstruction and interpolation for the radiation diffusion term, the Newton iteration method (for the IMEX scheme), and the strong stability preserving (SSP) high order time discretizations. The schemes can maintain conservation and uniformly high order accuracy both in space and time. The issue of positivity-preserving for the high order explicit Lagrangian scheme is also discussed. Various numerical tests for the high order Lagrangian schemes are provided to demonstrate the desired properties of the schemes such as high order accuracy, non-oscillation, and positivity-preserving.

Transformed implicit-explicit second derivative diagonally implicit multistage integration methods with strong stability preserving explicit part

In this paper, we discuss the construction of a class of implicit-explicit (IMEX) methods for systems of ordinary differential equations which their right hand side can be split into two parts; nonstiff or mildly stiff part and stiff part. The proposed methods treat the non-stiff part by an explicit second derivative diagonally implicit multistage integration method (SDIMSIM) and the stiff part by an implicit diagonally implicit multistage integration method (DIMSIM). The explicit part of these methods has strong stability preserving (SSP) property and the implicit part is A- and L-stable. We will construct methods with p=q=r=s and p=q+1=r=s up to order four with large SSP coefficients with respect to the large region of absolute stability, assuming that the implicit part of the method has Runge–Kutta stability (RKS) property together with A- and L-stability. These methods are tested on the linear advection-diffusion, advection-reaction and nonlinear shallow water equations, and the numerical results are presented conforming the efficiency and order of constructed methods.

Stability Analysis of an Alcoholism Model with Public Health Education and NSFD Scheme

In this paper, an alcoholism model of SEAR type with different susceptibilities due to public health education is investigated, with the form of continuous differential equations as well as discrete differential equations by applying the Mickens nonstandard finite difference (NSFD) scheme to the continuous equations. Threshold dynamics of the continuous model are performed by constructing Lyapunov functions. The analysis of a discrete model indicates that the alcohol-free equilibrium is globally asymptotically stable if the basic reproductive number R0<1, and conversely, the alcohol-present equilibrium is globally asymptotically stable if R0>1, revealing the consistency and efficiency of the discrete model to preserve the dynamical properties of the corresponding continuous model. In addition, stability preserving and the impact of the parameters related with public health education are conducted by numerical simulations.

A positivity-preserving Lagrangian Discontinuous Galerkin method for ideal magnetohydrodynamics equations in one-dimension

In this paper, we propose a conservative Lagrangian Discontinuous Galerkin (DG) scheme for solving the ideal compressible magnetohydrodynamics (MHD) equations in one-dimensional. This scheme can preserve positivity of physically positive variables such as density and thermal pressure. We first develop a Lagrangian HLLD approximate Riemann solver which can keep positivity-preserving property under some appropriate signal speeds. With this solver a first order positivity-preserving Lagrangian DG scheme can be constructed. Then we design a high order positivity-preserving and conservative Lagrangian DG scheme by using the strong stability preserving (SSP) high order time discretizations and the positivity-preserving scaling limiter. We adopt TVB minmod limiter with the local characteristic fields for the ideal MHD system to control spurious oscillations around the shock wave. Some numerical examples are presented to demonstrate the accuracy and positivity-preserving property of our scheme.

Strong stability preserving implicit–explicit transformed general linear methods

We consider the class of implicit–explicit (IMEX) general linear methods (GLMs) to construct methods where the explicit part has strong stability preserving (SSP) property, while the implicit part of the method has inherent Runge–Kutta stability (IRKS) property, and it is A-, or L-stable. We will also investigate the absolute stability of these methods when the implicit and explicit parts interact with each other. In particular, we will monitor the size of the region of absolute stability of the IMEX scheme, assuming that the implicit part of the method is A(α)-stable for α∈[0,π∕2]. Finally we furnish examples of SSP IMEX GLMs up to the order p=4 and stage order q=p with optimal SSP coefficients.

Matrix-free subcell residual distribution for Bernstein finite element discretizations of linear advection equations

In this work, we introduce a new residual distribution (RD) framework for the design of bound-preserving high-resolution finite element schemes. The continuous and discontinuous Galerkin discretizations of the linear advection equation are modified to construct local extremum diminishing (LED) approximations. To that end, we perform mass lumping and redistribute the element residuals in a manner which guarantees the LED property. The hierarchical correction procedure for high-order Bernstein finite element discretizations involves localization to subcells and definition of bound-preserving weights for subcell contributions. Using strong stability preserving (SSP) Runge–Kutta methods for time integration, we prove the validity of discrete maximum principles under CFL-like time step restrictions. The low-order version of our method has roughly the same accuracy as the one derived from a piecewise (multi)-linear approximation on a submesh with the same nodal points. In high-order extensions, we use an element-based flux-corrected transport (FCT) algorithm which can be interpreted as a nonlinear RD scheme. The proposed LED corrections are tailor-made for matrix-free implementations which avoid the rapidly growing cost of matrix assembly for high-order Bernstein elements. The results for 1D, 2D, and 3D test problems compare favorably to those obtained with the best matrix-based approaches.

Efficient SSP low-storage Runge–Kutta methods

In this paper we study the efficiency of Strong Stability Preserving (SSP) Runge–Kutta methods that can be implemented with a low number of registers using their Shu–Osher representation. SSP methods have been studied in the literature and stepsize restrictions that ensure numerical monotonicity have been found. However, for some problems, the observed stepsize restrictions are larger than the theoretical ones. Aiming at obtaining additional properties of the schemes that may explain their efficiency, in this paper we study the influence of the local error term in the observed stepsize restrictions. For this purpose, we consider the family of 5-stage third order SSP explicit Runge–Kutta methods, namely SSP(5,3), and the Buckley–Leverett equation. We deal with optimal SSP(5,3) schemes whose implementation requires at least 3 memory registers, and non-optimal 2-register SSP(5,3) schemes. The numerical experiments done show that small error constants improve the efficiency of the method in the sense that larger observed SSP coefficients are obtained.

Finite volume WENO schemes for nonlinear parabolic problems with degenerate diffusion on non-uniform meshes

We consider numerical approximation of the degenerate advection-diffusion equation, which is formally parabolic but may exhibit hyperbolic behavior. We develop both explicit and implicit finite volume weighted essentially non-oscillatory (WENO) schemes in multiple space dimensions on non-uniform computational meshes. The diffusion degeneracy is reformulated through the use of the Kirchhoff transformation. Space is discretized using WENO reconstructions with adaptive order (WENO-AO), which have several advantages, including the avoidance of negative linear weights and the ability to handle irregular computational meshes. A special two-stage WENO reconstruction procedure is developed to handle degenerate diffusion. Element averages of the solution are first reconstructed to give point values of the solution, and these point values are in turn used to reconstruct the Kirchhoff transform variable of the diffusive flux. Time is discretized using the method of lines and a Runge-Kutta time integrator. We use Strong Stability Preserving (SSP) Runge-Kutta methods for the explicit schemes, which have a severe parabolically scaled time step restriction to maintain stability. We also develop implicit Runge-Kutta methods. SSP methods are only conditionally stable, so we discuss the use of L-stable Runge-Kutta methods. We present in detail schemes that are third order in both space and time in one and two space dimensions using non-uniform meshes of intervals or quadrilaterals. Efficient implementation is described for computational meshes that are logically rectangular. Through a von Neumann (or Fourier mode) stability analysis, we show that smooth solutions to the linear problem are unconditionally L-stable on uniform computational meshes when using an implicit Radau IIA Runge-Kutta method. Computational results show the ability of the schemes to accurately approximate challenging test problems.

Strong Stability Preserving Integrating Factor Two-Step Runge–Kutta Methods

Problems with components that feature significantly different time scales, where the stiff time-step restriction comes from a linear component, implicit-explicit (IMEX) methods alleviate this restriction if the concern is linear stability. However, when nonlinear non-inner-product stability properties are of interest, such as in the evolution of hyperbolic partial differential equations with shocks or sharp gradients, linear inner-product stability is no longer sufficient for convergence, and so strong stability preserving (SSP) methods are often needed. Where the SSP property is needed, IMEX SSP Runge–Kutta (SSP-IMEX) methods have very restrictive time-steps. An alternative to SSP-IMEX schemes is to adopt an integrating factor approach to handle the linear component exactly and step the transformed problem forward using some time-evolution method. The strong stability properties of integrating factor Runge–Kutta methods were established in Isherwood et al. (SIAM J Numer Anal 56(6):3276–3307, 2018), where it was shown that it is possible to define explicit integrating factor Runge–Kutta methods that preserve strong stability properties satisfied by each of the two components when coupled with forward Euler time-stepping. It was proved that the solution will be SSP if the transformed problem is stepped forward with an explicit SSP Runge–Kutta method that has non-decreasing abscissas. However, explicit SSP Runge–Kutta methods have an order barrier of $$p=4$$, and sometimes higher order is desired. In this work we consider explicit SSP two-step Runge–Kutta integrating factor methods to raise the order. We show that strong stability is ensured if the two-step Runge–Kutta method used to evolve the transformed problem is SSP and has non-decreasing abscissas. We find such methods up to eighth order and present their SSP coefficients. Adding a step allows us to break the fourth order barrier on explicit SSP Runge–Kutta methods; furthermore, our explicit SSP two-step Runge–Kutta methods with non-decreasing abscissas typically have larger SSP coefficients than the corresponding one-step methods. A selection of our methods are tested for convergence and demonstrate the design order. We also show, for selected methods, that the SSP time-step predicted by the theory is a lower bound of the allowable time-step for linear and nonlinear problems that satisfy the total variation diminishing (TVD) condition. We compare some of the non-decreasing abscissa SSP two-step Runge–Kutta methods to previously found methods that do not satisfy this criterion on linear and nonlinear TVD test cases to show that this non-decreasing abscissa condition is indeed necessary in practice as well as theory. We also compare these results to our SSP integrating factor Runge–Kutta methods designed in Isherwood et al. (2018).

Strong Stability Preserving Second Derivative General Linear Methods

In this paper, we find sufficient strong stability preserving (SSP) conditions for second derivative general linear methods (SGLMs). Then we construct some optimal SSP SGLMs of order p up to eight and stage order $$1\le q\le p$$ with two external stages and $$2\le s\le 10$$ internal stages, which have larger effective Courant–Friedrichs–Lewy coefficients than the other existing class of the methods such as two derivative Runge–Kutta methods investigated by Christlieb et al., the class of two-step Runge–Kutta methods introduced by Ketcheson et al., the class of multistep multistage methods studied by Constantinescu and Sandu, and general linear methods investigated by Izzo and Jackiewicz. Some numerical experiments for scalar and systems of equations are provided which demonstrate that the constructed SSP SGLMs achieve the predicted order of convergence and preserve monotonicity.

Physical-constraints-preserving Lagrangian finite volume schemes for one- and two-dimensional special relativistic hydrodynamics

This paper studies the physical-constraints-preserving (PCP) Lagrangian finite volume schemes for one- and two-dimensional special relativistic hydrodynamic (RHD) equations. First, the PCP property (i.e. preserving the positivity of the rest-mass density and the pressure and the bound of the velocity) is proved for the first-order accurate Lagrangian scheme with the HLLC Riemann solver and forward Euler time discretization. The key is that the intermediate states in the HLLC Riemann solver are shown to be admissible or PCP when the HLLC wave speeds are estimated suitably. Then, the higher-order accurate schemes are proposed by using the high-order accurate strong stability preserving (SSP) time discretizations and the scaling PCP limiter as well as the WENO reconstruction. Finally, several one- and two-dimensional numerical experiments are conducted to demonstrate the accuracy and the effectiveness of the PCP Lagrangian schemes in solving the special RHD problems involving strong discontinuities, or large Lorentz factor, or low rest-mass density or low pressure, etc.