A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation S. Boscarino∗ L. Pareschi† G.Russo‡ January 17, 2017 7 1 0 2 Abstract n a In this paper we consider the development of Implicit-Explicit (IMEX) Runge-Kutta J schemes for hyperbolic systems with multiscale relaxation. In such systems the scaling 6 depends on an additional parameter which modifies the nature of the asymptotic behavior 1 which can be either hyperbolic or parabolic. Because of the multiple scalings, standard IMEX Runge-Kutta methods for hyperbolic systems with relaxation loose their efficiency ] A and a different approach should be adopted to guarantee asymptotic preservation in stiff regimes. We show that the proposed approach is capable to capture the correct asymptotic N limit of the system independently of the scaling used. Several numerical examples confirm . h our theoretical analysis. t a m Key words. IMEXRunge-Kutta methods, hyperbolicconservationlawswith sources, diffusion equations, hydrodynamic limits, stiff systems, asymptotic-preserving schemes. [ 1 AMS subject classification. 65C20, 65M06, 76D05, 82C40. v 0 7 Contents 3 4 0 1 Introduction 2 . 1 0 2 Previous IMEX Runge-Kutta approaches 4 7 2.1 Additive approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1 2.2 Partitioned approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 : v 2.3 Hybrid additive-partitioned approach . . . . . . . . . . . . . . . . . . . . . . . . . 5 i X r 3 A unified IMEX Runge-Kutta approach 6 a 3.1 Description of the method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.2 Extension to general IMEX Runge-Kutta schemes . . . . . . . . . . . . . . . . . 8 3.2.1 Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2.2 The unified IMEX-RK setting . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.3 Removing the parabolic stiffness . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4 Numerical applications 15 ∗Mathematics and Computer Science Department, University of Catania, Italy ([email protected]). †Mathematics Department, University of Ferrara, Italy ([email protected]). ‡Mathematics and Computer Science Department, University of Catania, Italy ([email protected]). 1 5 Conclusions 23 A Appendix: IMEX-RK schemes 23 1 Introduction Hyperbolic systems with relaxation often contain multiple space-time scales which may differ of several orders of magnitude due to the various physical parameters characterizing the model. This is the case, for example, of kinetic equations close to the hydrodynamic limit [4, 14, 12]. In such regimes these systems can be more conveniently described in terms of fluid-dynamic equations, when they are considered on a suitable space-time scale. As a prototype system, that we will use to illustrate the subsequent theory, we consider the following  ∂ u+∂ v = 0,  t x  (1) 1 1  ∂ v+ ∂ p(u) = − (v−f(u)), α ∈ [0,1] t ε2α x ε1+α (cid:112) wherep(cid:48)(u) > 0. System(1)ishyperbolicwithtwodistinctrealcharacteristicsspeeds± p(cid:48)(u)/εα. Note that the scaling introduced in (1) corresponds to the study of the limiting behavior of the solution for the usual hyperbolic system with a singular perturbation source  ∂ u+∂ V = 0,  τ ξ  (2) 1  ∂ V +∂ p(u) = − (V −F(u)), τ ξ ε under the rescaling t = εατ, ξ = x, v(x,t) = V(ξ,τ)/εα and f(u) = F(u)/εα. For α = 0, system reduces to the usual hyperbolic scaling (2). However, if α > 0, we are looking at larger microscopic times. In particular, for small values of ε, using the Chapman-Enskog expansion, the behavior of the solution to (1) is, at least formally, governed by the following nonlinear parabolic equation  v = f(u)−ε1−α∂ p(u)+ε1+αf(cid:48)(u)2∂ u+O(ε2),  x x    (cid:34)(cid:18)p(cid:48)(u) (cid:19) (cid:35) (3)  ∂tu+∂xf(u) = ε1+α∂x ε2α −f(cid:48)(u)2 ∂xu +O(ε2). Therefore, as ε → 0 when α ∈ [0,1) we obtain the scalar conservation law  v = f(u),  (4)  ∂ u+∂ f(u) = 0. t x Note that, the main stability condition [13, 28] for system (3) corresponds to p(cid:48)(u) f(cid:48)(u)2 < , (5) ε2α and it is always satisfied in the limit ε → 0 when α > 0, whereas for α = 0 it requires suitable assumptions on the functions f(u) and p(u). 2 In classical kinetic theory the space-time scaling just discussed leads the so-called hydrody- namical limits of the Boltzmann equation (see [12], Chapter 11). For α = 0 this corresponds to the compressible Euler limit, whereas for α ∈ (0,1) the incompressible Euler limit is obtained. Something special happens when α = 1. In this case, in fact, to leading order in ε, we obtain the parabolic equation  v = f(u)−∂ p(u),  x (6)  ∂ u+∂ f(u) = ∂ p(u). t x xx In other words, considering larger times than those typical for Euler dynamics (α = 1 instead of α ∈ [0,1)), dissipative effects become non-negligible. This behavior characterizes the incom- pressible Navier-Stokes limit in classical kinetic theory. The development of numerical methods to solve hyperbolic systems with stiff source terms in the case α = 0 has been an active area of research in the past three decades [10, 18, 22, 18, 31, 33, 34, 32, 8]. Another series of works is concerned with the construction of robust schemes for α = 1 when a diffusion limit is obtained [21, 23, 25, 26]. However, very few papers have consideredthegeneralmultiscaleproblemoftype(1)forthevariouspossiblevaluesofα[29,24]. The common goal of this general class methods, often referred to as asymptotic-preserving (AP) schemes, was to obtain the macroscopic behavior described by the equilibrium system by solving the original relaxation system (1) with coarse grids ∆t, ∆x (cid:29) ε, where ∆t and ∆x are respectively the time step and the mesh size. Note that, since the characteristic speeds of the hyperbolic part of system (1) are of order 1/εα, most of the popular methods [10, 22, 32], for the solution to hyperbolic conservation laws withstiffrelaxationpresentseverallimitationswhenconsideringthewholerangeofα ∈ [0,1]and fail to capture the right behavior of the limit equilibrium equation unless the small relaxation rateisnumericallyresolved,leadingtoastabilityconditionoftheform∆t ∼ εα∆x. Clearly,this hyperbolic stiffness becomes very restrictive when α > 0, and for α = 1 in the parabolic regime ε (cid:28) ∆x where for an explicit scheme a parabolic time step restriction of the type ∆t ∼ ∆x2 is expected. A special class of IMEX schemes with explicit flux and implicit relaxation is able to deal with the parabolic relaxation (α = 1), [8]. Such methods, however, converge to an explicit scheme for the limit parabolic equation thus requiring a penalization technique to remove the final parabolic stiffness. In the present paper, using a reformulation of the problem, we develop high-order IMEX Runge-Kutta schemes for a system like (1) in the stiff regime which work uniformly with respect to the scaling parameter α. In the parabolic regime, α = 1, our approach gives a scheme which is not only consistent with (6) without resolving the small ε scale, but is also capable to avoid the parabolic stiffness leading to the CFL condition ∆t ∼ ∆x2. In other limiting regimes, that is to say when α ∈ [0,1), the scheme maintains all the nice properties of the numerical schemes for hyperbolic conservation laws, such as the ability to capture shocks with high resolution. Here, although the final schemes we develop will work independently on ε, we will mainly concentrate on the study of the stiff regime for system (1) that is to say when ε (cid:28) 1. Finally we emphasize that from a physical point of view the problem we consider here is close in spirit to the description of the macroscopic incompressible Navier-Stokes equations of fluid-dynamics by the detailed kinetic equations [4, 14, 12]. Although, for the sake of simplicity, we develop our theory for one-dimensional 2×2, systems the results extend far beyond these models. The rest of the paper is organized as follows. In Section 2 we recall some of the most popular IMEX Runge-Kutta approaches and emphasize their limited applicability when α ranges on the whole [0,1] interval. Next, in Section 3 we introduce our new approach with the aim to avoid 3 the stiffness induced by the characteristic speeds of system (1). First we present the simple first order scheme and then, using the IMEX Runge-Kutta formalism, we construct high order methods. In the case α = 1, these methods give rise to an explicit approximation of the limiting parabolic problem. In Section 4 we modify the previous schemes in order to avoid the limiting parabolic stiffness. In these schemes the diffusion term in the limiting equation is integrated implicitly. Finally in Section 5 several numerical examples are presented showing the robustness of the present approach. Some final considerations are contained in the last section and an appendix is also included. 2 Previous IMEX Runge-Kutta approaches Inthissection,tomotivatethenewapproach,werecallbrieflyotherwaystotacklestiffproblems through an implicit-explicit partitioning of the differential system. We discretize time first, and then we discretize space on the time discrete scheme. The moti- vation for not adopting a method of line approach is that we can choose the space discretization which is more suitable for each term. For simplicity of notation, in the sequel we assume that α = 1, and consider the hyperbolic-to-parabolic relaxation. Similar conclusions are obtained for α ∈ (0,1). Considerations on the case α = 0 are reported at the end of each subsection. 2.1 Additive approach Let us consider the simple implicit-explicit Euler method applied to (1) based on taking the fluxes explicitly and the stiff source implicitly [32] un+1−un = −vn, ∆t x (7) ε2vn+1−vn = −p(un) −(cid:0)vn+1−f(un+1)(cid:1). x ∆t Solving the second equation for vn+1 one obtains vn+1 = ε2 vn− ∆t (cid:0)p(un) −f(un+1)(cid:1). (8) ε2+∆t ε2+∆t x Making use of this relation (replacing n by n−1) in the first equation we get un+1−un + ε2 vn−1 = ∆t (cid:0)p(un−1) −f(un) (cid:1). ∆t ε2+∆t x ε2+∆t xx x Therefore, in the limit ε → 0, we a two levels scheme (in time) for problem (6) un+1−un = p(un−1) −f(un) . (9) xx x ∆t Although Eq. (9) is a consistent time discretization of the limit convection-diffusion equation (6), the presence of the term un−1 degrades the accuracy of the first order scheme. Furthermore, since as ε → 0 the equilibrium state of Eq. (8) vn+1 = f(un+1)−p(un) x involves two time levels, in general additional conditions on the explicit and implicit schemes are necessary in order to obtain asymptotic preserving high order methods. We refer to [8] for more details. These drawbacks are also present for any value of α ∈ (0,1]. Of course, since the additive approach has been originally designed to deal with the case α = 0 there are no problems in the regime. 4 2.2 Partitioned approach A different way to apply the implicit-explicit Euler method to (1) is based on taking the first equation explicitly and the second implicitly [5, 6, 30] un+1−un = −vn, ∆t x (10) ε2vn+1−vn = −p(un+1) −(cid:0)vn+1−f(un+1)(cid:1). x ∆t Solving for vn+1 the second equation we get vn+1 = ε2 vn− ∆t (cid:0)p(un+1) −f(un+1)(cid:1), (11) ε2+∆t ε2+∆t x which substituted into the first equation results in the two level scheme un+1−un ε2 ∆t + vn−1 = (p(un) −f(un) ). ∆t ε2+∆t x ε2+∆t xx x Of course, since the method has been developed specifically for the case α = 1, in the limit ε → 0 it yields a consistent explicit scheme for problem (6) un+1−un = p(un) −f(un) . (12) xx x ∆t It is easy to verify that this approach gives a consistent explicit scheme also for any value of α ∈ [0,1) in the hyperbolic limit. However, for non vanishingly small values of ε, the discretization of the fluxes at different time levels poses several difficulties. For example, if we consider for simplicity p(cid:48)(u) = 1, introducing the diagonal variables u±εv system (10) reads (un+1+εvn+1)−(un+εvn) = −1(un+1+εvn) − 1 (cid:0)vn+1−f(un+1)(cid:1), x ∆t ε ε (13) (un+1−εvn+1)−(un−εvn) = +1(un+1−εvn) + 1 (cid:0)vn+1−f(un+1)(cid:1), x ∆t ε ε and therefore the space derivatives of diagonal variables in (13) are defined as a combination of two time levels and it is then not clear how to discretize the system in characteristic variables. Furthermore, most numerical methods based on conservative variables evaluate the flux at the same time level. We refer to [6] for a discussion on these aspects and extensions to schemes that avoid the parabolic stiffness in the relaxation limit for α = 1. 2.3 Hybrid additive-partitioned approach A method which combines the advantages of the previous approaches in the various regimes has been proposed in [25, 26]. The idea is to take a convex combination of the previous schemes in such a way that we have an additive scheme in hyperbolic regimes and a partitioned one in parabolic ones. This is achieved by taking un+1−un = −vn, ∆t x (14) ε2vn+1−vn = −φ(ε)p(un) −(1−φ(ε))p(un+1) −(cid:0)vn+1−f(un+1)(cid:1), x x ∆t 5 where 0 ≤ φ(ε) ≤ 1 is such that φ(ε) ≈ 1 for ε = O(1) and φ(ε) ≈ 0 when ε (cid:28) 1. For example φ(ε) = min{ε2,1} or the smoother approximation φ(ε) = tanh(ε2) were considered in [23, 25]. In the limit ε → 0 we have the asymptotic behavior of the partitioned approach, whereas for largervaluesofεwehavetheusualadditiveapproach. Themethodcanbenaturallyextendedto thegeneralmultiscalecaseandprovidesaconsistentdiscretizationalsoforanyvalueofα ∈ [0,1). We refer to [24, 29] for further details. One of the main advantages of this approach is that it results in a convex approximation at different times of the space derivative p(u) , therefore x one can use different space discretizations for the derivative appearing at time n and the one at time n+1. Typically, at time n one can choose a standard hyperbolic discretization that works for large values of ε whereas at time n+1 classical central difference schemes that works in the parabolic limit suffices to avoid a CFL condition of the type ∆t = O(ε). In this sense the function φ(ε) can be also interpreted as an interpolation parameter between different fluxes in the evaluation of the space derivatives in (10). The determination of the optimal expression of φ(ε) in terms of stability and accuracy is a delicate aspect which is beyond the scope of the present paper. 3 A unified IMEX Runge-Kutta approach In this section we present a different approach which overcomes some of the drawbacks of the above mentioned methods. 3.1 Description of the method Again let us consider initially, for simplicity of notation, the case α = 1. We now consider the following discretization for system (1) un+1−un = −vn+1, ∆t x (15) ε2vn+1−vn = −(cid:0)p(un) +vn+1−f(un)(cid:1). x ∆t Solving the second equation for vn+1 one obtains ε2 ∆t vn+1 = vn− (p(un) −f(un)). (16) ε2+∆t ε2+∆t x Making use of this relation in the first equation we get un+1−un ε2 ∆t + vn = (p(un) −f(un) ). ∆t ε2+∆t x ε2+∆t xx x Note that, at variance with all the previous approaches, the first equation now uses only two time levels and the space derivatives appear all at the same time level n. Therefore we can rewrite the scheme in the equivalent fully explicit form un+1−un ε2 ∆t ∆t + vn+ f(un) = p(un) , ∆t ε2+∆t x ε2+∆t x ε2+∆t xx (17) vn+1−vn 1 1 + p(un) = − (vn−f(un)). ∆t ε2+∆t x ε2+∆t 6 In particular, for small values of ∆t the scheme (17) corresponds to the system ε2 ∆t ∆t u + v + f(u) = p(u) +O(∆t), t ε2+∆t x ε2+∆t x ε2+∆t xx (18) 1 1 v + p(u) = − (v−f(u))+O(∆t), t ε2+∆t x ε2+∆t where we only used un+1−un vn+1−vn = u +O(∆t), = v +O(∆t) t t ∆t ∆t leaving all other terms. Note that the left part of system (18) is hyperbolic with characteristic speeds 1 (cid:32) (cid:115) (cid:33) ξ 4ε2 λ1(∆t,ε) = 1 c± c2+ , (19) ± 2 (∆t)2 with ξ = ∆t/(ε2+∆t) ∈ (0,1), c = f(cid:48)(u) and for simplicity we have set p(cid:48)(u) = 1. 1 Now, if ∆t → 0 for a fixed ε, system (18) converges to the original system (1) for α = 1 and by (19), the characteristics speeds converge to the usual ones, i.e. 1 λ1(0,ε) = ± . ± ε Ontheotherhand,forafixed∆t,thecharacteristicspeedsλ1 andλ1 in(19)arerespectively + − decreasing and increasing functions of ε and as ε → 0 they converge to 1 λ1(∆t,0) = (c±|c|). (20) ± 2 Therefore, if we denote by ∆x the space discretization parameter, we obtain the expected hyperbolic CFL condition ∆t ≤ ∆x/|c|. Concerning the second order derivative on the right hand side of (18) this induces a stability restriction of type ∆t ∼ (∆x)2/ξ . In particular, since the discrete system (17) as ε → 0 relaxes 1 towards un+1−un +f(un) = p(un) , (21) x xx ∆t in such a limit we have the natural stability condition which links the time step to the square of the mesh space ∆t ∼ (∆x)2. Similarly, in the case α ∈ [0,1), for small values of ∆t the scheme applied to (1) corresponds to the system ε1+α ∆t ∆tε1−α u + v + f(u) = p(u) +O(∆t), t ε1+α+∆t x ε1+α+∆t x ε1+α+∆t xx (22) ε1−α 1 v + p(u) = − (v−f(u))+O(∆t). t ε1+α+∆t x ε1+α+∆t The left-hand side in (22) now has characteristic speeds (cid:32) (cid:115) (cid:33) ξ 4ε2 λα(∆t,ε) = α c± c2+ , (23) ± 2 (∆t)2 1The subscript 1 in the next expression indicates α=1. 7 with ξ = ∆t/(ε1+α+∆t) ∈ [0,1). As before, if we fix ε and send ∆t → 0 we obtain the usual α characteristic speeds 1 λα(0,ε) = ± . ± εα The limit behavior of the discrete system now is given by un+1−un +f(un) = 0, (24) x ∆t and, similarly to the analysis for α = 1, we now observe that the characteristic speed do not diverge as ε → 0 since we have 1 λα(∆t,0) = (c±|c|). ± 2 Therefore, we obtain the expected hyperbolic CFL condition ∆t ≤ ∆x/|c|, coming from the hyperbolic part of the system. As before, the stability restriction coming from the parabolic term, is ∆t ∼ ∆x2/ξ . α 3.2 Extension to general IMEX Runge-Kutta schemes Now we extend the analysis just performed for the simple first order implicit-explicit Euler scheme to a general IMEX-RK scheme. 3.2.1 Notations An IMEX-RK scheme can be represented with a double Butcher tableau c˜ A˜ c A Explicit : Implicit : (25) ˜bT bT where A˜ = (a˜ ) is an s×s lower triangular matrix with zero diagonal entries for an explicit ij scheme, and, since computational efficiency in most cases is of paramount importance, usually thes×smatrixA = (a )foranimplicitschemeisrestrictedtotheparticularclassofdiagonally ij implicit Runge-Kutta (DIRK) methods, i.e., (a = 0, for j > i). In fact, the use of a DIRK ij schemeisenoughtoensurethattheexplicitpartoftheschemetermisalwaysevaluatedexplicitly (see [1], [11], [7]). The coefficients c˜ and c are given by the usual relation c˜ = (cid:80)i−1 a˜ , c = (cid:80)i a and i j=1 ij i j=1 ij the vectors ˜b = (˜b ) and b = (b ) provide the quadrature weights to combine internal i i=1···s i i=1···s stages of the RK method. The order conditions can be derived as usual matching the Taylor expansion of the exact and numerical solution, we refer to [1, 11, 32] for more details. Let us mention that from a practical viewpoint, coupling conditions becomes rather severe if one is interested in very high order schemes (say higher then third). Here we recall the order conditions for IMEX-RK schemes up to order p = 2 and p = 3, under the assumption that c˜= c. firstorder (consistency) ˜bTe = 1, bTe = 1. second order ˜bc = 1/2, bTc = 1/2. (26) third order ˜bc2 = 1/3, bTc2 = 1/3, bTA˜c = 1/6, ˜bTA˜c = 1/6, bTAc = 1/6, ˜bTAc = 1/6. 8 where we denote by eT = (1,··· ,1) ∈ Rs, and from the previous relaxation for c˜and c, we have A˜e = c˜and Ae = c. It is useful to characterize the different IMEX schemes we consider in this paper according to the structure of the DIRK method. Following [2] we have Definition 1 1. We call an IMEX-RK method of type I or type A (see [32]) if the matrix A ∈ Rs×s is invertible, or equivalently a (cid:54)= 0, i = 1,...,s. ii 2. We call an IMEX-RK method of type II or type CK (see [11]) if the matrix A can be written as (cid:18) (cid:19) 0 0 A = , (27) a Aˆ with a = (a ,...,a )T ∈ R(s−1) and the submatrix Aˆ ∈ R(s−1) × (s−1) is invertible, or 21 s1 equivalently a (cid:54)= 0, i = 2,...,s. In the special case a = 0, w = 0 the scheme is said to be ii 1 of type ARS (see [1]) and the DIRK method is reducible to a method using s−1 stages. Thefollowingdefinitionwillbealsousefultocharacterizethepropertiesofthemethod[6,8]. Definition 2 WecallanIMEX-RKmethod implicitlystifflyaccurate(ISA)ifthecorresponding DIRK method is stiffly accurate, namely a = b , i = 1,...,s. (28) si i If in addition the explicit methods satisfies a˜ =˜b , i = 1,...,s−1 (29) si i the IMEX-RK method is said to be globally stiffly accurate (GSA). ThedefinitionsofISAfollowsnaturallyfromthefactthats-stageimplicitRunge-Kuttamethods for which a = b for i = 1,··· ,ν are called stiffly accurate (see [19] for details). A ν-stage si i explicit Runge-Kutta method for which a˜ =˜b for i = 1,···s−1 is called FSAL (First Same si i As Last, see [20] for details). Note that FSAL methods have the advantage that they require only s−1 function evaluations per time step, because the last stage of step n coincides with the first stage of step n+1. Therefore, an IMEX-RK scheme is globally stiffly accurate if the implicit scheme is stiffly accurate and the explicit scheme is FSAL. We observe that this definition states also that the numerical solution of a GSA IMEX-RK scheme coincides exactly with the last internal stage of the scheme. 3.2.2 The unified IMEX-RK setting Let us consider system (1) and again for simplicity of notation, we set α = 1. The unified IMEX-RK approach that generalizes first order scheme (15) corresponds to compute first the internal stages U and V U = une−∆tAV x ∆t ∆t (30) V = vne− A˜(p(U) −f(U))− AV, ε2 x ε2 9 and then the numerical solution: un+1 = un−∆tbTV x ∆t ∆t (31) vn+1 = vn− ˜bT(p(U) −f(U))− bTV, ε2 x ε2 where f(U) and p(U) are the vectors with component f(U) = f(Ui) and p(U) = p(Ui) respec- i i tively. Solving the second equation for V in (30) with ζ = ε2/∆t, one obtains (cid:16) (cid:17) V = (ζI +A)−1 ζvne−A˜(p(U) −f(U)) . (32) x Substituting this relation in the first equation we have the resulting IMEX-RK scheme may be written as U −une +ζA(ζI +A)−1evn+A(ζI +A)−1A˜f(U) = A(ζI +A)−1A˜p(U) ∆t x x xx (33) V −vne 1 1 (cid:16) (cid:17) + A˜p(U) = − AV −A˜f(U) ∆t ε2 x ε2 and un+1−un +ζbT(ζI +A)−1evn+bT(ζI +A)−1A˜f(U) = bT(ζI +A)−1A˜p(U) ∆t x x xx (34) vn+1−vn 1 1 (cid:16) (cid:17) + ˜bTp(U) = − bTV −˜bTf(U) . ∆t ε2 x ε2 Then setting B = (ζI +A)−1 for small values of ∆t, from (34), we get the system u +ζbTBev +bTBA˜f(U) = bTBA˜p(U) +O(∆t) t x x xx 1 1 (cid:16) (cid:17) (35) v + ˜bTp(U) = − bTV −˜bTf(U) +O(∆t). t ε2 x ε2 Now,werewritef(U) = f(cid:48)(U)U andp(U)x = p(cid:48)(U)U ,wheref(cid:48)(U)andp(cid:48)(U)arediagonal x x x matricies with elements f(cid:48)(U) = f(cid:48)(Ui) and p(cid:48)(U) = p(cid:48)(Ui) respectively. Furthermore, for ii ii simplicity we assume f(cid:48)(Ui) = c and p(cid:48)(Ui) = 1. From the IMEX-RK stages (33) we have U = u e−ε2ABevn−∆tcABA˜U +∆tABA˜U , (36) n x x xx and using the fact that U = une+O(∆t) and bTe = 1 (consistency of the IMEX-RK scheme) we can write the hyperbolic part in (33) as u +ζbTBev +cbTBA˜eu = O(∆t) t x x 1 (37) v + u = O(∆t). t ε2 x By computing the eigenvalues of the hyperbolic part we obtain (cid:32) (cid:114) (cid:33) 1 (cid:16) (cid:17)2 bTBe Λ1(∆t,ε) = cbTBA˜e± cbTBA˜e +4 . (38) ± 2 ∆t Next, we show that the above characteristic speeds are limited. Here we need to assume that the scheme is GSA. In fact, (cid:40) 1, for ζ → 0 bTBe = bT(ζI +A)−1e = 1 (39) ∼ , for ζ → ∞ ζ 10