High order strong stability preserving multi-derivative implicit and IMEX Runge--Kutta methods with asymptotic preserving properties
HHIGH ORDER POSITIVITY PRESERVING AND ASYMPTOTICPRESERVING MULTI-DERIVATIVE METHODS
SIGAL GOTTLIEB ∗ , ZACHARY J. GRANT † , JINGWEI HU ‡ , AND
RUIWEN SHU § Abstract.
In this work we present multi-derivative implicit-explicit (IMEX) Runge–Kuttaschemes. We derive their order conditions up to third order, and show that such methods canpreserve positivity (and more generally strong stability) with a time-step restriction independent ofthe stiff term, under mild assumptions on the operators. We present sufficient conditions under whichsuch methods are positivity preserving and asymptotic preserving (AP) when applied to a range ofproblems, including a hyperbolic relaxation system, the Broadwell model, and the Bhatnagar-Gross-Krook (BGK) kinetic equation. Previous efforts to devise such methods [12, 11] have used an IMEXRunge–Kutta framework plus a second derivative final correction. In this work, we extend thisapproach to include derivative information at any stage of the computation. This multi-derivativeIMEX approach allowed us to find a second order AP and positivity preserving method that improvesupon the one in [11] in terms of the allowable time-step size. Furthermore, this approach producesa third order method that is AP and positivity preserving for a time-step independent of the stiffterm, a feature not possessed by any of the existing third-order IMEX schemes. We present numericalresults to support the theoretical results, on a variety of problems.
1. Introduction.
In this work, we are interested in the development of a highorder time integrator for the ODE system of the form:(1) d u d t = T ( u ) + 1 ε Q ( u ) , where the solution u ( t ) ∈ R N , and the operators T , Q : R N → R N , N ≥
2. Theparameter 0 < ε ≤ O (1) indicates the regime of the problem: ε = O (1) correspondsto the non-stiff regime; ε (cid:28) Very often the operator Q satisfies the following prop-erties: Q is “conservative” in the sense that there exists a linear operator R : R N → R n , n < N , s.t. R Q ( u ) = 0, ∀ u ; Q is dissipative and has a unique local equilibrium of theform E ( R u ), where E : R n → R N is some operator. Using these properties, applying ∗ Mathematics Department, University of Massachusetts Dartmouth, North Dartmouth, MA02747. Email: [email protected]. SG’s research was supported in part by AFOSR Grant No.FA9550-18-1-0383. † Multiscale Methods and Dynamics Group, Mathematics in Computation Section, Oak RidgeNational Laboratory, Oak Ridge, TN 37830. Email: [email protected]. This manuscript has beenauthored by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department ofEnergy (DOE). The US government retains and the publisher, by accepting the article for publication,acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide licenseto publish or reproduce the published form of this manuscript, or allow others to do so, for USgovernment purposes. DOE will provide public access to these results of federally sponsored researchin accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan). ‡ Department of Mathematics, Purdue University, West Lafayette, IN 47907. Email: [email protected]. JH’s research was supported in part by NSF CAREER grant DMS-1654152. § Department of Mathematics, University of Maryland, College Park, MD 20742. Email:[email protected]. 1 a r X i v : . [ m a t h . NA ] F e b to (1) yields(2) d ω d t = R T ( u ) , ω := R u, which is not a closed system. However, if ε →
0, (1) implies Q ( u ) →
0, hence u → E ( ω ). Substituting this u into (2) gives a closed (reduced) system:(3) d ω d t = R T ( E ( ω )) . The above simple analysis reveals that when ε →
0, (1) is not only stiff but alsopossesses a non-trivial asymptotic limit . (Recall that n < N and note that the originalvariable u is in R N while the reduced variable ω is in R n .)Systems of the form (1) arise (after the method of lines discretization of a PDE)from many physical problems in multi-scale modeling. A prominent example is theBoltzmann equation in kinetic theory [4]:(4) ∂ t f + v · ∇ x f = 1 ε Q ( f ) , x, v ∈ R d , where f = f ( t, x, v ) ≥ t , position x , andvelocity v . The term v · ∇ x f describes the particle transport, and Q ( f ) describes thecollisions between particles, which is a complicated nonlinear integral operator. Thedimensionless parameter ε , called the Knudsen number, is defined as the ratio of themean free path and characteristic length scale. When ε = O (1), the transport andcollision balance so the system is in the fully kinetic regime. When ε (cid:28)
1, the collisioneffect dominates, i.e., collisions happen so frequently that the overall system is closeto the local equilibrium or fluid regime. In this case, one can derive the limiting fluidequations (the compressible Euler equations) as ε → t , when ε →
0, the scheme for (1) automaticallyreduces to a high order time discretization for the limiting system (3). A numericalscheme with this property is called asymptotic preserving (AP) as initially coinedin [13]. To insure the AP property, the time step ∆ t should not be limited by thesmall parameter ε . This necessitates some implicit treatment of the stiff collision term ε Q ( u ). For this reason, we will consider implicit-explicit (IMEX) methods. Problems of the form (1) that are of interest to us have positivesolutions. It is preferable that the numerical solution will preserve this positivityproperty, for a time-step not dependent on ε . However, high order implicit treatmentof the collision term, as done in the IMEX methods, usually causes the numericalsolution to lose this ε -independent positivity property [6, 7].It should be pointed out that positivity is an important property when solvingkinetic equations. For example, the Bhatnagar-Gross-Krook (BGK) model (see equa-tion (23) below) requires the macroscopic quantities to be positive, and even smallnegative values of the solution f may cause some macroscopic quantities, especiallythe temperature, to fail to be well-defined. In such cases, the requirement that thenumerical solution remains positive for time-steps independent of ε is critical to thesuccess of the simulation. .3. Overview of this paper. The need for a high order numerical integratorthat is both asymptotic preserving and positivity preserving motivates the work in thispaper. Specifically, we will design asymptotic preserving high order time discretiza-tion methods (second and third order) that preserve the positivity of the solution forarbitrary ε . Designing a time-stepping scheme with both positivity and AP prop-erty has proven difficult. First order IMEX schemes with these properties exist, butmethods above first order may violate positivity unless the time-step is restricted by ε [8, 9].Second order IMEX schemes that preserve the AP property and positivity forarbitrary ε have been found, by incorporating a derivative correction term of theform ˙ Q ( u ) := Q (cid:48) ( u ) Q ( u ) at the final stage of each time-step, where Q (cid:48) is the Fr´echetderivative of Q . Such an approach was successfully considered in [12, 11] (note thatthe method in [12] only works for a special relaxation system and can preserve thepositivity of one component of the solution vector, while [11] works for a general classof equations and the scope is similar to what we consider in this work); however, thisstrategy failed to find methods of order three. By allowing the use of ˙ Q at every stage(not the just the final stage), we are able to obtain a third order IMEX method thatis AP and positivity preserving independent of ε .We begin this paper with a summary of the model equations and their propertiesin Section 2. We then present the multi-derivative IMEX approach in Section 3,and present sufficient conditions under which we can ensure the AP and positivitypreserving property of the method under a time-step independent of ε . We then derivethe order conditions for two-derivative IMEX methods. In Section 4 we present ournew second and third order methods. The second order method improves upon thepreviously presented method in [11], in the sense that by allowing ˙ Q at every stagewe obtain a larger allowable time-step. The third order method is the first method ofthis order that is both AP and positivity preserving under the assumptions we use.In Section 5 we demonstrate the numerical performance of these methods on sampleproblems. The paper is concluded in Section 6.
2. A summary of the models and properties.
We assume that the operators T and Q in (1) satisfy the following properties: Property 1
The operator T is conditionally positivity preserving under a forward Euler step: (5) u > ⇒ u + ∆ tT ( u ) > , ∀ ≤ ∆ t ≤ ∆ t FE , for some time step ∆ t FE > . Property 2
The operator Q is unconditionally positivity preserving under a backward Eulerstep: (6) u > , v = u + ∆ tQ ( v ) = ⇒ v > , ∀ ∆ t ≥ . We observe that the first two properties essentially concern the positivity preservingproperty of the operators T and Q in equation (1). Property 3
Conservation of Q : there exists a linear operator R : R N → R n , n < N , s.t. (7) R Q ( u ) = 0 , ∀ u. roperty 4 Equilibrium of Q : there exists an (possibly nonlinear) operator E : R n → R N , s.t. (8) Q ( u ) = 0 = ⇒ u = E ( R u ) . Moreover, E satisfies R E ( R u ) = R u , ∀ u . Note that Properties 3 and 4 together imply that (1) has a limiting system (3) as shownin the Introduction. Properties 1–4 are satisfied by a large class of kinetic equationsof the form (4), where the collision operator Q can be the full Boltzmann collisionoperator (an integral type operator), the kinetic Fokker-Planck operator (a diffusiontype operator), the BGK operator (a relaxation type operator), or its generalizedversion such as the ES-BGK operator. For more details about these operators, werefer the readers to [10]. Property 5
The Fr´echet derivative of Q satisfies (9) ˙ Q ( u ) := Q (cid:48) ( u ) Q ( u ) = − C R u Q ( u ) , where C R u is some positive function depending only on R u . The Fr´echet deriva-tive of Q at u is defined by (10) Q (cid:48) ( u ) v = lim δ → Q ( u + δv ) − Q ( u ) δ . Property 5 means the operator Q is dissipative in some sense. This property is notgeneric but it is satisfied by quite a few kinetic models including the BGK operatorand the Broadwell model. Some stiff ODE systems and hyperbolic relaxation sys-tems also satisfy this property, though for these problems positivity is usually not abig concern compared to the kinetic equations. Since our proposed multi-derivativemethods highly depend on Property 5, we list below a few examples. • An ODE model: (11) u (cid:48) = u ,u (cid:48) = 1 ε f ( u ) ( g ( u ) − u ) , where f and g are some functions of u , and f ( u ) >
0. Define(12) u = ( u , u ) T , T ( u ) = ( u , T , Q ( u ) = (0 , f ( u ) ( g ( u ) − u )) T , then (11) falls into the general form (1). It is easy to see that (11) has a limitas ε → u (cid:48) = g ( u ) . Indeed, one can just take R := (1 ,
0) and E ( R u ) = E ( u ) := ( u , g ( u )) T . Itcan also be verified by direct calculation that(14) ˙ Q ( u ) = − f ( u ) Q ( u ) . • A PDE model: the hyperbolic relaxation system [5]: (15) ∂ t u + ∂ x u = 0 ,∂ t u + ∂ x u = 1 ε ( F ( u ) − u ) , here F is some function of u . Equation (15) again has the form of (1) ifwe define u = ( u , u ) T , T ( u ) = − ( ∂ x u , ∂ x u ) T , Q ( u ) = (0 , F ( u ) − u ) T .Note that we abused the notation a bit: u , T and Q should be defined forthe system after spatial discretization. It is easy to see that (15) has a limitas ε → ∂ t u + ∂ x F ( u ) = 0 . Indeed, one can just take R := (1 ,
0) and E ( R u ) = E ( u ) := ( u , F ( u )) T .Similarly to the previous model, it can be verified that(17) ˙ Q ( u ) = − Q ( u ) . • The Broadwell model [3]:
The Broadwell model is a simple discrete ve-locity kinetic model:(18) ∂ t f + + ∂ x f + = 1 ε ( f − f + f − ) ,∂ t f = − ε ( f − f + f − ) ,∂ t f − − ∂ x f − = 1 ε ( f − f + f − ) , where f + = f + ( t, x ), f = f ( t, x ), and f − = f − ( t, x ) denote the densitiesof particles with speed 1, 0, and −
1, respectively. Define f = ( f + , f , f − ) T , T ( f ) = ( − ∂ x f + , , ∂ x f − ) T , and Q ( f ) = ( f − f + f − , − ( f − f + f − ) , f − f + f − ) T (again these should be defined for the system after spatial discretiza-tion). Then (18) falls into the general form (1). To see its limit as ε →
0, werewrite (18) using moment variables:(19) ∂ t ρ + ∂ x m = 0 ,∂ t m + ∂ x z = 0 ,∂ t z + ∂ x m = 12 ε ( ρ + m − ρz ) , where ρ := f + + 2 f + f − , m := f + − f − , and z := f + + f − . From (19), it isclear that when ε → z → ρ + m ρ . This, when substituted into the first twoequations, yields a closed hyperbolic system:(20) ∂ t ρ + ∂ x m = 0 ,∂ t m + ∂ x (cid:18) ρ + m ρ (cid:19) = 0 . Indeed, the operators R and E in Properties 3–4 can be taken as R := (cid:18) − (cid:19) , (21) E ( R f ) = E (( ρ, m ) T ) := (cid:18) ( ρ + m ) ρ , ρ − m ρ , ( ρ − m ) ρ (cid:19) T . Furthermore, it can be verified that(22) ˙ Q ( f ) = − ρQ ( f ) . The Bhatnagar-Gross-Krook (BGK) model [2]:
The BGK model is awidely used kinetic model introduced to mimic the full Boltzmann equation:(23) ∂ t f + v · ∇ x f = 1 ε ( M − f ) , x, v ∈ R d , where f = f ( t, x, v ) is the probability density function and M is the so-calledMaxwellian given by(24) M ( t, x, v ) = ρ ( t, x )(2 πT ( t, x )) d/ exp (cid:18) − | v − u ( t, x ) | T ( t, x ) (cid:19) , where the density ρ , bulk velocity u and temperature T are given by themoments of f :(25) ρ = (cid:90) R d f d v, ρu = (cid:90) R d f v d v, ρdT = 12 (cid:90) R d f | v − u | d v. To see its asymptotic limit, we multiply (23) by (1 , v, | v | / T and integratew.r.t. v to obtain ∂ t ρ + ∇ x · (cid:90) R d vf d v = 0 ,∂ t ( ρu ) + ∇ x · (cid:90) R d v ⊗ vf d v = 0 ,∂ t E + ∇ x · (cid:90) R d v | v | f d v = 0 , (26)where E = ρu + ρdT is the total energy. This system is not closed.However, if ε →
0, (23) implies f → M . Substituting this f into (26), we canget a closed system ∂ t ρ + ∇ x · ( ρu ) = 0 ,∂ t ( ρu ) + ∇ x · ( ρu ⊗ u + pI ) = 0 ,∂ t E + ∇ x · (( E + p ) u ) = 0 , (27)where I is the identity matrix and p = ρT is the pressure. Equation (27) isnothing but the compressible Euler equations. To write the BGK model intothe form (1), we define T ( f ) = − v · ∇ x f and Q ( f ) = M − f (these should bedefined for (23) after spatial and velocity discretization). Moreover, R is theoperator (cid:82) · (1 , v, | v | / T d v , and(28) E ( R f ) = E (( ρ, ρu, E ) T ) = M. Furthermore, it can be verified that(29) ˙ Q ( f ) = − Q ( f ) . To summarize, we have introduced four different models (including both ODE andPDEs) which all satisfy Properties 3–5. For the Broadwell model and BGK model,one can check that they also satisfy the positivity-preserving Properties 1–2 provideda positivity preserving spatial discretization is used for the transport/convection term,see [11] for more details. . Multi-derivative IMEX methods. For problem (1), we propose an s -stagemulti-derivative IMEX method, written in the Shu-Osher formulation, as follows u ( i ) = r i u n + i − (cid:88) j =1 p ij u ( j ) + i − (cid:88) j =1 w ij (cid:18) u ( j ) + ∆ tr T ( u ( j ) ) (cid:19) (30a) + ∆ tε d ii Q ( u ( i ) ) + ∆ t ε ˙ d ii ˙ Q ( u ( i ) ) , i = 1 , ..., s,u n +1 = u ( s ) . (30b)where the terms Q and ˙ Q appear only implicitly.The intermediate stages can be conveniently written in a matrix form: U = Re u n + P U + W (cid:18) U + ∆ tr T ( U ) (cid:19) + ∆ tε D Q ( U ) + ∆ t ε ˙D ˙ Q ( U ) , (31)where P , W , and R = I − P − W are s × s matrices, r i are the i th row sum of R , D and ˙D are s × s diagonal matrices, and e is a vector of ones. The numerical solution u n +1 is then given by the final element of the vector U .The Shu-Osher form allows us to easily observe the positivity preserving propertiesof the method, which will be established in the next subsection. In this section wepresent sufficient conditions under which the method (30) is AP and preserves thepositivity of the solution for time-steps that are independent of ε . Proposition Assume that the problem (1) satisfies the Properties 1–5 listedin Section 2. Then the scheme (30) that satisfies the inequalities (element-wise) Re ≥ , P ≥ , W ≥ , D ≥ , ˙ D ≤ , (32) will preserve the positivity of the solution for all ∆ t ≤ r ∆ t FE . Furthermore, if thestrict inequality d ii + (cid:12)(cid:12)(cid:12) ˙ d ii (cid:12)(cid:12)(cid:12) > , for all i = 1 , . . . , s, (33) is also satisfied, then the scheme is AP: when ∆ t is fixed, as ε → , (30) automaticallyreduces to an explicit Runge-Kutta scheme, with the same order as the original scheme,applied to the limiting system (3).Proof. We consider each stage of (30), u ( i ) = r i u n + i − (cid:88) j =1 p ij u ( j ) + i − (cid:88) j =1 w ij (cid:18) u ( j ) + ∆ tr T ( u ( j ) ) (cid:19) (34) + ∆ tε d ii Q ( u ( i ) ) + ∆ t ε ˙ d ii ˙ Q ( u ( i ) )= r i u n + i − (cid:88) j =1 p ij u ( j ) + i − (cid:88) j =1 w ij (cid:18) u ( j ) + ∆ tr T ( u ( j ) ) (cid:19) + (cid:18) ∆ tε d ii − ∆ t ε ˙ d ii C R u ( i ) (cid:19) Q ( u ( i ) ) , where we applied Property 5 to the last term ˙ Q ( u ( i ) ). t the first stage, we have u (1) = u n + (cid:18) ∆ tε d − ∆ t ε ˙ d C R u (1) (cid:19) Q ( u (1) ) . Given a positive u n , and since d ≥
0, ˙ d ≤
0, and C R u (1) >
0, using Property 2 weobtain u (1) > u n and positive stages u ( j ) for j < i , Property 1 gives usthe positivity of the explicit terms (cid:18) u ( j ) + ∆ tr T ( u ( j ) ) (cid:19) > , for all ∆ tr ≤ ∆ t FE . Consequently, the non-negativity of r i , p ij , and w ij , together with the fact that r i + (cid:80) i − j =1 ( p ij + w ij ) = 1 ensures the positivity of the explicit terms in u ( i ) r i u n + i − (cid:88) j =1 p ij u ( j ) + i − (cid:88) j =1 w ij (cid:18) u ( j ) + ∆ tr T ( u ( j ) ) (cid:19) > . Finally, since d ii ≥
0, ˙ d ii ≤
0, and C R u ( i ) >
0, Property 2 assures that u ( i ) > R to (34) to obtain (define ω n = R u n , ω ( i ) = R u ( i ) ) ω ( i ) = r i ω n + i − (cid:88) j =1 p ij ω ( j ) + i − (cid:88) j =1 w ij (cid:18) ω ( j ) + ∆ tr R T ( u ( j ) ) (cid:19) , (35)where the collision terms are gone due to Property 3. On the other hand, when ∆ t is fixed and ε →
0, since d ii + | ˙ d ii | > C R u ( i ) >
0, we have from (34) that Q ( u ( i ) ) →
0, hence u ( i ) → E ( ω ( i ) ) by Property 4. Note that this holds for every i = 1 , . . . , s . Replacing u ( j ) by E ( ω ( j ) ) in (35) yields ω ( i ) = r i ω n + i − (cid:88) j =1 p ij ω ( j ) + i − (cid:88) j =1 w ij (cid:18) ω ( j ) + ∆ tr R T ( E ( ω ( j ) )) (cid:19) , i = 1 , . . . , s ;together with ω n +1 = ω ( s ) , this is a high order explicit Runge-Kutta scheme appliedto the limiting system (3). In fact, it is the explicit part of (30) applied to (3). Remark The proof above focuses on preserving the strict positivity of the so-lution. Under more relaxed conditions, this proof can be adapted to preserve non-negativity, or other types of more general convex functional conditions. We addressthis in Appendix A.
Remark It is worth noting that many AP schemes require the initial data tobe consistent, i.e., one should start from the equilibrium u = E ( ω ) ; otherwise, therewould be an order reduction even in the limit ε → (for example, in the commonlyused ARS type IMEX Runge-Kutta schemes [1, 7]). In our scheme, since d ii + (cid:12)(cid:12)(cid:12) ˙ d ii (cid:12)(cid:12)(cid:12) > for all i , we are solving an implicit collision step at every stage of the scheme, henceany initial condition is allowed to guarantee the AP property. In Section 4 we show that it is indeed possible to find methods that satisfy therequirements in Proposition 1. We present second and third order methods that areprovably AP and positivity preserving. .2. Formulating order conditions. The order conditions for a method (30)are generally easier to formulate if the method is written in its Butcher form: u ( i ) = u n + ∆ t i − (cid:88) j =1 ˆ a ij T ( u ( j ) ) + ∆ tε i (cid:88) j =1 a ij Q ( u ( j ) ) + ∆ t ε i (cid:88) j =1 ˙ a ij ˙ Q ( u ( j ) ) , (36a) i = 1 , ..., s,u n +1 = u n + ∆ t i − (cid:88) j =1 ˆ b j T ( u ( j ) ) + ∆ tε i (cid:88) j =1 b j Q ( u ( j ) ) + ∆ t ε i (cid:88) j =1 ˙ b j ˙ Q ( u ( j ) ) . (36b)To be consistent with (30), we require that u n +1 = u ( s ) , so that ˆ b j = ˆ a sj , b j = a sj , ˙ b j = ˙ a sj . The intermediate stages of this method can be written in a matrixform:(37) U = e u n + ∆ t (cid:98) A T ( U ) + ∆ tε A Q ( U ) + ∆ t ε ˙ A ˙ Q ( U ) . The numerical solution u n +1 is then given by the final element of the vector U .The conversion between the two formulations (31) and (37) is given by: (cid:98) A = 1 r R − W , A = R − D , ˙ A = R − ˙ D . (38)The vectors (cid:98) b , b , and ˙ b are given by the last row of (cid:98) A , A , and ˙ A , respectively.The vectors c = Ae , ˙ c = ˙ Ae , and (cid:98) c = (cid:98) Ae define the time-levels at which the stagesare happening; these values are known as the abscissas. The order conditions formethods of this form are:For p ≥ b t e = 1 (cid:98) b t e = 1For p ≥ b t c + ˙ b t e = b t (cid:98) c = (cid:98) b t c = (cid:98) b t (cid:98) c = For p ≥ b t Ac + ˙ b t c + b t ˙ c = b t A (cid:98) c + ˙ b t (cid:98) c = b t (cid:98) Ac = b t (cid:98) A (cid:98) c = For p ≥ (cid:98) b t Ac + (cid:98) b t ˙ c = (continued) (cid:98) b t A (cid:98) c = (cid:98) b t (cid:98) Ac = (cid:98) b t (cid:98) A (cid:98) c = b t ( c · c ) + 2 ˙ b t c = b t ( c · (cid:98) c ) + ˙ b t (cid:98) c = b t ( (cid:98) c · (cid:98) c ) = (cid:98) b t ( c · c ) = (cid:98) b t ( c · (cid:98) c ) = (cid:98) b t ( (cid:98) c · (cid:98) c ) = Remark Note that these order conditions do not guarantee that we will notobserve order reduction. When ε = O (1) we expect to see the design accuracy predictedby the order conditions. When ε (cid:28) design accuracy may not be evident due to theorder reduction phenomenon. However, the AP property allows us to recover fullaccuracy in the asymptotic limit ε → .
4. New positivity preserving and asymptotic preserving methods. .1. Second order method. We begin with a method that has Shu-Oshercoefficients W = / , P = / , Re = , and diag ( D ) = , diag ( ˙ D ) = − , with r = 1.In Butcher form, these become (cid:98) A = / / , A = / / / / , ˙ A = − / − / . The benefit of this method over the one in [11] is that the positivity preserving co-efficient r = 1 for this method is larger than the positivity preserving coefficient r = 0 . We found a third order method of this form, aswell. This method has r = 0 . diag ( D ) = . . . , diag ( ˙ D ) = − . . . . Note that d ii + (cid:12)(cid:12)(cid:12) ˙ d ii (cid:12)(cid:12)(cid:12) > i . W = . . . . . . . , P = . . . . . . . , e = . . . We have the Butcher form coefficient matrices: ˆA = . . . . . . . . . . , A = . . . . . . . . . . . . . , and ˙A = − . . . . . . . . . . . . . . To achieve a third order method we required six stages and correction terms atmultiple stages (the first two and the final two). However, this allowed us to design athird order method that is both AP and positivity preserving.
5. Numerical results.
In this section, we verify the accuracy of the proposedsecond and third order methods in Section 4 on the ODE model, the Broadwell model,and the BGK equation. We will see that the methods exhibit the design accuracy inthe kinetic regime ε = O (1) as well as the fluid regime ε (cid:28)
1. This latter behavior isexactly due to the AP property of the methods. For completeness, we also report theresults of the methods in the intermediate regime (i.e., ε lies between 0 and 1), wherethe methods may exhibit some order reduction as expected. A careful study of thisbehavior is beyond the scope of the current work and left for future work. We consider the ODE model (11) with(39) f ( u ) = 1 + u , g ( u ) = sin u . We take the initial data as u (0) = (2 , T (which is an inconsistent initial data), andsolve (11) by the second and third order methods in Section 4, up to final time T = 1,with various ε and ∆ t . To calculate the error of a numerical solution U = [ U , U ] T ,we compare with a reference solution U ref obtained by the MATLAB solver ode15s with relative tolerance RelTol = 1 e −
13 and absolute tolerance
AbsTol = 1 e − | U ( T ) − U ref1 ( T ) | + | U ( T ) − U ref2 ( T ) | . he results are shown in Figures 1 and 2. For both methods, one can see the designorder accuracy in the kinetic regime ( ε = O (1) and ∆ t is relatively small) and thefluid regime ( ε (cid:28) t is not very small), while in the intermediate regime (when ε and ∆ t are comparable) one can see some order reduction. In Figure 2 with ε = 1(and similar for ε = 0 . , e − t decreaseswhen ∆ t is less than 5 e −
5, and this is a consequence of the accumulation of round-offerrors. -7 -6 -5 -4 -3 -2 -1 t -12 -10 -8 -6 -4 -2 e rr o r =1=0.01=0.0001=1e-06=1e-08=1e-10slope=2 Fig. 1 . Accuracy test of the new secondorder IMEX scheme for an ODE model. -7 -6 -5 -4 -3 -2 -1 t -12 -10 -8 -6 -4 -2 e rr o r =1=0.01=0.0001=1e-06=1e-08=1e-10slope=3 Fig. 2 . Accuracy test of the new third orderIMEX scheme for an ODE model.
We consider the Broadwell model (18) on thedomain x ∈ [0 ,
2] with periodic boundary condition, with inconsistent initial data f + (0 , · ) = 1 + 0 . . πx )) , f − (0 , · ) = exp(0 . πx )) ,f (0 , · ) = 11 + 0 . πx ) . We discretize in space by the fifth order finite volume WENO scheme, and the collisionoperator Q is evaluated pointwise on the Gauss quadrature points in each cell, asdescribed in Section 3.3.2 of [11]. We fix the the CFL number as ∆ t = ∆ x , andsolve (18) by the second and third order methods in Section 4 up to final time T = 0 . L norm of the difference between the numerical solutionand one with a refined mesh.The results are shown in Figures 3 and 4, and one can see similar behavior as inthe previous subsection. We consider the 1D BGK model (23) on the physicaldomain x ∈ [0 ,
2] with periodic boundary condition, and inconsistent initial data givenby(41) f (0 , x, v ) = 0 . M [˜ ρ ( x ) , ˜ u ( x ) , ˜ T ( x )]( v ) + 0 . M [˜ ρ ( x ) , − . u ( x ) , ˜ T ( x )]( v ) , with(42) ˜ ρ ( x ) = 1 + 0 . πx ) , ˜ u ( x ) = 1 , ˜ T ( x ) = 11 + 0 . πx ) . The velocity domain is truncated into [ − v max , v max ] with v max = 15 and discretizedwith N v = 150 grid points, and the physical space is discretized in the same way as -4 -3 -2 -1 t -12 -10 -8 -6 -4 -2 e rr o r =1=0.01=0.0001=1e-06=1e-08=1e-10slope=2 Fig. 3 . Accuracy test of the new secondorder IMEX scheme for the Broadwell model. -4 -3 -2 -1 t -12 -10 -8 -6 -4 -2 e rr o r =1=0.01=0.0001=1e-06=1e-08=1e-10slope=3 Fig. 4 . Accuracy test of the new third orderIMEX scheme for the Broadwell model. the previous subsection. We fix the the CFL number as ∆ t =
12 ∆ xv max , and solve (23)by the second and third order methods in Section 4 up to final time T = 0 .
1. Theerror is computed by the L norm (in the ( x, v ) space) of the difference between thenumerical solution and one with a refined mesh. -4 -3 -2 t -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 e rr o r =1=0.01=0.0001=1e-06=1e-08=1e-10slope=2 Fig. 5 . Accuracy test of the new secondorder IMEX scheme for the BGK model. -4 -3 -2 t -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 e rr o r =1=0.01=0.0001=1e-06=1e-08=1e-10slope=3 Fig. 6 . Accuracy test of the new third orderIMEX scheme for the BGK model.
The results are shown in Figures 5 and 6. For the second order scheme, one cansee clearly the second order accuracy when ∆ t is small enough (so that the temporalerror dominates) for both ε = O (1) and ε (cid:28)
1, and order reduction is observed in theintermediate regime. For the third order scheme, when ε = O (1) or ε (cid:28)
1, the errorconverges at a higher than expected rate even for the smallest ∆ t in the simulation,which suggests that the spatial error is still dominating. By comparing with theresults of the second order scheme we see that the third order scheme indeed gives amuch smaller error under the same time-step size.
6. Conclusions.
In this work, we extended the work in [11] which derived pos-itivity preserving and asymptotic preserving IMEX Runge–Kutta methods with aderivative correction term. We presented a class of multi-derivative implicit-explicit(IMEX) Runge–Kutta schemes that include derivative information at any stage ofthe computation, and derived order conditions up to third order. The appearanceof the derivative term at any stage allows us to obtain order p > r i , w ij , p ij and d ii are all non-negative, ˙ d ii is non-positive, and ither d ii or ˙ d ii must be nonzero at every stage u ( i ) . These conditions ensure that themethod is positivity preserving and asymptotic preserving when applied to problemsthat satisfy the five properties in Section 2. Examples of equations that satisfy theseproperties include a hyperbolic relaxation model, the Broadwell model, and the BGKkinetic equation.The multi-derivative IMEX Runge–Kutta approach we developed in this workallowed us to find a second order asymptotic preserving and positivity preservingmethod that improves upon the one in [11] in the sense that ∆ t ≤ r ∆ t FE where thestep size coefficient is increased from r = 0 . r = 1, a 23% relative increasein efficiency. Furthermore, we were able to obtain a third order method that isasymptotic preserving and positivity preserving with a time step restriction that isindependent of ε . These properties are not satisfied by any of the existing third-orderIMEX schemes.It is important to note that while the results in this paper focus on the positiv-ity preserving and asymptotic preserving properties, the methods we devised are ofbroader use. In particular, when r i , w ij , p ij and d ij are all non-negative, and ˙ d ij isnon-positive, these methods are strong stability preserving (SSP) for a broader classof equations, as we show in the appendix. Appendix A. Strong stability preservation.
The methods we presented inSection 4 are asymptotic preserving and positivity preserving. It is important to notethat under certain mild assumptions, such methods are also strong stability preserv-ing (SSP) and the time-step restriction depends only on the term that is handledexplicitly. Thus, these methods apply to a much broader class of problems. Thesemulti-derivative IMEX methods are significant in that they are the first such meth-ods of order p > u t = F ( u ) + G ( u ) where F and G preservesome nonlinear stability propeties under a convex functional (cid:107) · (cid:107) : Condition 1: (cid:107) u + ∆ tF ( u ) (cid:107) ≤ (cid:107) u (cid:107) for all ∆ t ≤ ∆ t FE , for some ∆ t FE >
0, and
Condition 2: (cid:107) u + ∆ tG ( u ) (cid:107) ≤ (cid:107) u (cid:107) for all ∆ t ≤ k ∆ t FE . Generally, we consider the case where the time-step restriction ∆ t FE coming from F is of a reasonable size, but the time-step restriction coming from G is very small(0 < k (cid:28) G . However, when we consider more general norms, semi-norms, orconvex functionals, the use of IMEX schemes does not result in the removal of thetime step restriction caused by the operator G . However, if G satisfies an additionalcondition Condition 3: (cid:107) u − ∆ t ˙ G ( u ) (cid:107) ≤ (cid:107) u (cid:107) for all ∆ t ≤ ˙ k ∆ t , (where ˙ k > G . Thefollowing Theorem formally states this result. heorem Given operators F and G that satisfy Conditions 1, 2, and 3, withvalues ∆ t FE > , k > , ˙ k > , for some convex functional (cid:107) · (cid:107) , and if the methodgiven by (31) with r > satisfies the conditions Re ≥ , P ≥ , W ≥ , D ≥ , ˙ D ≤ , (43) then it preserves the strong stability property (cid:107) u n +1 (cid:107) ≤ (cid:107) u n (cid:107) under the time-step condition ∆ t ≤ r ∆ t FE . Proof.
Each stage of the method is u ( i ) = r i u n + i − (cid:88) j =1 p ij u ( j ) + i − (cid:88) j =1 w ij (cid:18) u ( j ) + ∆ tr F ( u ( j ) ) (cid:19) + (cid:16) ∆ td ii G ( u ( i ) ) + ∆ t ˙ d ii ˙ G ( u ( i ) ) (cid:17) . First, let’s consider the implicit terms, u ( i ) = u + ∆ td ii G ( u ( i ) ) + ∆ t ˙ d ii ˙ G ( u ( i ) ) . Using Conditions 2 and 3 we can show that (cid:107) u ( i ) (cid:107) ≤ (cid:107) u (cid:107) , whenever d ii ≥ d ii ≤
0. To see this add ( α + β ) u ( i ) to both sides and rearrange u ( i ) = u α + β + α α + β (cid:18) u ( i ) + 1 α ∆ td ii G ( u ( i ) ) (cid:19) + β α + β (cid:18) u ( i ) − β ∆ t | ˙ d ii | ˙ G ( u ( i ) ) (cid:19) . Assuming that α ≥ β ≥ (cid:107) u ( i ) (cid:107) ≤
11 + α + β (cid:107) u (cid:107) + α α + β (cid:13)(cid:13)(cid:13) u ( i ) (cid:13)(cid:13)(cid:13) + β α + β (cid:13)(cid:13)(cid:13) u ( i ) (cid:13)(cid:13)(cid:13) , hence (cid:107) u ( i ) (cid:107) ≤ (cid:107) u (cid:107) , for any ∆ t such that α ∆ td ii ≤ k ∆ t FE and β | ˙ d ii | ∆ t ≤ ˙ k ∆ t . Since we can choose α and β to be arbitrarily large then this is true for any ∆ t .Now, we let u be the explicit part that depends only on F , u = r i u n + i − (cid:88) j =1 p ij u ( j ) + i − (cid:88) j =1 w ij (cid:18) u ( j ) + ∆ tr F ( u ( j ) ) (cid:19) . Given the non-negativity of all the coefficients (43) we have (cid:107) u (cid:107) = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) r i u n + i − (cid:88) j =1 p ij u ( j ) + i − (cid:88) j =1 w ij (cid:18) u ( j ) + ∆ tr F ( u ( j ) ) (cid:19)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) ≤ r i (cid:107) u n (cid:107) + i − (cid:88) j =1 p ij (cid:107) u ( j ) (cid:107) + i − (cid:88) j =1 w ij (cid:13)(cid:13)(cid:13)(cid:13) u ( j ) + ∆ tr F ( u ( j ) ) (cid:13)(cid:13)(cid:13)(cid:13) ≤ r i (cid:107) u n (cid:107) + i − (cid:88) j =1 p ij (cid:107) u ( j ) (cid:107) + i − (cid:88) j =1 w ij (cid:13)(cid:13)(cid:13) u ( j ) (cid:13)(cid:13)(cid:13) , or all ∆ t ≤ r ∆ t FE . Now, assuming that (cid:107) u ( j ) (cid:107) ≤ (cid:107) u n (cid:107) for j < i , we obtain (cid:107) u (cid:107) ≤ (cid:107) u n (cid:107) , from the condition R + W + P = I .Putting together the explicit and implicit parts, we obtain (cid:107) u ( i ) (cid:107) ≤ (cid:107) u n (cid:107) underthe time-step ∆ t ≤ r ∆ t FE . Remark In the case of the Broadwell model and BGK equation, the abovetheorem can be used to prove the discrete entropy decay property of the numericalmethod. Taking the following 1D BGK equation as an example, (44) ∂ t f + v∂ x f = 1 ε ( M − f ) . We set G to be the BGK operator and F be the transport operator discretized by thefirst order upwind method ( k is the spatial index): (45) ( v∂ x f ) k = v + | v | f k − f k − ∆ x + v − | v | f k +1 − f k ∆ x , together with the periodic or compactly supported boundary condition. The convexfunctional (cid:107) · (cid:107) is taken as the discrete entropy (46) S [ f ] = ∆ x (cid:88) k (cid:90) f k log f k d v. Then it can be verified that F and G satisfy the Conditions 1–3 (for more details see[11]). Therefore, the numerical solution obtained by method (31) satisfies (47) S [ f n +1 ] ≤ S [ f n ] , under the conditions listed in Theorem 1. REFERENCES[1]
U. Ascher, S. Ruuth, and R. Spiteri , Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations , Appl. Numer. Math., 25 (1997), pp. 151–167.[2]
P. Bhatnagar, E. Gross, and M. Krook , A model for collision processes in gases. I. Smallamplitude processes in charged and neutral one-component systems , Phys. Rev., 94 (1954),pp. 511–525.[3]
J. Broadwell , Shock structure in a simple discrete velocity gas , Phys. Fluids, 7 (1964),pp. 1013–1037.[4]
C. Cercignani , The Boltzmann Equation and Its Applications , Springer-Verlag, New York,1988.[5]
G.-Q. Chen, C. D. Levermore, and T.-P. Liu , Hyperbolic conservation laws with stiff relax-ation terms and entropy , Commun. Pure Appl. Math., XLVII (1994), pp. 787–830.[6]
S. Conde, S. Gottlieb, Z. Grant, and J. Shadid , Implicit and implicit-explicit strong sta-bility preserving Runge–Kutta methods with high linear order , Journal of Scientific Com-puting, 73(2) (2017), pp. 667–690.[7]
G. Dimarco and L. Pareschi , Asymptotic preserving implicit-explicit Runge-Kutta methodsfor nonlinear kinetic equations , SIAM J. Numer. Anal., 51 (2013), pp. 1064–1087.[8]
I. Higueras , Strong stability for additive Runge-Kutta methods , SIAM J. Numer. Anal., 44(2006), pp. 1735–1758.[9]
I. Higueras and T. Roldan , Positivity-preserving and entropy-decaying IMEX methods ,Monografias del Seminario Matematico Garcia de Galdeano, 33 (2006), pp. 129–136.1610]
J. Hu and R. Shu , A second-order asymptotic-preserving and positivity-preserving exponentialRunge-Kutta method for a class of stiff kinetic equations , Multiscale Model. Simul., 17(2019), pp. 1123–1146.[11]
J. Hu, R. Shu, and X. Zhang , Asymptotic-preserving and positivity-preserving implicit-explicitschemes for the stiff BGK equation , SIAM J. Numer. Anal., 56 (2018), pp. 942–973.[12]
J. Huang and C.-W. Shu , A second-order asymptotic-preserving and positivity-preservingdiscontinuous Galerkin scheme for the Kerr-Debye model , Math. Models Methods Appl.Sci., 27 (2017), pp. 549–579.[13]