EFFICIENT INEQUALITY-PRESERVING INTEGRATORS FOR DIFFERENTIAL EQUATIONS SATISFYING FORWARD EULER CONDITIONS

. Developing explicit, high-order accurate, and stable algorithms for nonlinear differential equations remains an exceedingly difficult task. In this work, a systematic approach is proposed to develop high-order, large time-stepping schemes that can preserve inequality structures shared by a class of differential equations satisfying forward Euler conditions. Strong-stability-preserving (SSP) methods are popular and effective for solving equations of this type. However, few methods can deal with the situation when the time-step size is larger than that allowed by SSP methods. By adopting time-step-dependent stabilization and taking advantage of integrating factor methods in the Shu–Osher form, we propose enforcing the inequality structure preservation by approximating the exponential function using a novel recurrent approximation without harming the convergence. We define sufficient conditions for the obtained parametric Runge–Kutta (pRK) schemes to preserve inequality structures for any time-step size, namely, the underlying Shu–Osher coefficients are non-negative. To remove the requirement of a large stabilization term caused by stiff linear operators, we further develop inequality-preserving parametric integrating factor Runge–Kutta (pIFRK) schemes by incorporating the pRK with an integrating factor related to the stiff term, and enforcing the non-decreasing of abscissas. The only free parameter can be determined a priori based on the SSP coefficient, the time-step size, and the forward Euler condition. We demonstrate that the parametric methods developed here offer an effective and unified approach to study problems that satisfy forward Euler conditions, and cover a wide range of well-known models. Finally, numerical experiments reflect the high-order accuracy, efficiency, and inequality-preserving properties of the proposed schemes.

(1.2) Definition 1.2 (Positivity preservation). A method is positive if whenever 0 ≥ 0, it guarantees that ≥ 0, ∀ > 0 under the assumption that ∃ FE > 0, such that + ( , ) ≥ 0, ∀0 < ≤ FE , ∀ ∈ [0, ], ∀ ≥ 0, (1.3) where the inequalities are component-wise. Typical problems which highly depend on above inequality-preserving integrators include the diminishing of total variation in fluid dynamics [5], the positivity of density and pressure [6], the range boundedness of solutions in hyperbolic equations [7,8], and the maximum principle of order parameter in phase field models [4]. As the inequalities in Definitions 1.1-1.4 are characterized by forward Euler conditions (1.2)-(1.6), they can be preserved using the traditional SSP methods in most situations, as presented in the research of Hundsdorfer et al. [9], Ferracina and Spijker [10,11], Higueras [12,13], Gottlieb et al. [1]. Consider the preservation of strong stability in Definition 1.1. In the last three decades, significant efforts have been made to develop methods that reach the highest-order and admit the largest possible SSP coefficient for prescribed stage or step numbers, such that a method can preserve strong stability under the condition ≤ FE . Nevertheless, there are restrictions on the attainable value of for a given method, and SSP Runge-Kutta methods with non-negative coefficients have order barrier ≤ 4 for explicit schemes, ≤ 6 for implicit schemes [1,2,14]. Although, SSP multi-step schemes with non-negative coefficients have no order barrier [15,16], they require many steps to reach a given high-order. For studies on conditionally SSP methods, see [1,5,14,15,[17][18][19], and the references therein.
To remove restrictions on the time-step size, various studies have been focused on the development of unconditionally inequality-preserving schemes. Although the implicit Euler method can preserve strong stability without a time-step restriction (p. 10 in [9]), it has been proven that high-order implicit Runge-Kutta (RK) or linear multi-step (MS) methods cannot achieve this goal [9,20]. Therefore, many ideas have been proposed beyond the traditional framework. The diagonally split Runge-Kutta methods (DSRK), which were originally introduced to unconditionally preserve the contractivity [21,22], were adopted by Macdonald et al. [23] to preserve the strong stability. Unfortunately, it was shown that although second-and third-order unconditionally contractive DSRK methods preserved the strong stability for all time-step sizes, they suffered from order reduction at large step sizes. By incorporating negative coefficients and downwind-biased spatial discretization, Ketcheson [24] developed a class of second-order RK methods with an arbitrarily large SSP coefficient. Bonaventura and Della Rocca [3] proposed two unconditionally SSP extensions of the TR-BDF2 (trapezoidal rule-backward differentiation formula 2) method based on hybridization with the unconditionally SSP implicit Euler method, which can be activated using a sensor detecting violations of relevant functional bounds. Recently, by enforcing the backward derivative condition, Gottlieb et al. [25] proposed up to the fourth-order unconditionally SSP implicit two-derivative RK methods. Considering differential equations that have a special graph Laplacian structure, Blanes et al. [26] developed second-order methods that preserve positivity unconditionally and a third-order method that preserves positivity under very mild conditions. Other significant work in preserving the inequality structures was performed using the (modified) Patankar RK schemes [27][28][29].
To avoid solving implicit systems, explicit large time-stepping SSP schemes are desirable for numerically solving ODEs, and it is not surprising that such methods are difficult to construct. In this spirit, consider adding an ( 2 ) term ( − +1 ) to the forward Euler formula +1 = + ( , ). Then we construct a first-order implicit-explicit (IMEX) Euler scheme (1.7) Assuming that ≥ max {︁ We can see that if ≤ FE , with the diminishing of , the IMEX Euler scheme (1.7) becomes the forward Euler scheme, thus retaining the accuracy of the original scheme. When > FE , the introduction of the stabilization term ( − +1 ) stabilizes the scheme and preserves the strong stability by rescaling the time-step size from to 1+ ≤ FE . Although the construction of this first-order SSP method for any time-step size by adding an ( 2 ) term with a time-step-dependent (TSD) parameter is direct, severe time delay exists because of the time step rescaling, and the extension of this idea to higher order schemes is not straightforward.
In addition to the maximum-principle preservation for (1.1) in Definition 1.3, the preservation of the maximum principle for stiff ODE systems in the form 8) which are obtained by space discretization of Allen-Cahn-type (AC-type) equations, has attracted much attention in recent years. Du et al. [4] investigated sufficient conditions for the semilinear problem (1.8) to have a maximum principle by developing an abstract framework.
The nonlinear operator acts as a vector-valued function induced by a given continuously differentiable function 0 : R → R, that is ( ) = [ 0 ( 0 ), 0 ( 1 ), . . . , 0 ( −1 )] ∈ R , and there exists a constant > 0 such that Because of the change of sign on both sides of zero, a forward Euler condition on ( ) [30] holds, Unlike studies on the inequality-preserving methods in Definitions 1.1-1.4, the maximum-principle preservation of (1.8) has been independently investigated from another perspective. By taking advantage of the forward Euler condition (1.5), Tang and Yang [31] proposed the first maximum-principle-preserving (MPP) scheme by using the IMEX Euler formula with the TSD stabilization technique. Later, by adopting the central finite difference for the diffusion term and upwind discretization for the advection term, Shen et al. [32] investigated an unconditionally MPP semi-implicit scheme with the TSD stabilization for the generalized AC equation. Recently, by using a fourth-order spatial finite-difference discretization, Shen and Zhang [33] proved that the stabilized IMEX scheme was MPP under a suitable mesh-size constraint. However, traditional high-order temporal integrators fail to unconditionally preserve the maximum principle because of a lack of strategies for dealing with the stiff term and the stabilization term.
When considering the preservation of strong stability of systems in the form (1.8), Isherwood et al. [34] noted that a forward Euler condition on the linear term with respect to ‖ · ‖, that is, implies a contraction property ‖e ‖ ≤ ‖ ‖, ∀ > 0. Based on this finding, Isherwood et al. developed SSP integrating factor Runge-Kutta (IFRK) [34,35] and integrating factor two-step Runge-Kutta (IFTSRK) [36] schemes that only require the time-step size to satisfy a forward Euler condition on the nonlinear term. The IFRK approach has also been investigated for the maximum-principle preservation of (1.8) in [30,37]. To remove the time step restriction in the preservation of the maximum principle, a circle condition has been frequently incorporated. Let = 1 ≥ 1 FE , then the forward Euler condition (1.9) is equivalent to This so-called circle condition (1.10) [1,2,20] means that ( ) is bounded by a 'circle' measured by ‖ · ‖ ∞ that is centered at − ∈ R with radius . By utilizing exponential time difference (ETD) methods and the circle condition, Du et al. [38] developed two unconditionally MPP schemes for the non-local AC equation, that is, the stabilized first-order ETD1 scheme and the second-order ETD2 scheme. The ETD schemes have become popular in recent years. Du et al. [4] further proved that stabilized ETD schemes can preserve the maximum-principle unconditionally for a class of semilinear parabolic equations, and Ju et al. [39] showed that the ETD2 scheme incorporated with the generalized SAV and stabilization techniques unconditionally dissipates the energy of AC-type gradient flows. The success of stabilized ETD schemes lies in the fact that the contraction semi-group property of can be directly incorporated with the stabilization parameter in the temporal integrators, i.e., ‖e ( − ) ‖ ∞ ≤ e − ‖ ‖ ∞ . This stabilization technique has been introduced to the IFRK framework, in which a class of up to third-order unconditionally MPP stabilized IFRK schemes was proposed by Li et al. [40]. Zhang et al. proposed to approximate exponential functions in the stabilized IFRK, IFTSRK (integrating factor two-step Runge-Kutta), and IFMS (integrating factor multi-step) schemes using different Taylor series [41,42], recurrent approximations or combinations of exponential and linear functions [43,44]. The resulting up to the  fourth-order parametric IFRK, eighth-order parametric IFTSRK, and arbitrarily high-order parametric IFMS  schemes not only preserve the maximum-principle unconditionally but also conserve the mass for conservative  AC-type equations. Moreover, the stabilization technique was introduced to scalar hyperbolic equations with stiff source terms by Huang and Shu [8]. Starting from traditional th-order stabilized integrating factor Runge-Kutta formulas in the Shu-Osher form, they proposed replacing exponential functions with th-or ( +1)th-degree Taylor series to weaken the exponential effects produced by the stabilizing exponential term without destroying the convergence. Thus, the obtained high-order modified exponential SSP RK schemes only require the time-step size to satisfy the forward Euler condition on the nonlinear term, that is, ≤˜F E . In [45], Du et al. developed Taylor series approximation for stabilized integrating factor multi-step schemes. They constructed up to the third-order conservative, and bound-preserving schemes for stiff multispecies detonation. Later, Du et al. [46,47] proposed third-order modified exponential RK schemes in the Shu-Osher form by using conservative approximations comprising combinations of linear and exponential functions. To reduce the large truncation errors brought by the stabilization, Yang et al. [48] applied the standard RK scheme to the nonstiff equation and the modified exponential RK scheme to the stiff equation. If equations with different stiffness were encountered, they also considered applying the exponential RK method with a different stabilization parameter. This greatly improved the numerical accuracy. Still, the stabilization did not appear in the forward Euler condition of the nonlinear term, and the time-step size restriction existed.
In addition to the above studies on MPP and bound-preserving schemes using the stabilization technique, the idea of adding and subtracting a linear term has also been developed to improve the stability of numerical integrators, as demonstrated by Shen and Yang [49], as well as in earlier papers [50][51][52][53][54][55]. Examples include the introduction of ( − ) to AC-type equations [49], (∆ − ∆ ) to the mean-curvature or Cahn-Hilliard equations [51,54,56], to a surface diffusion equation [52] and a fourth-order image inpainting model [57], and ( − ) wtih = ( , ) ( ) to general nonlinear problem [58]. Such stabilization techniques were also referred to as Douglas splitting by Hundsdorfer [59], and the explicit-implicit-null (EIN) method by Duchemin and Eggers [60]. The first-and second-order EIN integrators were successfully incorporated with local discontinuous Galerkin discretizations to develop unconditionally stable schemes for nonlinear diffusion problems by Wang et al. [61]. To the best of our knowledge, there is still no stabilization method higher than first-order with a TSD stabilizer.
Since explicit, high-order accurate and stable schemes for both (1.1) and (1.8) are desirable, and a larger stabilization term introduces more truncation errors and time delay [48,62,63], the objective of this paper is to develop high-order inequality-preserving methods with TSD stabilization for both (1.1) and (1.8). The purpose is to increase the order (up to the fourth) and, simultaneously reduce additional truncation errors introduced by stabilization. To the best of our knowledge, there has been few study on explicit, high-order accurate and stable schemes to preserve inequality structures for (1.1), nor on direct incorporation of TSD stabilization in high-order schemes. Compared with the existing literature, the new contributions of this work include: -Recurrent approximations for exponential functions are proposed to develop up to the fourth-order, large time-stepping and inequality-preserving parametric RK schemes with TSD stabilization. -To tackle problems with both stiff and nonlinear terms, a family of explicit parametric integrating factor RK schemes that need small stabilization parameters are constructed to preserve the inequality structures. -The inequality-preservation of the pRK schemes has relatively mild requirements on the stabilization parameter and underlying schemes, that is, the stabilizing parameter ≥ max{ 1 FE − , 0} and all coefficients are non-negative. The parametric integrating factor RK schemes further require that all abscissas are nondecreasing and ≥ max{ 1 FE − , 0}. The rest of this paper is organized as follows. In Section 2, we provide sufficient conditions for (1.1) to admit inequality structures in Definitions 1.1-1.4. In Section 3, we start with a fixed-point-preserving improvement over the traditional integrating factor method, and illustrate main idea to construct up to the fourth-order inequality-preserving single-step schemes with TSD stabilization. The linear stability analysis and convergence estimate are carried out in Section 4. In Section 5, various experiments are performed to demonstrate the effectiveness and advantages of the proposed schemes. Some concluding remarks are presented in Section 6.

Preliminaries
Before constructing high-order schemes that can preserve the inequalities in Definitions 1.1-1.4 at time steps larger than that allowed by SSP integrator, we first provide sufficient conditions for (1.1) to have corresponding structures. By extending the analytical framework developed by Du et al. [4], the following result can be derived. Strong stability : Range boundedness : Maximum principle : Here, the symbol ‖ · ‖ * can be any norm, e.g., the infinity-norm or the 2-norm. Then (1.1) has a unique solution ∈ (︀ [0, ]; R )︀ and admits the corresponding inequalities: Strong stability : Rangeboundedness : where the inequalities in (2.7) and (2.8) are component-wise, and ( ) in (2.10) is the solution to (1.1) with initial condition 0 .
The proofs of (2.6)-(2.10) essentially follow that of Theorem 2.3 in [4]. For the completeness, we present it in Appendix A.

Inequality-preserving integrators with a time-step-dependent parameter
To construct efficient, accurate and stable time integration methods for general stiff, nonlinear differential equations, we follow previous studies [57,[64][65][66][67] and set out three key design principles. First, fixed points of the system are preserved. Second, the nonlinear term is handled simply and inexpensively. Third, the time-step size is selected to reflect the accuracy requirement, rather than restricted by the stability constraint. In this section, we develop a family of inequality-preserving parametric schemes that satisfy the above principles and are remarkably easy to implement.
, and powers of vectors are component-wise.
Assume the initial condition 0 = * is a steady state to (3.1), i.e., ( , * ) = 0. It is expected that a numerical method will preserve the steady state.
An -stage, th-order explicit RK scheme for (3.1) can be presented in the Shu-Osher form [68], , = 1, and non-negative coefficients , , , are constrained by certain accuracy requirements [68,69]. The explicit Shu-Osher formulation can be written in the Butcher form with Butcher coefficients given by [14], For convenience, we present second-and third-order conditions [1] using Shu-Osher coefficients, and fourth-order conditions [70] using Butcher coefficients in Table 1.
For an SSP RK scheme, we assume a given , is zero only if its corresponding , is zero. If we denote the SSP coefficient by = min , To remove the restriction on the time-step size, let us apply the Lawson transformation [71] to the unknown ( ), i.e., ( ) = e ( ). Then we obtain an equivalent system of the stabilized formulation (3.1) Applying the Shu-Osher formulation (3.2) to problem (3.5) and using the inverse transformation ( ) = e − ( ) yield the Shu-Osher IFRK method: Note that the denominator e will not equal the numerator ∑︀ −1 =0 e ( , + , ) unless = 0. Therefore, IFRK (3.6) cannot preserve fixed points unless * is a zero vector. To solve this problem, let 0 ( ) := 1, we propose approximating the exponential functions e using Then we get the modification in the Shu-Osher form whereˆ, We denote (3.9) [equivalent to (3.8)] as the parametric RK method (pRK), and the parametric abscissas are Clearly, pRK preserves fixed points and are explicit.
By using the equivalence between the Shu-Osher IFRK formulation (3.6) and the following Butcher IFRK formulation we rewrite the parametric RK scheme in the Butcher form The parametric Butcher tableau of pRK is given bŷˆ= We then study the accuracy of pRK schemes. To prove the th order convergence of pRK, it is sufficient to prove consistency of order (see, e.g., Theorem 2.3.4 of [72]).  Table 2, then the pRK (3.9) as applied to (1.1) has a ( +1)th order truncation error, that is, Proof. To derive order conditions for the pRK scheme, we consider its Butcher form (3.12) on an ODE = ( , ), and apply the Taylor expansions to obtain 14) Substituting (3.14) and (3.15) into (3.16) gives the defects where , = ( , , , ) − ( , , ( , )) are the right-hand-side stage errors. We assume an expansion for , as a power series in , Denoting the "error factors" of the defects ∆ , (3.17) as , = ! − 1 Hence the pRK method is consistent of order if We then relate , recursively to , . Note that We obtain the expansion of , as Hence, the order conditions contain both andˆ. Let us definẽ By using the rooted trees theory [70], we present the second-to fourth-order conditions in Table 2. Noting that the Butcher modification can be transformed to the Shu-Osher modification, we also give the order conditions of pRK(2, 2) and pRK (3,3) in Table 2 using the parametric Shu-Osher coefficients.
In the error expansion, each of these residuals is multiplied by the ( ) term in the "Term" column of Table 2, then the overall truncation error is ( +1 ∑︀ =0 ), which is ( +1 ) after absorbing . For an -dimensional system of ODEs, the truncation error is derived using the same approach by transition to multidimensional linear operators.
Next, we consider the pRK(2, 2), pRK (3,3) and pRK (4,4) schemes to show that the pRK schemes retain the original convergence orders of underlying RK schemes. By applying conditions . . , , and the transformation (3.4), we compute the approximations to exponential functions e , = 1, 2, 3 using Shu-Osher coefficients and e 4 using Butcher coefficients as below.
For fixed value of , verification of second-order conditions for pRK(2, 2) is done as follows: Other conditions containing 1,0 in Table 2 can be verified similarly.
Verification of third-order conditions for pRK (3,3) is done as follows: Other conditions containing 1,0 and 2 in Table 2 can be verified similarly.
Verification of the fourth-order conditions with = 0 for pRK (4,4) is carried out using the Butcher coefficients and done as follows: For the conditionsˆ˜2 = 1 3 + ( 2 2 ), = 0, 1, 2 of pRK (4,4) and all other conditions, we verify them for each Butcher tableau appearing in this work in repository [73]. It shows that the pRK scheme retains the original convergence order of underlying RK scheme, but the truncation error depends on and . Thus the smaller is , the better is the accuracy. We point out that this is a common phenomenon of exponential integrators [74], although the exponential functions in IFRK are replaced by polynomial functions in pRK. Assume ‖ , ‖ ≤ ‖ ‖, = 0, . . . , − 1; by using (3.21) and applying the forward Euler condition (1.2) and definitions of ( ), we obtain Thus ‖ +1 ‖ = ‖ , ‖ ≤ ‖ ‖. Other properties in Definitions 1.1-1.4 can be proved similarly.
Note that there is no explicit four-stage, fourth-order SSP RK scheme with > 0 [2,14,75], and no explicit RK scheme higher than the fourth order has positive SSP coefficient [2,75]. The unique four-stage, fourth-order accurate explicit RK scheme with non-negative Butcher coefficients is the classical RK(4, 4) scheme with = 0 (p. 521 in [2]). Therefore, the explicit pRK approach can give at most fourth-order accurate inequality-preserving integrators for any positive time step.
We present some non-negative RK formulas up to the fourth order with SSP coefficients below.  Since we interpret the solution , of (3.12) as an approximation to ( + ) and the proposed approximation guarantees thatˆ= ), ∀ ≥ 0, the time inconsistency is inevitable unless = 0. Thus the essence of the stabilization method is to trade consistency for stability (inequality-preservation) when > FE , which is an effective and mainstream strategy in the construction of high-order stable schemes for stiff and nonlinear problems [4,8,46,49]. Nevertheless, the pRK can reach the fourth-order convergence in preserving inequality structures and the TSD stabilization = max {︁ }︁ makes good use of the original SSP coefficient and has no influence on the accuracy when ≤ FE , hence pRK successfully reduces the time inconsistency comparing with existing stabilization schemes that use a constant parameter ≥ 1 FE (presented in Sec. 3.3). In the literature, Ju et al. [39] pointed out that larger time step size leads to less accurate solutions in stabilization ETD-SAV schemes, and Li et al. [76] noted that stabilization would slow down the convergence rate. Although pRK is inequalitypreserving for any time step, the time step size should be restricted by the desired accuracy instead of being arbitrarily large. Therefore, we call the proposed pRK schemes as efficient inequality-preserving rather than unconditionally inequality-preserving.

Parametric methods for systems with both stiff and nonlinear components
Now let us re-consider the system with both stiff and nonlinear terms in the form = + ( , ), ∀ ∈ (0, ], (3.24) where the stiff constant coefficient linear term , ∈ R × and the nonlinear term ( , ) satisfy some forward Euler conditions, and the allowable time step for the forward Euler condition on the linear component can be significantly smaller than that for the nonlinear component. We first show that some forward Euler conditions on the linear term indicate that e can preserve corresponding inequalities for any > 0.
Lemma 3.5. Assume that the linear term , ∈ R satisfies one or more of the following forward Euler conditions: Positivity : Then, it holds that where the inequalities in Since can be arbitrarily large and ≥ 0, we have e ≥ 0, ∀ > 0. Other results can be proved in the similar way.

35)
Proof. Let˜= 1 ≥ 1 FE ,¯= 1 ≥ 1 FE · By multiplying (1.9) and (3.28) with˜and¯, respectively, we obtain Then we obtain which is equivalent to By using Lemma 3.7 and Theorem 3.3, we present the following corollary without proof: Note that the strong stiffness of the linear term usually forces the forward Euler time-step size¯F E to be small; thus, a large stabilization parameter will be needed to ensure the inequality-preserving properties. To reduce truncation errors introduced by a large , we consider taking advantage of properties of the operator e in (3.30)-(3.34).
By treating e − as the integrating factor, the SSP integrating factor approach for (3.24) has been extensively studied by Isherwood et al. [34][35][36]. Their studies showed that by applying the contraction semi-group property ‖e ‖ ≤ ‖ ‖, the integrating factor approach only required the time step to satisfy the forward Euler condition on the nonlinear term for preserving strong stability, i.e., ≤˜F E , which greatly relaxed the requirements ≤˜F E¯FẼ FE+¯FE on explicit SSP schemes. By incorporating the integrating factor method (integrating factor being e − ) to (3.24) with underlying parametric single-step schemes pRK (3.9), we obtain Although pIFRK may not preserve fixed points of (3.1) because of the introduction of the integrating factor, as compensation, we have the preservation of inequalities with smaller for pIFRK. Proof. Consider the preservation of strong stability as an example. We prove this by induction. Assuming that ‖ , ‖ ≤ ‖ ‖, = 0, . . . , − 1, noting that , , + , ≤˜F E , by applying the forward Euler condition on ( , ), the contraction semi-group of (3.30) and definitions of ( ) ≥ 0, ∀ ≥ 0, we can derive This completes the proof.
In Theorem 3.9, a critical restriction on underlying RK schemes is the non-decreasing abscissas. Such explicit SSP RK schemes have been studied by Isherwood et al. [34], and shown to be at most fourth-order. We present some third-and fourth-order SSP RK schemes with non-decreasing abscissas (denoted by a superscript " + ") and optimal SSP coefficients (w.r.t. the stage number) below. We mention that, when = 0, the pIFRK schemes reduce to the IFRK schemes developed in [34]. ]︀ :

Inequality-preserving exponential time difference schemes
Noting that the existing higher-than-first-order unconditionally SSP [3,[23][24][25] and positivity-preserving [26,28] integrators mentioned in the introduction are specifically designed for particular problems, and have not been proven to preserve inequality structures for general problems. For comparison, we present the well-known ETD integrators in this section. Ostermann and co-authors [77,78] proved that the order of positive ETD Runge-Kutta and multi-step schemes cannot exceed two. Therefore, we only show that the ETD1 and ETD2 schemes preserve inequality structures for any time step size by adopting a stabilization technique.
The proofs for preservation of other structures in Definitions 1.2-1.4 can be performed similarly.
Remark 3.13. The preservation of maximum principles for Allen-Cahn-type parabolic equations using stabilization ETD1/2 has been proven by Du et al. [4,38]. The results there demonstrate the efficiency and stability of stabilization methods. When considering the system (1.1), the corresponding stabilization ETD1/2 schemes have the form Since the proof of inequality-preservation can be similarly carried out as Theorem 3.12 when ≥ 1 FE for any > 0, we omit the details here. Because of the introduced stabilization parameter, the time inconsistency of ETD1/2 also exists. This issue will be illustrated in the numerical experiments. Proof. We prove this by mathematical induction. The first stage solution is calculated by
To analyze the stability of proposed schemes, we present a useful lemma. Then a property of proposed schemes (4.2) is given as follows: , it follows by symmetry that ⃒ ⃒ (︀ 1 2 , i )︀⃒ ⃒ = 1 for all real . Using the fact that the roots of the polynomial (− ), ≤ 4 satisfy Re( ) > 0, = 1, . . . , . Then As increases from 0 to 1.25, we obtain the boundaries of the stability regions for the proposed pRK schemes as | ( , )| = 1. The stability boundary curves for some pRK schemes are presented in Figure 1. Clearly, we observe that for each scheme, the stability region becomes larger as increases. When equals zero, the stability regions of pRK reduce to those of classical RK schemes. For = 1, 2, when ≥ 1 2 , the stability regions contain the whole left-half complex plane, see Figure 1(a) and (e). For s = 3-5, when > 1 2 , the stability regions exclude part of the left-half complex plain, see Figure 1(f, g, d, h). We remark that a larger may increase the stability but reduce the accuracy, as noted by Ju et al. [62] for stabilization integrating factor methods.

Numerical experiments
In this section, we underline our theoretical analysis with numerical experiments. For the convenience of discussion, we consider the preservation of strong stability for linear advection equations and the maximum principle for AC-type semilinear parabolic equations. First, we test the proposed pRK and pIFRK schemes to verify the convergence. Then we study the inequality-preservation of proposed schemes in terms of their allowable stabilization parameter for a given time-step size . When adopting parametric integrating factor and ETD1/2 schemes, we compute the products of matrix exponentials with vectors via the fast Fourier transform to accelerate computations [38].

Convergence order verification
We study the convergence of the proposed pRK (3.9) and pIFRK (3.36) schemes on ODE systems obtained from spatial-discretizations of a linear advection equation (5.1) and a nonlinear AC equation (5.3), respectively. where¯is the reference solution computed with a refined time-step size, and is the numerical solution obtained by pRK, pIFRK or ETD1/2. For the temporal convergence, the grid number is set to = 256, and the reference solution is computed with = 2 −18 using the fourth-order RK(5, 4) scheme (3.23). The time-step size required by the forward Euler condition with respect to the total variation is ≤ FE = ℎ | | = ℎ· Since the theoretical convergence orders are verified using a fixed value of , we choose the stabilization parameter = 1 FE that satisfies the requirement (3.20) for all time-step sizes. Figure 2(a) presents the numerical errors for the linear advection equation (5.1) at = 1 obtained by pRK schemes with underlying RK(1, 1), RK(2, 2), RK (3,3), RK (4,4), and RK(5, 4) Butcher tableaux in Section 3.1, and ETD1/2 in Section 3.3. Clearly, we observe that each of the proposed schemes converges with the theoretical order of accuracy. Although pRK(1, 1) and pRK (2,2) are less accurate than ETD1 and ETD2, respectively, the numerical errors can be greatly reduced by using the high-order pRK (3,3), pRK (4,4) and pRK (5,4) schemes. This shows the advantages of using high-order pRK schemes. with periodic boundary conditions. By discretizing the Laplace operator using the standard second-order central finite difference method, and denoting the obtained differentiation matrix by 2 , we obtain a semi-discrete system in the form = + ( ), where is given by 2 2 , ( ) = − 3 , and the powers of vectors are component-wise. Since the central finite difference discretization of the Laplace operator in one-, two-and three-dimensional domains with periodic boundary conditions produces a diagonally dominant matrix with negative diagonal entries, it can be verified that satisfies the condition that ‖e ‖ ∞ ≤ 1.
We set = 0.01 and solve this problem to time = 2.0. For the temporal convergence, the grid number is set to = 256, and the reference solution in (5.2) is computed with = 2 −12 using the fourth-order pIFRK + (5, 4) scheme (3.37) with = 0. Note that ( ) satisfies the forward Euler condition (1.9) for = 1 and˜F E = 1 2 (Theo. 1 in [31]). Thus we take two different stabilization parameters, = 1 and = 1 FE = 2, to show the influence of the stabilization parameter.
The numerical errors (5.2) computed using ETD1/2, and pIFRK with underlying RK(1, 1), RK(2, 2), RK + (3, 3), RK (4,4) and RK + (5,4) are demonstrated in Figure 2(b) and (c). As expected, we observe perfect firstto fourth-order convergence in the time direction for corresponding schemes, and the numerical errors do not blow up at large time step = 1. As increases from 1 to 2, the numerical errors become larger. Therefore, to reduce additional splitting errors, it is important to develop high-order schemes that can be used with small stabilization parameters. Example 5.3. We test the strong-stability preservation of the proposed schemes on a linear advection equation with discontinuous initial data + + = 0, ∈ (0, 1),

Strong-stability-preserving tests
and periodic boundary conditions.
We denote the semi-discrete system obtained by the upwind scheme as = + ( ), where = − and ( ) = − . We set = 1000,˜F E = 1 1000 , = 1 FE for ETD1/2, = max }︁ for pIFRK, and solve this problem for 100 time steps using existing IFRK [34], ETD1/2, and the proposed pIFRK schemes (3.36). When = 0, the IFRK and pIFRK schemes become the RK and pRK schemes, respectively. The results of the maximum rise in TV obtained by second-to fourth-order RK-type schemes with = 0, 5, 10 are presented in Figure 3. It can be observed that the allowable time-step size of each IFRK scheme becomes larger as increases. This is caused by the significant damping of exponential terms [34], thus reducing the oscillations and associated rise in total variation. The results computed by pRK, pIFRK and ETD1/2 with the prescribed values of are noteworthy. As an improvement over the existing RK and IFRK schemes, we can clearly observe that all pRK, pIFRK and ETD1/2 schemes (dashed lines with markers in Fig. 3) preserve the strong stability for any value of by introducing the stabilizing parameter.
Next, we compare the solutions of pIFRK schemes with those obtained by IFRK schemes [34], the unconditionally SSP backward Euler (BE) scheme, and ETD1/2 schemes at large time steps. In Figure 4, we present solution profiles computed by different schemes with = 0.1 at = 0.18. In Figure 4 = 0.00151. The corresponding ratios˜F E are: 1.13, 1.355, 1.51. Since the time steps are larger than˜F E for all IFRK schemes, oscillations appear in the solution profiles obtained by these schemes. With the introduction of stabilization parameters, ETD1/2 and pIFRK schemes produce stable profiles, and the time delays of pIFRK schemes are weaker than those of ETD1/2. Moreover, the solutions of pIFRK(2, 2) and pIFRK + (3, 3) are more accurate than that of BE. Noting that IFRK (4,4) has SSP coefficient = 0 and can not guarantee the strong stability for any > 0, thus the time delay of pIFRK (4,4) is obvious because the stabilization parameter is  larger than those for pIFRK(2, 2) and pIFRK + (3,3). This is what we call stabilization trades consistency for stability.
In Figure 4(b), we choose a critical large time-step size = 0.0024 >˜F E for all IFRK schemes. It demonstrates that the BE scheme introduces severe smoothing effect to the solution because of its low order, and the standard IFRK + (5, 4) generates oscillations for this time step. Since the solutions of IFRK(2, 2), IFRK + (3, 3) and IFRK (4,4) blow up at this time step size, we do not present their solutions. On the contrary, the pIFRK(2, 2) scheme with a proper makes the solution stable, but large time delay appears because of introduced truncation errors. With the increasing of temporal convergence order, the solution delay of pIFRK becomes weaker, and the fourth-order pIFRK + (5, 4) provides more accurate solution than the backward Euler scheme. This shows the superiority of the fourth-order pIFRK + (5, 4) scheme over lower-order schemes, which is one main  contribution of this work. If we further increase the time step to a very large value, we point out that pIFRK + (5, 4) will suffer from stronger time delay and may become less accurate than BE method, thus the time step of parametric schemes should be restricted by the prescribed accuracy. To further reduce the time delay at large time steps, the fourth-order pIFRK method with more stages [34] could be used. The ratios of this time step and the allowable time steps for the pIFRK schemes, as well as the CPU time costs are presented in Table 3.
Because of the explicitness, the pIFRK schemes are more efficient than the backward Euler method. This is what we mean the pIFRK schemes are efficient inequality-preserving.  where (·, ·) denotes the 2 inner product, i.e., ( , ) = ∫︀ Ω d , and the energy functional is given by

Maximum-principle-preserving tests
We test the third-and fourth-order integrators on the second-order central finite difference semi-discrete system. For a -dimensional space, we denote the 2 -inner product by ⟨ , ⟩ := ℎ , 1⟩. The solution profiles, evolutions of ‖ ‖ ∞ , and energy profiles are presented in Figure 5. Figure 5 column (a) presents the results of the AC equation computed with underlying RK + (3, 3) scheme. We can observe that when the time step size increases from 0.1 to 1.5, IFRK + (3, 3) produces oscillations in the solution (first row) and energy (third row) profiles, finally breaking the maximum principle (second row). By introducing a stabilization term with = 1 FE − = 1.5, the solution computed by pIFRK + (3, 3) with = 1.5 is not only accurate but also preserves the maximum-principle, and dissipates the energy. When adopting the fourth-order underlying RK(4, 4) scheme, although Figure 5 column (b) shows that the solution computed using IFRK(4, 4) with = 2.0 does not break the maximum principle, the solution profile is no longer accurate (first row), and oscillations appear in the energy profile. With the introduction of = 2.0, pIFRK (4,4) gives an accurate solution and well preserves such inequality structures. The results computed by the fourth-order schemes with underlying RK + (5, 4) Butcher tableau are presented in column (c) of Figure 5. Similar as those in column (b), the solution of IFRK + (5, 4) with = 3.0 has oscillations. With suitable approximations of the exponential functions, we can see that the pIFRK + (5, 4) scheme with = 3.0 and = 1 FE − ≈ 1.55 gives an accurate solution, preserves the maximum principle, and dissipates the energy. In addition, there is little difference between the computed solution and the reference solution obtained with = 0.1.   with periodic boundary conditions. Here and below, the rand() function is uniformly distributed in the interval (0, 1).
By choosing = 0.01, and discretizing the domain using = 128 2 grid points. We solve the AC equation to final time = 400 using the fourth-order schemes with underlying RK + (5, 4) Butcher tableau. Figure 6 presents the initial profile and the zero-level set snapshots at = 0, 15, 90, 240, 390 computed by IFRK + (5, 4) and pIFRK + (5, 4). The evolutions of ‖ ‖ ∞ and discrete energy are presented in Figure 7. As increases from 0.1 to 3.0, oscillations appear in the solution profile (second row of Figs. 6 and 7(a)), and the energy profile deviates from the reference one computed with = 0.1. Contrarily, by introducing the stabilization term with ≈ 1.55, the pIFRK + (5, 4) scheme with = 3.0 preserves the maximum principle in Figure 7(a) and dissipates the energy in Figure 7(b). In Figure 6, the pIFRK + (5, 4) solutions (bottom row) closely follow the reference ones in the top row. This clearly demonstrates the superiority of the proposed pIFRK schemes.

Concluding remarks
In this paper we proposed a unified approach to develop up to the fourth-order inequality-preserving schemes for differential equations satisfying forward Euler conditions. One distinctive feature of the proposed schemes is the treatment of the stabilization term. By introducing a time-step-dependent stabilization parameter and taking advantage of the integrating factor method, we proposed approximating the exponential terms containing the stabilization parameter using novel recurrent approximations. Consequently, high-order fixed-point-preserving schemes were constructed. The presented order conditions showed that the pRK schemes retain the original convergence orders of underlying RK schemes used in this paper. Together with forward Euler conditions, we derived that the proposed schemes preserve the corresponding inequalities when ≥ max {︁ This restriction on is weaker than those required in existing ETD1/2 schemes, i.e., ≥ 1 FE , thus reducing additional truncation errors. To remove the requirement of a large stabilization parameter caused by the stiff linear term, we further developed pIFRK schemes by introducing the stiff operator in exponential functions. Linear stability analysis showed that proposed pRK schemes have good stability when suitable parameters are chosen. Numerical experiments on some typical advection equations and the Allen-Cahn equation were also carried out in this work, which demonstrated the robustness of proposed schemes, and verified the theoretical preservation of strong stability and maximum principle.
Because of the particularity of selected problems, the introduced stabilization term is of the form ( − ), it is also meaningful to investigate the pRK approach for general nonlinear problems [58], e.g., the mean curvature equation [57] or the Cahn-Hilliard equation [54], where different stabilization terms are introduced to allow large time-step sizes. Furthermore, as has been recognized in the references [8,46,48,63,85], the technique of adding and subtracting a stabilization term − , > 0 to and from a system indeed trades consistency for stability by rescaling the time steps, and consequently time-delay appears as a by-product. Since it is a nontrivial task to remove the rescaling effects or estimate the correct rescaled time steps of pIFRK and ETD1/2, we will deal with such issues in a future work. The properties in (2.7)-(2.10) can be proved similarly.