Journal of Scientific Computing, Vol. 25, Nos. 1/2, November 2005 (© 2005) DOI: 10.1007/s10915-004-4636-4 Implicit–Explicit Runge–Kutta Schemes and Applications to Hyperbolic Systems with Relaxation Lorenzo Pareschi1 and Giovanni Russo2 Received October 27, 2003; accepted (in revised form) February 27, 2004 We consider new implicit–explicit (IMEX) Runge–Kutta methods for hyperbolic systems of conservation laws with stiff relaxation terms. The explicit part is treated by a strong-stability-preserving (SSP) scheme, and the implicit part is treated by an L-stable diagonally implicit Runge–Kutta method (DIRK). The schemes proposed are asymptotic preserving (AP) in the zero relaxation limit. High accuracy in space is obtained by Weighted Essentially Non Oscillatory (WENO) reconstruction. After a description of the mathematical properties of the schemes, several applications will be presented. KEY WORDS: Runge–Kutta methods; hyperbolic systems with relaxation; stiff systems; high order shock capturing schemes. AMS SUBJECT CLASSIFICATION: 65C20; 82D25. 1. INTRODUCTION Several physical phenomena of great importance for applications are described by stiff systems of differential equations in the form 1 ∂t U = F(U ) + R(U ), ε (1) where U = U (t) ∈ RN , F, R : RN → RN and ε > 0 is the stiffness parameter. System (1) may represent a system of N ODE’s or a discretization of a system of PDE’s, such as, for example, convection–diffusion equations, reaction–diffusion equations and hyperbolic systems with relaxation. 1 Department of Mathematics, University of Ferrara, Via Machiavelli 35, I-44100 Ferrara, Italy. E-mail: pareschi@dm.unife.it 2 Department of Mathematics and Computer Science, University of Catania, Via A.Doria 6, 95125 Catania, Italy. E-mail: russo@dmi.unict.it 129 0885-7474/05/1100-0129/0 © 2005 Springer Science+Business Media, Inc. 130 Pareschi and Russo In this work we consider the latter case, which in recent years has been a very active field of research, due to its great impact on applied sciences [8, 22]. For example, we mention that hyperbolic systems with relaxation appears in kinetic theory of rarefied gases, hydrodynamical models for semiconductors, viscoelasticity, multiphase flows and phase transitions, radiation hydrodynamics, traffic flows, shallow waters, etc. In one space dimension these systems have the form 1 ∂t U + ∂x F (U ) = R(U ), ε x ∈ R, (2) that corresponds to (1) for F(U ) = −∂x F (U ). In (2) the Jacobian matrix F ′ (U ) has real eigenvalues and admits a basis of eigenvectors ∀ U ∈ RN and ε is called relaxation parameter. The development of efficient numerical schemes for such systems is challenging, since in many applications the relaxation time varies from values of order one to very small values if compared to the time scale determined by the characteristic speeds of the system. In this second case the hyperbolic system with relaxation is said to be stiff and typically its solutions are well approximated by solutions of a suitable reduced set of conservation laws called equilibrium system [8]. Usually it is extremely difficult, if not impossible, to split the problem in separate regimes and to use different solvers in the stiff and non stiff regions. Thus one has to use the original relaxation system in the whole computational domain. The construction of schemes that work for all ranges of the relaxation time, using coarse grids that do not resolve the small relaxation time, has been studied mainly in the context of upwind methods using a method of lines approach combined with suitable operator splitting techniques [7, 17] and more recently in the context of central schemes [21,24]. Splitting methods have been widely used for such problems. They are attractive because of their simplicity and robustness. Strang splitting provides second order accuracy if each step is at least second order accurate [33]. This property is maintained under fairly mild assumptions even for stiff problems [15]. However, Strang splitting applied to hyperbolic systems with relaxation reduces to first order accuracy when the problem becomes stiff. The reason is that the kernel of the relaxation operator is non trivial, which corresponds to a singular matrix in the linear case, and therefore the assumptions in [15] are not satisfied. Furthermore with a splitting strategy it is difficult to obtain higher order accuracy even in non stiff regimes (high order splitting schemes can be constructed, see [9], but they are seldom used because of stability problems). Implicit–Explicit Runge–Kutta Schemes and Applications 131 Recently developed Runge–Kutta schemes overcome this difficulties, providing basically the same advantages of the splitting schemes, without the drawback of the order restriction [7, 17,36]. In this paper we will present some recent results on the development of high order, underresolved Runge–Kutta time discretization for such systems. In particular, using the formalism of implicit-explicit (IMEX) Runge-Kutta schemes [3, 4, 19,25,26,36] we derive new IMEX schemes up to order 3 that are strong-stability-preserving (SSP) for the limiting system of conservation laws (such methods were originally referred to as total variation diminishing (TVD) methods, see [10,11,32]). To this aim, we derive general conditions that guarantee the asymptotic preserving property, i.e. the consistency of the scheme with the equilibrium system and the asymptotic accuracy, i.e. the order of accuracy is maintained in the stiff limit. For a stability analysis of IMEX schemes we refer to [27]. The rest of the paper is organized as follows. In Sec. 2 we introduce the general structure of IMEX Runge–Kutta schemes. Sec. 3 is devoted to IMEX Runge–Kutta schemes applied to hyperbolic systems with relaxation. In Sec. 4 we describe space discretization obtained by conservative finite-volume and finite difference schemes. In Sec. 5 we present some numerical results concerning the accuracy of IMEX schemes when applied to a prototype hyperbolic system with relaxation. Finally in Sec. 6 we investigate the performance of the schemes in several realistic applications to shallow waters, traffic flows and granular gases. 2. IMEX RUNGE–KUTTA SCHEMES An IMEX Runge–Kutta scheme consists of applying an implicit discretization to the source terms and an explicit one to the nonstiff term. When applied to system (1) it takes the form U (i) = U n − ∆t i−1
ãij ∂x F (U (j ) ) + ∆t U n+1 = U n − ∆t ν
w̃i ∂x F (U (i) ) + ∆t j =1 i=1 ν
1 aij R(U (j ) ), ε (3) 1 wi R(U (i) ). ε (4) j =1 ν
i=1 The matrices à = (ãij ), ãij = 0 for j
i and A = (aij ) are ν × ν matrices such that the resulting scheme is explicit in F , and implicit in R. An IMEX Runge–Kutta scheme is characterized by these two matrices and the coefficient vectors w̃ = (w̃1 , . . . , w̃ν )T , w = (w1 , . . . , wν )T . 132 Pareschi and Russo Since the simplicity and efficiency of solving the algebraic equations corresponding to the implicit part of the discretization at each step is of paramount importance, it is natural to consider diagonally implicit Runge–Kutta (DIRK) schemes [14] for the source terms (aij = 0, for j > i). The IMEX Runge–Kutta schemes can be represented by a double tableau in the usual Butcher notation, c̃ Ã c A w̃ T wT , where the coefficients c̃ and c used for the treatment of non autonomous systems, are given by the usual relation c̃i = i−1
j =1 ãij , ci = i
aij . (5) j =1 The use of a DIRK scheme for R is a sufficient condition to guarantee that F is always evaluated explicitly. In the case of system (2), in order to obtain a numerical scheme, a suitable space discretization of equations (3) and (4) is required. This discretization can be performed at the level of the original system (2) in which case one has a system of ODEs that is then solved in time by the IMEX scheme (method of lines). Alternatively one can apply a suitable space discretization directly to the time discrete equations (3) and (4). Finally we remark that previously developed schemes, such as the Additive semi-implicit Runge–Kutta methods of Zhong [36], and the splitting Runge–Kutta methods of Jin et al. [17], [7] can be rewritten as IMEXRK schemes [27]. 2.1. Order Conditions The general technique to derive order conditions for Runge–Kutta schemes is based on the Taylor expansion of the exact and numerical solution. In particular, conditions for schemes of order p are obtained by imposing that the solution of system (2) at time t = t0 +
t, with a given initial condition at time t0 , agrees with the numerical solution obtained by one step of a Runge–Kutta scheme with the same initial condition, up to order ∆t p+1 . Implicit–Explicit Runge–Kutta Schemes and Applications 133 Here we report the order conditions for IMEX Runge–Kutta schemes up to order p = 3, which is already considered high order for PDE problems. We apply schemes (3) and (4) to system (2), with ε = 1. We assume that the coefficients c̃i , ci , ãij , aij satisfy conditions (5). Then the order conditions are the following First Order ν
i=1 w̃i = 1, ν
i=1 wi = 1. (6) Second Order
w̃i c̃i = 1/2,
wi ci = 1/2, (7)
w̃i ci = 1/2,
wi c̃i = 1/2. (8)
w̃i ãij c̃j = 1/6,
i i i i Third Order ij
ij  ij  wi aij cj = 1/6, w̃i ãij cj = 1/6,   ij i
i w̃i c̃i c̃i = 1/3, (9) wi ci ci = 1/3, w̃i aij c̃j = 1/6,   ij w̃i aij cj = 1/6, wi ãij cj = 1/6, (10) ij wi aij c̃j = 1/6, ij wi ãij c̃j = 1/6,    i w̃i ci ci = 1/3, i w̃i c̃i ci = 1/3, i wi c̃i c̃i = 1/3, i wi c̃i ci = 1/3. ij  Conditions (6), (7), (9) are the standard order conditions for the two tableau, each of them taken separately. Conditions (8) and (10) are new conditions that arise because of the coupling of the two schemes. The order conditions will simplify a lot if c̃ = c. For this reason only such schemes are considered in [3]. In particular, we observe that, if the two tableau differ only for the value of the matrices A, i.e. if c̃i = ci and w̃i = wi , then the standard order conditions for the two schemes are enough to ensure that the combined scheme is third order. Note, however, that this is true only for schemes up to third order. 134 Pareschi and Russo Table I. Number of Coupling Conditions in IMEX Runge–Kutta Schemes IMEX-RK order 1 2 3 4 5 6 Number of coupling conditions General case w̃i = wi c̃ = c 0 2 12 56 252 1128 0 0 3 21 110 528 0 0 2 12 54 218 c̃ = c and w̃i = wi 0 0 0 2 15 78 Higher order conditions can be derived as well using a generalization of Butcher 1-trees to 2-trees, see [12,13,19]. However the number of coupling conditions increase dramatically with the order of the schemes. The relation between coupling conditions and accuracy of the schemes is reported in Table I. 3. APPLICATIONS TO HYPERBOLIC SYSTEMS WITH RELAXATION In this section we give sufficient conditions for asymptotic preserving and asymptotic accuracy properties of IMEX schemes. This properties are strongly related to L-stability of the implicit part of the scheme. 3.1. Zero Relaxation Limit Let us consider here one-dimensional hyperbolic systems with relaxation of the form (2). The operator R : RN → RN is called a relaxation operator, and consequently (2) defines a relaxation system, if there exists a constant n × N matrix Q with rank(Q) = n < N such that QR(U ) = 0 ∀ U ∈ RN . (11) This gives n independent conserved quantities u = QU . Moreover we assume that equation R(U ) = 0 can be uniquely solved in terms of u, i.e. U = E(u) such that R(E(u)) = 0. (12) The image of E represents the manifold of local equilibria of the relaxation operator R. Implicit–Explicit Runge–Kutta Schemes and Applications 135 Using (11) in (2) we obtain a system of n conservation laws which is satisfied by every solution of (2) ∂t (QU ) + ∂x (QF (U )) = 0. (13) For vanishingly small values of the relaxation parameter ε, from (2) we get R(U ) = 0 which by (12) implies U = E(u). In this case system (2) is well approximated by the equilibrium system [8] ∂t u + ∂x G(u) = 0, (14) where G(u) = QF (E(u)). System (14) is the formal limit of system (1) as ε → 0. The solution u(x, t) of the system will be the limit of QU , with U solution of system (1), provided a suitable condition on the characteristic velocities of systems (1) and (14) is satisfied (the so called subcharacteristic condition, see [8, 35]). 3.2. Asymptotic Properties of IMEX Schemes We start with the following Definition 3.1. We say that an IMEX scheme for system (2) in the form (3) and (4) is asymptotic preserving (AP) if in the limit ε → 0 the scheme becomes a consistent discretization of the limit equation (14). Note that this definition does not imply that the scheme preserves the order of accuracy in t in the stiff limit ǫ → 0. In the latter case the scheme is said asymptotically accurate. In order to give sufficient conditions for the AP and asymptotically accurate property, we make use of the following simple Lemma 3.1. If all diagonal elements of the triangular coefficient matrix A that characterize the DIRK scheme are non zero, then lim R(U (i) ) = 0. ǫ→0 (15) Proof. In the limit ǫ → 0 from (3) we have i
j =1 aij R(U (j ) ) = 0, i = 1, . . . , ν. Since the matrix A is nonsingular, this implies R(U (i) ) = 0, i = 1, . . . , ν. 136 Pareschi and Russo In order to apply Lemma 3.1, the vectors of c and c̃ cannot be equal. In fact c̃1 = 0 whereas c1 = 0. Note that if c1 = 0 but aii = 0 for i > 1, then we still have limǫ→0 R(U (i) ) = 0 for i > 1 but limǫ→0 R(U (1) ) = 0 in general. The corresponding scheme may be inaccurate if the initial condition is not “well prepared” (R(U0 ) = 0). In this case the scheme is not able to treat the so called “initial layer” problem and degradation of accuracy in the stiff limit is expected (see Sec. 5 and references [7, 24,25].) On the other hand, if the initial condition is “well prepared” (R(U (0) ) = 0), then relation (15), i = 1, . . . , ν holds even if a11 = c1 = 0. Next we can state the following Theorem 3.1. If det A = 0, in the limit ǫ → 0, the IMEX scheme (3) and (4) applied to system (2) becomes the explicit RK scheme characterized by (Ã, w̃, c̃) applied to the limit equation (14). Proof. As in the case of the continuous system, multiplying equations (3) and (4) by the matrix Q, and making use of the relation QR(U ) = 0 ∀ U , we obtain u (i) = u0 + h i−1
ãij QF (U (j ) ) (16) w̃i QF (U (i) ), (17) j =1 and u1 = u0 + h ν
i=1 where u(i) = QU (i) , i = 1, . . . ν, u0 = QU0 , u1 = QU1 . From the previous lemma, in the stiff limit one has R(U (i) ) = 0, i = 1, . . . , ν, and, from property (12) of the relaxation operator R, this implies U (i) = E(u(i) ), i = 1, . . . , ν. Substituting this expression in (16)-(17) one obtains u(i) = u0 + h i−1
j =1 ãij G(u(j ) ) (18) Implicit–Explicit Runge–Kutta Schemes and Applications 137 and u1 = u0 + h ν
w̃i G(u(i) ), (19) i=1 where G(u) = QF (E(u)). Clearly one may claim that if the implicit part of the IMEX scheme is A-stable or L-stable the previous theorem is satisfied. Note however that this is true only if the tableau of the implicit integrator does not contain any column of zeros that makes it reducible to a simpler A-stable or L-stable form. Some remarks are in order. Remarks. (i) There is a close analogy between hyperbolic systems with stiff relaxation and differential algebraic equations (DAE) [2]. The limit system as ǫ → 0 is the analog of an index 1 DAE, in which the algebraic equation is explicitly solved in terms of the differential variable. In the context of DAE, the initial condition that we called “well prepared” is called “consistent”. (ii) This result does not guarantee the accuracy of the solution for the N − n nonconserved quantities. In fact, since the very last step in the scheme is not a projection towards the local equilibrium, a final layer effect occurs. The use of stiffly accurate schemes (i.e. schemes for which aνj = wj , j = 1, . . . , ν) in the implicit step may serve as a remedy to this problem. (iii) The theorem guarantees that in the stiff limit the numerical scheme becomes the explicit RK scheme applied to the equilibrium system, and therefore the order of accuracy of the limiting scheme for the conserved variables, is greater or equal to the order of accuracy of the original IMEX scheme. When constructing numerical schemes for conservation laws, one has to take a great care in order to avoid spurious numerical oscillations arising near discontinuities of the solution. This is avoided by a suitable choice of space discretization (see Sec. 4) and time discretization. Solution of scalar conservation equations, and equations with a dissipative source have some norm that decreases in time. It would be desirable that such property is maintained at a discrete level by the numerical method. If U n represents a vector of solution values (for example obtained from a method of lines approach in solving (14)) we recall the following [32]. 138 Pareschi and Russo Definition 3.2. A sequence {U n } is said to be strongly stable in a given norm || · || provided that ||U n+1 ||  ||U n || for all n
0. The most commonly used norms are the T V -norm and the infinity norm. A numerical scheme that maintains strong stability at discrete level is called SSP. Here we review some basic facts about RK-SSP schemes. First, it has been shown [11], under fairly general conditions, that high order SSP schemes are necessarily explicit. Second, observe that a generic explicit RK scheme can be written as U (0) = U n , i−1
(i) U = (αik U (k) + ∆tβik L(U (k) ), k=0 (ν) U n+1 = U i = 1, . . . , ν, (20) , where αik
0 and αik = 0 only if βik = 0. This representation of RK schemes (which is not unique) can be converted to a standard Butcher form i−1 in a straightforward manner. Observe that for consistency, one has k=0 αik = 1. It follows that if the scheme can be written in the form (20) with non negative coefficients βik then it is a convex combination of Forward Euler steps, with step sizes βik /αik ∆t. A consequence of this is that if Forward Euler is SSP for ∆t  ∆t ∗ , then the RK scheme is also SSP for ∆t  c∆t ∗ , with c = minik (αik /βik ) [11,29]. The constant c is a measure of the efficiency of the SSP-RK scheme, therefore for the applications it is important to have c as large as possible. For a detailed description of optimal SSP schemes and their properties see [32]. By point iii) of the above remarks it follows that if the explicit part of the IMEX scheme is SSP then, in the stiff limit, we will obtain an SSP method for the limiting conservation law. This property is essential to avoid spurious oscillations in the limit scheme for Eq. (14). In this section we give the Butcher tableau for the new second and third order IMEX schemes that satisfy the conditions of Theorem 3.1. In all these schemes the implicit tableau corresponds to an L-stable scheme, that is w T A−1 e = 1, e being a vector whose components are all equal to 1 [14], whereas the explicit tableau is SSPk, where k denotes the order of the SSP scheme. Several examples of asymptotically SSP schemes are reported in Tables (II)–(VI). We shall use the notation SSPk(s, σ, p), where the triplet (s, σ, p) characterizes the number s of stages of the implicit Implicit–Explicit Runge–Kutta Schemes and Applications Tableau for the Explicit (left) Implicit (right) IMEXSSP2(2,2,2) L-Stable Scheme Table II. 0 0 1 1 139 0 0 γ γ 0 1 − γ 1 − 2γ γ 1/2 1/2 1/2 γ = 1 − √1 2 1/2 Table III. Tableau for the Explicit (left) Implicit (right) IMEX-SSP2(3,2,2) Stiffly Accurate Scheme 0 0 0 0 0 0 1 0 1 1/2 1/2 0 0 0 −1/2 1/2 0 0 1/2 1/2 1 0 0 0 0 1/2 1/2 0 1/2 1/2 Table IV. Tableau for the Explicit (left) Implicit (right) IMEX-SSP2(3,3,2) Stiffly Accurate Scheme 0 0 0 0 1/2 1/2 0 0 1 1/2 1/2 0 1/4 1/4 0 0 1/4 0 1/4 0 1 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 scheme, the number σ of stages of the explicit scheme and the order p of the IMEX scheme. 4. IMEX-WENO SCHEMES For simplicity we consider the case of the single scalar equation 1 ut + f (u)x = r(u). ε (21) 140 Pareschi and Russo Table V. Tableau for the Explicit (left) Implicit (right) IMEXSSP3(3,3,2) L-Stable scheme 0 0 0 0 1 1 0 0 1/2 1/4 1/4 0 γ γ 0 1 − γ 1 − 2γ γ 1/2 1/2 − γ 0 1/6 1/6 2/3 Table VI. 0 0 1 1/2 0 0 0 0 1/6 0 0 γ γ = 1 − √1 2 1/6 2/3 Tableau for the Explicit (left) Implicit (right) IMEXSSP3(4,3,3) L-Stable scheme 0 0 1 1/4 0 0 0 1/4 0 0 0 0 α 0 1 1/2 α 0 0 −α α 0 0 1−α α β η 1/2 − β − η − α 0 1/6 1/6 2/3 0 α = 0.24169426078821, 0.12915286960590 β 1/6 = 1/6 0 0 0 α 2/3 0.06042356519705 η = We have to distinguish between schemes based on cell averages (finite volume approach, widely used for conservation laws) and schemes based on point values (finite difference approach). Let ∆x and
t be the mesh widths. We introduce the grid points xj = j ∆x, 1 xj +1/2 = xj + ∆x, 2 j = . . . , −2, −1, 0, 1, 2, . . . and use the standard notations unj = u(xj , t n ), ūnj = 1 ∆x  xj +1/2 u(x, t n ) dx. xj −1/2 4.1. Finite Volumes Integrating equation (21) on Ij = [xj −1/2 , xj +1/2 ] and dividing by
x we obtain d ūj 1 1 =− [f (u(xj +1/2 , t)) − f (u(xj −1/2 , t))] + r(u)j . dt ∆x ε (22) In order to convert this expression into a numerical scheme, one has to approximate on right hand side with a function of the cell averages {ū(t)}j , which are the basic unknowns of the problem. Implicit–Explicit Runge–Kutta Schemes and Applications 141 The first step is to perform a reconstruction of the unknown function u(x, t) by a piecewise polynomial, starting from cell averages ūj (t). Given {ūj }, compute a piecewise polynomial reconstruction
R(x) = Pj (x)χj (x), j where Pj (x) is a polynomial satisfying some accuracy and non oscillatory property, and χj (x) is the indicator function of cell Ij . For first order schemes, the reconstruction is piecewise constant, while second order schemes can be obtained by a piecewise linear reconstruction. Higher order schemes are obtained by higher order polynomials. The reconstruction step is very delicate, because the function u(x, t) may be discontinuous. To this goal, suitable techniques, such as Weighted Essentially NonOscillatory (WENO), have been developed. The reader can consult the review by Shu [30] and references therein for a more detailed description of high order non oscillatory reconstructions. The flux function at the edge of the cell can be computed by using a suitable numerical flux function, consistent with the analytical flux, + f (u(xj +1/2 , t)) ≈ F (u− j +1/2 , uj +1/2 ), where the quantities u± j +1/2 are obtained from the reconstruction, as ± uj +1/2 = limx→x ± R(x). j +1/2 Example of numerical flux functions are the Godunov flux + + ∗ − F (u− j +1/2 , uj +1/2 ) = f (u (uj +1/2 , uj +1/2 )), + where u∗ (u− j +1/2 , uj +1/2 ) is the solution of the Riemann problem at xj +1/2 , + corresponding to the states u− j +1/2 and uj +1/2 , and the Local Lax Friedrichs flux (also known as Rusanov flux), + − + + − 1 F (u− j +1/2 , uj +1/2 ) = 2 [(f (uj +1/2 ) + f (uj +1/2 )) − α(uj +1/2 − uj +1/2 )], where α = maxw |f ′ (w)|, and the maximum is taken over the relevant range of w. The two examples constitute two extreme cases of numerical fluxes: the Godunov flux is the most accurate and the one that produces the best results for a given grid size, but it is very expensive, since it requires the solution to the Riemann problem. Local Lax-Friedrichs flux, on the other hand, is less accurate, but much cheaper. This latter is the numerical flux that has been used throughout all the calculations performed for 142 Pareschi and Russo this paper. The difference in resolution provided by the various numerical fluxes becomes less relevant with the increase in the order of accuracy of the method. The right hand side of Eq.(22) contains the average of the source term r(u) instead of the source term evaluated at the average of u, r(ū). The two quantities agree within second order accuracy r(u)j = r(ūj ) + O(∆x 2 ). This approximation can be used to construct schemes up to second order. First order (in space) semidiscrete schemes can be obtained using the numerical flux function F (ūj , ūj +1 ) in place of f (u(xj +1/2 , t)), F (ūj , ūj +1 ) − F (ūj −1 , ūj ) 1 d ūj + r(ūj ). =− ∆x dt ε (23) Second order schemes are obtained by using a piecewise linear reconstruction in each cell, and evaluating the numerical flux on the two sides of the interface − + − + F (uj +1/2 , uj +1/2 ) − F (uj −1/2 , uj −1/2 ) 1 d ū =− + r(ūj ). dt ∆x ε The quantities at cell edges are computed by piecewise linear reconstruction. For example, u− j +1/2 = ūj + ∆x ′ u , 2 j where the slope u′j is a first order approximation of the space derivative of u(x, t), and can be computed by suitable slope limiters (see, for example, [20] for a discussion on TVD slope limiters.) For schemes of order higher than second, a suitable quadrature formula is required to approximate g(u)j . For example, for third and fourth order schemes, one can use Simpson’s rule   − r(u)j ≈ 16 r(u+ ) + 4r(u ) + r(u ) j j −1/2 j +1/2 , − where the pointwise values u+ j −1/2 , uj , uj +1/2 are obtained from the reconstruction. For a general problem, this has the effect that the source term couples the cell averages of different cells, thus making almost impractical the use of finite volume methods for high order schemes applied to stiff sources. Implicit–Explicit Runge–Kutta Schemes and Applications 143 Note, however, that in many relevant cases of hyperbolic systems with relaxation the implicit step, thanks to the conservation properties of the system, can be explicitly solved, and finite volume methods can be successfully used. We mention here all relaxation approximation of Jin– Xin type [18], some simple discrete velocity models, such as Carlemann and Broadwell models, monoatomic gas in Extended Thermodynamics, semiconductor models, and shallow water equations. 4.2. Finite Differences In a finite difference scheme the basic unknown is the pointwise value of the function, rather than its cell average. Osher and Shu observed that it is possible to write a finite difference scheme in conservative form [31]. Let us consider the equation ∂u ∂f 1 + = r(u) ∂t ∂x ε and write ∂f fˆ(u(x + (∆x/2), t)) − fˆ(u(x − (∆x/2), t)) (u(x, t)) = . ∂x ∆x The relation between f and fˆ is the following. Let us consider the sliding average operator 1 ū(x, t) = u x ≡ ∆x  x+(∆x/2) u(ξ, t) dξ . x−(∆x/2) Differentiating with respect to x one has      ∂ ū 1 ∆x ∆x = ,t −u x − ,t . u x+ ∂x ∆x 2 2 Therefore the relation between f and fˆ is the same that exists between ū(x, t) and u(x, t), namely, flux function f is the cell average of the function fˆ. This also suggests a way to compute the flux function. The technique that is used to compute pointwise values of u(x, t) at the edge of the cell from cell averages of u can be used to compute fˆ(u(xj +1/2 , t)) from f (u(xj , t)). This means that in the finite difference method it is the flux function which is computed at xj and then reconstructed at xj +1/2 . But the reconstruction at xj +1/2 may be discontinuous. Which value should 144 Pareschi and Russo one use? A general answer to this question can be given if one considers flux functions that can be split f (u) = f + (u) + f − (u) (24) with the condition that df + (u)
0, du df − (u)  0. du (25) There is a close analogy between flux splitting and numerical flux functions. In fact, if a flux can be split as (24), then F (a, b) = f + (a) + f − (b) will define a monotone consistent flux, provided condition (25) is satisfied. Together with non oscillatory reconstructions (such as WENO) and SSP time discretization, the monotonicity condition will ensure that the overall scheme will not produce spurious numerical oscillations (see, for example, [20,30].) This is the case, for example, of the local Lax–Friedrichs flux. A finite difference scheme therefore takes the following form duj 1 1 =− [F̂j +1/2 − F̂j −1/2 ] + g(uj ), dt ∆x ε ˆ− + F̂j +1/2 = fˆ+ (u− j +1/2 ) + f (uj +1/2 ), fˆ+ (u− j +1/2 ) is obtained by • • computing f + (ul ) and interpret it as cell average of fˆ+ ; performing pointwise reconstruction of fˆ+ in cell j , and evaluate it in xj +1/2 . fˆ− (u+ j +1/2 ) is obtained by • • computing f − (ul ), interpret as cell average of fˆ− ; performing pointwise reconstruction of fˆ− in cell j + 1, and evaluate it in xj +1/2 . A detailed account on high order finite difference schemes can be found in [30]. At variance with finite volume schemes, since the source is evaluated pointwise, finite difference schemes do not couple the cells. This is a big Implicit–Explicit Runge–Kutta Schemes and Applications 145 advantage in the cases in which the implicit step can not be explicitly solved. Remarks. (i) Finite difference can be used only with uniform (or smoothly varying) mesh. In this respect finite volume are more flexible, since they can be used even with unstructured grids. (ii) The treatment presented for the scalar equation can be extended to systems, with a minor change of notation. In particular, the parameter α appearing in the local Lax-Friedrichs flux will be computed using the spectral radius of the Jacobian matrix. (iii) For schemes applied to systems, better results are usually obtained if one uses characteristic variables rather than conservative variables in the reconstruction step. The use of conservative variables may result in the appearance of small spurious oscillations in the numerical solution. For a treatment of this effect see, for example, [28]. 5. NUMERICAL TESTS In this section we investigate numerically the convergence rate and the zero relaxation limit behavior of the schemes. To this aim we apply the IMEX-WENO schemes to the Broadwell equations of rarefied gas dynamics [7, 17,21,24]. In all the computations presented in this paper we used finite difference WENO schemes with local Lax-Friedrichs flux and conservative variables. Of course the sharpness of the resolution of the numerical results can be improved using a less dissipative flux. As a comparison, together with the new IMEX-SSP schemes, we have considered the second order ARS(2,2,2) method presented in [3], which was the one recommended by the authors. For this scheme, since c1 = 0, there will be a degradation of accuracy due to the initial layer. In all figures, the value of N represents the number of grid points in space. The kinetic model is characterized by a hyperbolic system with relaxation of the form (2) with U = (ρ, m, z), F (U ) = (m, z, m),   R(U ) = 0, 0, 21 (ρ 2 + m2 − 2ρz) . Here ε represents the mean free path of particles. The only conserved quantities are the density ρ and the momentum m. 146 Pareschi and Russo In the fluid-dynamic limit ε → 0 we have z = zE ≡ ρ 2 + m2 , 2ρ (26) and the Broadwell system is well approximated by the reduced system (14) with   m u = (ρ, ρv), G(u) = ρv, 21 (ρ + ρv 2 ) , v= , ρ which represents the corresponding Euler equations of fluid-dynamics. Convergence Rates We have considered a periodic smooth solution with initial data as in [7, 21] given by ρ(x, 0) = 1 + aρ sin 2πx L , 2 2 v(x, 0) = 21 + av sin 2πLx , +m(x,0) z(x, 0) = az ρ(x,0) . 2ρ(x,0) (27) In our computations we used the parameters aρ = 0.3, av = 0.1, az = 0.2 (initial layer), az = 1.0 (no initial layer) and L = 20, (28) and we integrate the equations for t ∈ [0, 30]. A Courant number
t/
x = 0.6 has been used. The plots of the relative error are given in Fig. 1. Notice how, in absence of an initial layer, all schemes tested have the prescribed order of accuracy both in the non stiff and in the stiff limit, with some degradation of the accuracy at intermediate regimes. As expected, scheme ARS(2,2,2), shows a degradation of the accuracy when an initial layer is present. Next we test the shock capturing properties of the schemes in the case of non smooth solutions characterized by the following two Riemann problems [7] ρl = 2, ρr = 1, ρl = 1, ml = 1, mr = 0.13962, ρr = 0.2, ml = 0, mr = 0, zl = 1, zr = 1, zl = 1, zr = 1, x < 0.2, x > 0.2, x < 0, x > 0. (29) (30) Implicit–Explicit Runge–Kutta Schemes and Applications Relative error in density vs N, ε = 1, no initial layer –2 10 10 –3 10 10 –4 10 Errρ Errρ 10 –5 10 –6 10 SSP2–222 SSP2–322 SSP2–332 SSP3–332 SSP3–433 ARS–222 –7 10 –8 10 10 Errρ 10 10 10 10 –2 10 10 10 3 10 –4 –5 10 SSP2–222 SSP2–322 SSP2–332 SSP3–332 SSP3–433 ARS–222 –6 1 10 2 10 10 –4 10 1 10 10 N 2 10 3 –8 1 10 2 10 3 N Relative error in density vs N, ε = 10–3, with initial layer SSP2–222 SSP2–322 SSP2–332 SSP3–332 SSP3–433 ARS–222 –6 10 10 SSP2–222 SSP2–322 SSP2–332 SSP3–332 SSP3–433 ARS–222 –7 –5 10 –7 –6 –7 10 –3 –6 –5 –4 10 SSP2–222 SSP2–322 SSP2–332 SSP3–332 SSP3–433 ARS–222 –4 –3 10 3 –5 –3 –2 N Relative error in density vs N, –6 ε = 10 , no initial layer –2 10 10 Relative error in density vs N, ε = 1, with initial layer –2 10 Errρ Errρ 10 2 10 10 10 10 –3 –7 10 10 N Relative error in density vs N, ε = 10–3, no initial layer 10 10 10 Errρ 10 1 10 147 1 2 10 10 3 N Relative error in density vs N, –6 ε = 10 , with initial layer –2 –3 –4 –5 SSP2–222 SSP2–322 SSP2–332 SSP3–332 SSP3–433 ARS–222 –6 –7 10 10 1 10 2 10 3 N Fig. 1. Relative errors for density ρ in the Broadwell equations with initial data (28). Left column az = 1.0 (no initial layer), right column az = 0.2 (initial layer). From top to bottom, ε = 1.0, 10−3 , 10−6 . 148 Pareschi and Russo For brevity we report the numerical results obtained with the second order IMEX-SSP2(2,2,2) and third order IMEX-SSP3(4,3,3) schemes that we will refer to as IMEX-SSP2-WENO and IMEX-SSP3-WENO respectively. The result are shown in Figs. 2 and 3 for a Courant number
t/
x = 0.5. Both schemes, as expected, give an accurate description of the solution in all different regimes even using coarse meshes that do not resolve the small scales. In particular the shock formation in the fluid limit is well captured without spurious oscillations. We refer to [1, 7, 17,21,24] for a comparison of the present results with previous ones. 6. APPLICATIONS Finally we present some numerical results obtained with IMEX-SSP2WENO and IMEX-SSP3-WENO concerning situations in which hyperbolic systems with relaxation play a major role in applications. The results have been obtained with N = 200 grid points. As usual the reference solution is computed on a much finer grid. 6.1. Shallow Water First we consider a simple model of shallow water flow [17] ∂t h + ∂x (hv) = 0,     1 h h ∂t (hv) + ∂x h + h2 = −v , 2 ǫ 2 (31) where h is the water height with respect to the bottom and hv the flux. The zero relaxation limit of this model is given by the inviscid Burgers equation. The initial data we have considered is [17] h = 1 + 0.2 sin(8π x), hv = h2 , 2 (32) with x ∈ [0, 1]. The solution at t = 0.5 in the stiff regime ǫ = 10−8 using periodic boundary conditions is given in Fig. 4. For IMEX-SSP2-WENO the dissipative effect due to the use of the Lax-Friedrichs flux is very pronounced. As expected, this effect becomes less relevant with the increase of the order of accuracy. We refer to [17] for a comparison with the present results. Implicit–Explicit Runge–Kutta Schemes and Applications 2.5 IMEX–SSP2–WENO, ε=1, t=0.5, N=200 2.5 1.5 1 0.5 0 –1 –0.8–0.6–0.4–0.2 0 0.2 0.4 0.6 0.8 1 x 2.5 IMEX–SSP2–WENO, ε=0.02, t=0.5, N=200 1 0.5 2.5 1.5 1 0.5 0 –1 –0.8–0.6–0.4–0.2 0 0.2 0.4 0.6 0.8 1 x 1.5 1 0.5 0 –1 –0.8–0.6–0.4–0.2 0 0.2 0.4 0.6 0.8 1 x –8 –8 IMEX–SSP2–WENO, ε=10 , t=0.5, N=200 2.5 1 0.5 0 –1 –0.8–0.6–0.4–0.2 0 0.2 0.4 0.6 0.8 1 x ρ(x,t), m(x,t), z(x,t) 1.5 IMEX–SSP3–WENO, ε=10 , t=0.5, N=200 2 2 ρ(x,t), m(x,t), z(x,t) IMEX–SSP3–WENO, ε=0.02, t=0.5, N=200 2 ρ(x,t), m(x,t), z(x,t) ρ(x,t), m(x,t), z(x,t) 1.5 0 –1 –0.8–0.6–0.4–0.2 0 0.2 0.4 0.6 0.8 1 x 2 2.5 IMEX–SSP3–WENO, ε=1, t=0.5, N=200 2 ρ(x,t), m(x,t), z(x,t) ρ(x,t), m(x,t), z(x,t) 2 149 1.5 1 0.5 0 –1 –0.8–0.6–0.4–0.2 0 0.2 0.4 0.6 0.8 1 x Fig. 2. Numerical solution of the Broadwell equations with initial data (29) for ρ(◦), m(∗) and z(+) at time t = 0.5. Left column IMEX-SSP2-WENO scheme, right column IMEXSSP3-WENO scheme. From top to bottom, ε = 1.0, 0.02, 10−8 . 6.2. Traffic Flows In [5] a new macroscopic model of vehicular traffic has been presented. The model consists of a continuity equation for the density ρ of vehicles together with an additional velocity equation that describes the 150 Pareschi and Russo –8 –8 1.2 IMEX–SSP2–WENO, ε=10 , t=0.25, N=200 1.2 1 ρ(x,t), m(x,t), z(x,t) 1 ρ(x,t), m(x,t), z(x,t) IMEX–SSP3–WENO, ε=10 , t=0.25, N=200 0.8 0.6 0.4 0.2 0.8 0.6 0.4 0.2 0 0 –0.2 –1 –0.8–0.6–0.4–0.2 0 0.2 0.4 0.6 0.8 1 x –0.2 –1 –0.8–0.6–0.4–0.2 0 0.2 0.4 0.6 0.8 1 x Fig. 3. Numerical solution of the Broadwell equations with initial data (30) for ρ(◦), m(∗) and z(+) at time t = 0.25 for ε = 10−8 . Left IMEX-SSP2-WENO scheme, right IMEX-SSP3WENO scheme. –8 IMEX–SSP2–WENO, ε=1e–8, t=0.5, N=200 1.3 1.2 1.2 1.1 1.1 1 h(x,t), hv(x,t) h(x,t), hv(x,t) 1.3 0.9 0.8 0.7 0.6 1 0.9 0.8 0.7 0.6 0.5 0.5 0.4 0.4 0.3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x 1 IMEX–SSP3–WENO, ε=10 , t=0.5, N=200 0.3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x 1 Fig. 4. Numerical solution of the shallow water model with initial data (32) for h(◦) and hv(∗) at time t = 0.5 for ε = 10−8 . Left IMEX-SSP2-WENO scheme, right IMEX-SSP3-WENO scheme. mass flux variations due to the road conditions in front of the driver. The model can be written in conservative form as follows ∂t ρ + ∂x (ρv) = 0, (33) ρ ∂t (ρw) + ∂x (vρw) = A (V (ρ) − v), T where w = v + P (ρ) with P (ρ) a given function describing the anticipation of road conditions in front of the drivers and V (ρ) describes the dependence of the velocity with respect to the density for an equilibrium situation. The parameter T is the relaxation time and A > 0 is a positive constant. Implicit–Explicit Runge–Kutta Schemes and Applications 151 If the relaxation time goes to zero, under the subcharacteristic condition −P ′ (ρ)  V ′ (ρ)  0, ρ > 0, we obtain the Lighthill–Whitham [35] model ∂t ρ + ∂x (ρV (ρ)) = 0. (34) A typical choice for the function P (ρ) is given by ⎧  γ ⎨ cv ρ γ > 0, γ ρm  P (ρ) = ⎩ cv ln ρ γ = 0, ρm where ρm is a given maximal density and cv a constant with dimension of velocity. In our numerical results we assume A = 1 and an equilibrium velocity V (ρ) fitting to experimental data [6] V (ρ) = vm π/2 + arc tan (α(ρ/ρm − β)/(ρ/ρm − 1)) π/2 + arc tan (αβ) with α = 11, β = 0.22 and vm a maximal speed. We consider γ = 0 and, in order to fulfill the subcharacteristic condition, assume cv = 2. All quantities are normalized so that vm = 1 and ρm = 1. We consider a Riemann problem centered at x = 0 with left and right states ρL = 0.05, vL = 0.05, ρR = 0.05, vR = 0.5. (35) The solution at t = 1 for T = 0.2 is given in Fig. 5. The figure shows the development of the density of the vehicles. Both schemes gives very similar results. Again, in the second order scheme the shock is smeared out if compared to the third order case. See [6] for more numerical results. 6.3. Granular Gases We consider the continuum equations of Euler type for a granular gas [16,34]. These equations have ben derived for a dense gas composed of inelastic hard spheres. The model reads ρt + (ρu)x = 0, (ρu)t + (ρu2 + p)x = ρg, (36)     (1 − e2 ) 2 3 2 3/2 1 3 1 3 G(ρ)ρ T , 2 ρu + 2 ρT t + 2 ρu + 2 uρT + pu x = − ǫ 152 Pareschi and Russo IMEX–SSP2–WENO, ε=0.2, t=1, N=200 0.052 0.05 0.05 0.048 0.048 ρ(x,t),ρ v(x,t) ρ(x,t),ρ v(x,t) 0.052 0.046 0.044 0.046 0.044 0.042 0.042 0.04 0.04 0.038 –2.5 –2 –1.5 –1 –0.5 x 0 0.5 1 1.5 IMEX–SSP2–WENO, ε=0.2, t=1, N=200 0.038 –2.5 –2 –1.5 –1 –0.5 x 0 0.5 1 1.5 Fig. 5. Numerical solution of the traffic model with initial data (35) for ρ(◦) and ρv(∗) at time t = 1 for ε = 0.2. Left IMEX-SSP2-WENO scheme, right IMEX-SSP3-WENO scheme. where e is the coefficient of restitution, g the acceleration due to gravity, ǫ a relaxation time, p is the pressure given by p = ρT (1 + 2(1 + e)G(ρ)) and G(ρ) is the statistical correlation function. In our experiments we assume  ν G(ρ) = ν 1 − νM 4/3νM −1 , where ν = σ 3 ρπ/6 is the volume fraction, σ is the diameter of a particle, and νM = 0.64994 is 3D random close-packed constant. We consider the following initial data [23] on the interval [0, 10] ρ = 34.37746770, v = 18, P = 1589.2685472, (37) which corresponds to a supersonic flow at Mach number Ma = 7 (the ratio of the mean fluid speed to the speed of sound). Zero-flux boundary condition have been used on the bottom (right) boundary whereas on the top (left) we have an ingoing flow characterized by (37). The values of the restitution coefficient and the particle diameter have been taken e = 0.97 and σ = 0.1. We report the solution at t = 0.2 with ǫ = 0.01 in Fig. 6 (see [23] for similar results). Due to the nonlinearity of the source term, the implicit solver has been efficiently solved using Newton’s method in each cell thanks to the use of finite difference space discretization. Implicit–Explicit Runge–Kutta Schemes and Applications IMEX–SSP2–WENO, ε=0.01, t=0.2, N=200 0.5 0.4 0.4 0.35 0.35 0.3 ν(x,t) ν(x,t) 0.3 0.25 0.25 0.2 0.2 0.15 0.15 0.1 0.1 0.05 0.05 0 2 3 4 250 200 200 150 150 ρ u(x,t) ρ u(x,t) 1 0 5 6 7 8 9 10 x IMEX–SSP2–WENO, ε=0.01, t=0.2, N=200 250 100 50 0 0 1 2 3 4 5 6 7 8 9 10 x IMEX–SSP3–WENO, ε=0.01, t=0.2, N=200 100 50 0 1 2 3 4 6 5 x 6 7 8 9 0 10 0 1 2 3 4 5 6 7 8 9 10 x 6 x 10 IMEX–SSP3–WENO, ε=0.01, t=0.2, N=200 2.5 x 10 IMEX–SSP2–WENO, ε=0.01, t=0.2, N=200 2.5 2 2 1.5 1.5 p(x,t) p(x,t) IMEX–SSP3–WENO, ε=0.01, t=0.2, N=200 0.5 0.45 0.45 0 153 1 1 0.5 0.5 0 0 –0.5 0 1 2 3 4 5 x 6 7 8 9 10 –0.5 0 1 2 3 4 5 x 6 7 8 9 10 Fig. 6. Numerical solution of the hydrodynamical model of a granular gas with initial data (37). Left column IMEX-SSP2-WENO scheme, right column IMEX-SSP3-WENO scheme. From top to bottom, volume fraction ν, momentum density ρu, and pressure p. Both methods provide a good description of the shock that propagates backward after the particles impact with the bottom. Note that the second order method provides excessive smearing of the layer at the right boundary. This problem is not present in the third order scheme. However, due to the use of conservative variables, we can observe the presence of small spurious oscillations in the pressure profile. 154 Pareschi and Russo ACKNOWLEDGMENTS The authors would like to thank the unknown referees for their useful remarks and for pointing out reference [15]. This work was supported by the European network HYKE, funded by the EC as contract HPRN-CT2002-00282. REFERENCES 1. Arora, M., and Roe, P. L. (1998). Issues and strategies for hyperbolic problems with stiff source terms in Barriers and challenges in computational fluid dynamics, Hampton, VA, 1996, Kluwer Academic Publication, Dordrecht, pp. 139–154. 2. Ascher, U., and Petzold, L. (1998). Computer Methods for Ordinary Differential Equations, and Differential Algebraic Equations, SIAM, Philadelphia. 3. Ascher, U., Ruuth, S., and Spiteri, R. J. (1997). Implicit-explicit Runge–Kutta methods for time dependent partial differential equations. Appl. Numer. Math. 25, 151–167. 4. Ascher, U., Ruuth, S., and Wetton, B. (1995). Implicit-explicit methods for time dependent PDE’s, SIAM J. Numer. Anal. 32, 797–823. 5. Aw, A., and Rascle, M. (2000). Resurrection of second order models of traffic flow? SIAM. J. Appl. Math. 60, 916–938. 6. Aw, A., Klar, A., Materne, T., and Rascle, M. (2002). Derivation of continuum traffic flow models from microscopic follow-the-leader models, SIAM J. Appl. Math. 63, 259–278. 7. Caflisch, R. E., Jin, S., and Russo, G. (1997). Uniformly accurate schemes for hyperbolic systems with relaxation. SIAM J. Numer. Anal. 34, 246–281. 8. Chen, G. Q., Levermore, D., and Liu, T. P. (1994). Hyperbolic conservations laws with stiff relaxation terms and entropy. Comm. Pure Appl. Math. 47, 787–830. 9. Dia, B. O., and Schatzman, M. (1996). Commutateur de certains semi-groupes holomorphes et applications aux directions alternées. Math. Modelling Num. Anal. 30, 343–383. 10. Gottlieb, S., and Shu, C. -W. (1998). Total variation diminishing Runge–Kutta schemes. Math. Comp. 67, 73–85. 11. Gottlieb, S., Shu, C. -W., and Tadmor, E. (2001). Strong-stability-preserving high order time discretization methods. SIAM Review 43, 89–112. 12. Hairer, E. (1981). Order conditions for numerical methods for partitioned ordinary differential equations. Numerische Mathematik 36, 431–445. 13. Hairer, E., Nørsett, S. P., and Wanner, G. (1987). Solving Ordinary Differential Equations, Vol.1 Nonstiff problems, Springer-Verlag, New York. 14. Hairer, E., and Wanner, G. (1987). Solving Ordinary Differential Equations, Vol.2 Stiff and Differential-algebraic Problems, Springer-Verlag, New York. 15. Jahnke, T., and Lubich, C. (2000). Error bounds for exponential operator splitting, BIT. 735–744. 16. Jenkins, J., and Richman, M. (1985). Grad’s 13-moment system for a dense gas of inelastic spheres. Arch. Rat. Mech. Anal. 87, 355–377. 17. Jin, S. (1995). Runge–Kutta methods for hyperbolic systems with stiff relaxation terms. J. Comput. Phys. 122, 51–67. 18. Jin, S., and Xin, Z. P. (1995). The relaxation schemes for systems of conservation laws in arbitrary space dimensions. Comm. Pure Appl. Math. 48(3), 235–276. Implicit–Explicit Runge–Kutta Schemes and Applications 155 19. Kennedy, C. A., and Carpenter, M. H. (2003). Additive Runge–Kutta schemes for convection-diffusion-reaction equations. Appl. Numer. Math. 44, 139–181. 20. LeVeque, R. J. (1992). Numerical Methods for Conservation Laws, Lectures in Mathematics, Birkhauser Verlag, Basel. 21. Liotta, S. F., Romano, V., and Russo, G. (2000). Central schemes for balance laws of relaxation type. SIAM J. Numer. Anal. 38, 1337–1356. 22. Liu, T. P. (1987). Hyperbolic conservation laws with relaxation. Comm. Math. Phys. 108, 153–175. 23. Marquina, A., and Serna, S. (2004). Capturing Shock Waves in Inelastic Granular Gases, UCLA-CAM Report, 04–04. 24. Pareschi, L. (2001). Central differencing based numerical schemes for hyperbolic conservation laws with stiff relaxation terms. SIAM J. Num. Anal. 39, 1395–1417. 25. Pareschi, L., and Russo, G. (2001). Implicit-explicit Runge-Kutta schemes for stiff systems of differential equations. Adv. Theory Comput Math. 3, 269–289. 26. Pareschi, L., and Russo, G. (2003). High order asymptotically strong-stability-preserving methods for hyperbolic systems with stiff relaxation, Proceedings HYP2002, Pasadena USA, Springer, Berlin, 241–255. 27. Pareschi, L., and Russo, G. (2004). Stability analysis of implicit–explicit Runge–Kutta schemes. preprint. 28. Qiu, J. and Shu, C.-W. (2002). On the construction, comparison, and local characteristic decomposition for high-order central WENO schemes. J. Comput. Phys. 183(1), 187–209. 29. Shu, C.-W. (1988). Total variation diminishing time discretizations, SIAM J. Sci. Stat. Comput. 9, 1073–1084. 30. Shu, C.-W. (2000). Essentially Non Oscillatory and Weighted Essentially Non Oscillatory Schemes for Hyperbolic Conservation Laws, in Advanced numerical approximation of nonlinear hyperbolic equations, Lecture Notes in Mathematics, 1697, 325–432. 31. Shu, C.-W. and Osher, S. (1988). Efficient implementation of essentially nonoscillatory shock-capturing schemes. J. Comput. Phys. 77(2), 439–471. 32. Spiteri, R. J., and Ruuth, S. J. (2002). A new class of optimal strong-stability-preserving time discretization methods. SIAM. J. Num. Anal. 40(2), 469–491. 33. Strang, G. (1968). On the construction and comparison of difference schemes. SIAM J. Numer. Anal. 5, 505–517. 34. Toscani, G., Kinetic and hydrodinamic models of nearly elastic granular flows. Monatsch. Math. 142 (1–2), 179–192. 35. Whitham, G. B. (1974). Linear and Nonlinear waves, Wiley, New York. 36. Zhong, X. (1996). Additive semi-implicit Runge–Kutta methods for computing high speed nonequilibrium reactive flows. J. Comp. Phys. 128, 19–31.