A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation
AA unified IMEX Runge-Kutta approach for hyperbolic systemswith multiscale relaxation
S. Boscarino ∗ L. Pareschi † G.Russo ‡ January 17, 2017
Abstract
In this paper we consider the development of Implicit-Explicit (IMEX) Runge-Kuttaschemes for hyperbolic systems with multiscale relaxation. In such systems the scalingdepends on an additional parameter which modifies the nature of the asymptotic behaviorwhich can be either hyperbolic or parabolic. Because of the multiple scalings, standardIMEX Runge-Kutta methods for hyperbolic systems with relaxation loose their efficiencyand a different approach should be adopted to guarantee asymptotic preservation in stiffregimes. We show that the proposed approach is capable to capture the correct asymptoticlimit of the system independently of the scaling used. Several numerical examples confirmour theoretical analysis.
Key words.
IMEX Runge-Kutta methods, hyperbolic conservation laws with sources, diffusionequations, hydrodynamic limits, stiff systems, asymptotic-preserving schemes.
AMS subject classification.
Contents ∗ 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] ). a r X i v : . [ m a t h . NA ] J a n Conclusions 23A Appendix: IMEX-RK schemes 23
Hyperbolic systems with relaxation often contain multiple space-time scales which may differof 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-dynamicequations, 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 thefollowing ∂ t u + ∂ x v = 0 ,∂ t v + 1 ε α ∂ x p ( u ) = − ε α ( v − f ( u )) , α ∈ [0 ,
1] (1)where p (cid:48) ( u ) >
0. System (1) is hyperbolic with two distinct real characteristics speeds ± (cid:112) p (cid:48) ( u ) /ε α .Note that the scaling introduced in (1) corresponds to the study of the limiting behavior ofthe solution for the usual hyperbolic system with a singular perturbation source ∂ τ u + ∂ ξ V = 0 ,∂ τ V + ∂ ξ p ( u ) = − ε ( V − F ( u )) , (2)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 arelooking at larger microscopic times.In particular, for small values of ε , using the Chapman-Enskog expansion, the behavior ofthe solution to (1) is, at least formally, governed by the following nonlinear parabolic equation v = f ( u ) − ε − α ∂ x p ( u ) + ε α f (cid:48) ( u ) ∂ x u + O ( ε ) ,∂ t u + ∂ x f ( u ) = ε α ∂ x (cid:34) (cid:18) p (cid:48) ( u ) ε α − f (cid:48) ( u ) (cid:19) ∂ x u (cid:35) + O ( ε ) . (3)Therefore, as ε → α ∈ [0 ,
1) we obtain the scalar conservation law v = f ( u ) ,∂ t u + ∂ x f ( u ) = 0 . (4)Note that, the main stability condition [13, 28] for system (3) corresponds to f (cid:48) ( u ) < p (cid:48) ( u ) ε α , (5)and it is always satisfied in the limit ε → α >
0, whereas for α = 0 it requires suitableassumptions on the functions f ( u ) and p ( u ). 2n 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 tothe 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 obtainthe parabolic equation v = f ( u ) − ∂ x p ( u ) ,∂ t u + ∂ x f ( u ) = ∂ xx p ( u ) . (6)In other words, considering larger times than those typical for Euler dynamics ( α = 1 insteadof α ∈ [0 , incom-pressible Navier-Stokes limit in classical kinetic theory.The development of numerical methods to solve hyperbolic systems with stiff source termsin 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 schemesfor α = 1 when a diffusion limit is obtained [21, 23, 25, 26]. However, very few papers haveconsidered the general multiscale problem of type (1) for the various possible values of α [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 bysolving the original relaxation system (1) with coarse grids ∆ t, ∆ x (cid:29) ε , where ∆ t and ∆ x arerespectively the time step and the mesh size.Note that, since the characteristic speeds of the hyperbolic part of system (1) are of order1 /ε α , most of the popular methods [10, 22, 32], for the solution to hyperbolic conservation lawswith stiff relaxation present several limitations when considering the whole range of α ∈ [0 ,
1] andfail to capture the right behavior of the limit equilibrium equation unless the small relaxationrate is numerically resolved, leading to a stability condition of the form ∆ 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 ∼ ∆ x isexpected. A special class of IMEX schemes with explicit flux and implicit relaxation is able todeal with the parabolic relaxation ( α = 1), [8]. Such methods, however, converge to an explicitscheme for the limit parabolic equation thus requiring a penalization technique to remove thefinal parabolic stiffness.In the present paper, using a reformulation of the problem, we develop high-order IMEXRunge-Kutta schemes for a system like (1) in the stiff regime which work uniformly with respectto the scaling parameter α . In the parabolic regime, α = 1, our approach gives a scheme whichis not only consistent with (6) without resolving the small ε scale, but is also capable to avoidthe parabolic stiffness leading to the CFL condition ∆ t ∼ ∆ x . In other limiting regimes, thatis to say when α ∈ [0 , ε , we will mainlyconcentrate on the study of the stiff regime for system (1) that is to say when ε (cid:28) ×
2, systems the results extend far beyond thesemodels.The rest of the paper is organized as follows. In Section 2 we recall some of the most popularIMEX Runge-Kutta approaches and emphasize their limited applicability when α ranges on thewhole [0 ,
1] interval. Next, in Section 3 we introduce our new approach with the aim to avoid3he stiffness induced by the characteristic speeds of system (1). First we present the simplefirst order scheme and then, using the IMEX Runge-Kutta formalism, we construct high ordermethods. In the case α = 1, these methods give rise to an explicit approximation of the limitingparabolic problem. In Section 4 we modify the previous schemes in order to avoid the limitingparabolic stiffness. In these schemes the diffusion term in the limiting equation is integratedimplicitly. Finally in Section 5 several numerical examples are presented showing the robustnessof the present approach. Some final considerations are contained in the last section and anappendix is also included. In this section, to motivate the new approach, we recall briefly other ways to tackle stiff problemsthrough 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 discretizationwhich 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 , α = 0 are reported at the end of each subsection. Let us consider the simple implicit-explicit Euler method applied to (1) based on taking thefluxes explicitly and the stiff source implicitly [32] u n +1 − u n ∆ t = − v nx ,ε v n +1 − v n ∆ t = − p ( u n ) x − (cid:0) v n +1 − f ( u n +1 ) (cid:1) . (7)Solving the second equation for v n +1 one obtains v n +1 = ε ε + ∆ t v n − ∆ tε + ∆ t (cid:0) p ( u n ) x − f ( u n +1 ) (cid:1) . (8)Making use of this relation (replacing n by n −
1) in the first equation we get u n +1 − u n ∆ t + ε ε + ∆ t v n − x = ∆ tε + ∆ t (cid:0) p ( u n − ) xx − f ( u n ) x (cid:1) . Therefore, in the limit ε →
0, we a two levels scheme (in time) for problem (6) u n +1 − u n ∆ t = p ( u n − ) xx − f ( u n ) x . (9)Although Eq. (9) is a consistent time discretization of the limit convection-diffusion equation(6), the presence of the term u n − degrades the accuracy of the first order scheme.Furthermore, since as ε → v n +1 = f ( u n +1 ) − p ( u n ) x involves two time levels, in general additional conditions on the explicit and implicit schemesare necessary in order to obtain asymptotic preserving high order methods. We refer to [8] formore details. These drawbacks are also present for any value of α ∈ (0 , α = 0 there are noproblems in the regime. 4 .2 Partitioned approach A different way to apply the implicit-explicit Euler method to (1) is based on taking the firstequation explicitly and the second implicitly [5, 6, 30] u n +1 − u n ∆ t = − v nx ,ε v n +1 − v n ∆ t = − p ( u n +1 ) x − (cid:0) v n +1 − f ( u n +1 ) (cid:1) . (10)Solving for v n +1 the second equation we get v n +1 = ε ε + ∆ t v n − ∆ tε + ∆ t (cid:0) p ( u n +1 ) x − f ( u n +1 ) (cid:1) , (11)which substituted into the first equation results in the two level scheme u n +1 − u n ∆ t + ε ε + ∆ t v n − x = ∆ tε + ∆ t ( p ( u n ) xx − f ( u n ) x ) . Of course, since the method has been developed specifically for the case α = 1, in the limit ε → u n +1 − u n ∆ t = p ( u n ) xx − f ( u n ) x . (12)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 discretizationof the fluxes at different time levels poses several difficulties. For example, if we consider forsimplicity p (cid:48) ( u ) = 1, introducing the diagonal variables u ± εv system (10) reads( u n +1 + εv n +1 ) − ( u n + εv n )∆ t = − ε ( u n +1 + εv n ) x − ε (cid:0) v n +1 − f ( u n +1 ) (cid:1) , ( u n +1 − εv n +1 ) − ( u n − εv n )∆ t = + 1 ε ( u n +1 − εv n ) x + 1 ε (cid:0) v n +1 − f ( u n +1 ) (cid:1) , (13)and therefore the space derivatives of diagonal variables in (13) are defined as a combination oftwo 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 thesame time level. We refer to [6] for a discussion on these aspects and extensions to schemes thatavoid the parabolic stiffness in the relaxation limit for α = 1. A method which combines the advantages of the previous approaches in the various regimes hasbeen proposed in [25, 26]. The idea is to take a convex combination of the previous schemesin such a way that we have an additive scheme in hyperbolic regimes and a partitioned one inparabolic ones. This is achieved by taking u n +1 − u n ∆ t = − v nx ,ε v n +1 − v n ∆ t = − φ ( ε ) p ( u n ) x − (1 − φ ( ε )) p ( u n +1 ) x − (cid:0) v n +1 − f ( u n +1 ) (cid:1) , (14)5here 0 ≤ φ ( ε ) ≤ φ ( ε ) ≈ ε = O (1) and φ ( ε ) ≈ ε (cid:28)
1. For example φ ( ε ) = min { ε , } or the smoother approximation φ ( ε ) = tanh( ε ) were considered in [23, 25].In the limit ε → ε we have the usual additive approach. The method can be naturally extended tothe general multiscale case and provides a consistent discretization also for any value of α ∈ [0 , p ( u ) x , thereforeone can use different space discretizations for the derivative appearing at time n and the oneat time n + 1. Typically, at time n one can choose a standard hyperbolic discretization thatworks for large values of ε whereas at time n + 1 classical central difference schemes that worksin the parabolic limit suffices to avoid a CFL condition of the type ∆ t = O ( ε ). In this sensethe function φ ( ε ) can be also interpreted as an interpolation parameter between different fluxesin the evaluation of the space derivatives in (10). The determination of the optimal expressionof φ ( ε ) in terms of stability and accuracy is a delicate aspect which is beyond the scope of thepresent paper. In this section we present a different approach which overcomes some of the drawbacks of theabove mentioned methods.
Again let us consider initially, for simplicity of notation, the case α = 1.We now consider the following discretization for system (1) u n +1 − u n ∆ t = − v n +1 x ,ε v n +1 − v n ∆ t = − (cid:0) p ( u n ) x + v n +1 − f ( u n ) (cid:1) . (15)Solving the second equation for v n +1 one obtains v n +1 = ε ε + ∆ t v n − ∆ tε + ∆ t ( p ( u n ) x − f ( u n )) . (16)Making use of this relation in the first equation we get u n +1 − u n ∆ t + ε ε + ∆ t v nx = ∆ tε + ∆ t ( p ( u n ) xx − f ( u n ) x ) . Note that, at variance with all the previous approaches, the first equation now uses only twotime levels and the space derivatives appear all at the same time level n . Therefore we canrewrite the scheme in the equivalent fully explicit form u n +1 − u n ∆ t + ε ε + ∆ t v nx + ∆ tε + ∆ t f ( u n ) x = ∆ tε + ∆ t p ( u n ) xx ,v n +1 − v n ∆ t + 1 ε + ∆ t p ( u n ) x = − ε + ∆ t ( v n − f ( u n )) . (17)6n particular, for small values of ∆ t the scheme (17) corresponds to the system u t + ε ε + ∆ t v x + ∆ tε + ∆ t f ( u ) x = ∆ tε + ∆ t p ( u ) xx + O (∆ t ) ,v t + 1 ε + ∆ t p ( u ) x = − ε + ∆ t ( v − f ( u )) + O (∆ t ) , (18)where we only used u n +1 − u n ∆ t = u t + O (∆ t ) , v n +1 − v n ∆ t = v t + O (∆ t )leaving all other terms. Note that the left part of system (18) is hyperbolic with characteristicspeeds λ ± (∆ t, ε ) = ξ (cid:32) c ± (cid:115) c + 4 ε (∆ t ) (cid:33) , (19)with ξ = ∆ t/ ( ε + ∆ t ) ∈ (0 , c = f (cid:48) ( u ) and for simplicity we have set p (cid:48) ( u ) = 1.Now, if ∆ t → ε , system (18) converges to the original system (1) for α = 1 andby (19), the characteristics speeds converge to the usual ones, i.e. λ ± (0 , ε ) = ± ε . On the other hand, for a fixed ∆ t , the characteristic speeds λ and λ − in (19) are respectivelydecreasing and increasing functions of ε and as ε → λ ± (∆ t,
0) = 12 ( c ± | c | ) . (20)Therefore, if we denote by ∆ x the space discretization parameter, we obtain the expectedhyperbolic CFL condition ∆ t ≤ ∆ x/ | c | .Concerning the second order derivative on the right hand side of (18) this induces a stabilityrestriction of type ∆ t ∼ (∆ x ) /ξ . In particular, since the discrete system (17) as ε → u n +1 − u n ∆ t + f ( u n ) x = p ( u n ) xx , (21)in such a limit we have the natural stability condition which links the time step to the squareof the mesh space ∆ t ∼ (∆ x ) .Similarly, in the case α ∈ [0 , t the scheme applied to (1) correspondsto the system u t + ε α ε α + ∆ t v x + ∆ tε α + ∆ t f ( u ) x = ∆ t ε − α ε α + ∆ t p ( u ) xx + O (∆ t ) ,v t + ε − α ε α + ∆ t p ( u ) x = − ε α + ∆ t ( v − f ( u )) + O (∆ t ) . (22)The left-hand side in (22) now has characteristic speeds λ α ± (∆ t, ε ) = ξ α (cid:32) c ± (cid:115) c + 4 ε (∆ t ) (cid:33) , (23) The subscript 1 in the next expression indicates α = 1. ξ α = ∆ t/ ( ε α + ∆ t ) ∈ [0 , ε and send ∆ t → λ α ± (0 , ε ) = ± ε α . The limit behavior of the discrete system now is given by u n +1 − u n ∆ t + f ( u n ) x = 0 , (24)and, similarly to the analysis for α = 1, we now observe that the characteristic speed do notdiverge as ε → λ α ± (∆ t,
0) = 12 ( c ± | c | ) . Therefore, we obtain the expected hyperbolic CFL condition ∆ t ≤ ∆ x/ | c | , coming from thehyperbolic part of the system. As before, the stability restriction coming from the parabolicterm, is ∆ t ∼ ∆ x /ξ α . Now we extend the analysis just performed for the simple first order implicit-explicit Eulerscheme to a general IMEX-RK scheme.
An IMEX-RK scheme can be represented with a double Butcher tableauExplicit : ˜ c ˜ A ˜ b T Implicit : c Ab T (25)where ˜ A = (˜ a ij ) is an s × s lower triangular matrix with zero diagonal entries for an explicitscheme, and, since computational efficiency in most cases is of paramount importance, usuallythe s × s matrix A = ( a ij ) for an implicit scheme is restricted to the particular class of diagonallyimplicit Runge-Kutta (DIRK) methods, i.e., ( a ij = 0 , for j > i ). In fact, the use of a DIRKscheme is enough to ensure that the explicit part of the scheme term is always evaluated explicitly(see [1], [11], [7]).The coefficients ˜ c and c are given by the usual relation ˜ c i = (cid:80) i − j =1 ˜ a ij , c i = (cid:80) ij =1 a ij andthe vectors ˜ b = (˜ b i ) i =1 ··· s and b = ( b i ) i =1 ··· s provide the quadrature weights to combine internalstages of the RK method.The order conditions can be derived as usual matching the Taylor expansion of the exactand numerical solution, we refer to [1, 11, 32] for more details. Let us mention that from apractical viewpoint, coupling conditions becomes rather severe if one is interested in very highorder 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) ˜ b T e = 1 , b T e = 1 . second order ˜ bc = 1 / , b T c = 1 / . third order ˜ bc = 1 / , b T c = 1 / ,b T ˜ Ac = 1 / , ˜ b T ˜ Ac = 1 / , b T Ac = 1 / , ˜ b T Ac = 1 / . (26)8here we denote by e T = (1 , · · · , ∈ R s , and from the previous relaxation for ˜ c and c , we have˜ Ae = ˜ c and Ae = c .It is useful to characterize the different IMEX schemes we consider in this paper accordingto 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 ∈ R s × s isinvertible, or equivalently a ii (cid:54) = 0 , i = 1 , . . . , s .2. We call an IMEX-RK method of type II or type CK (see [11]) if the matrix A can bewritten as A = (cid:18) a ˆ A (cid:19) , (27) with a = ( a , . . . , a s ) T ∈ R ( s − and the submatrix ˆ A ∈ R ( s − × ( s − is invertible, orequivalently a ii (cid:54) = 0 , i = 2 , . . . , s . In the special case a = 0 , w = 0 the scheme is said to beof type ARS (see [1]) and the DIRK method is reducible to a method using s − stages. The following definition will be also useful to characterize the properties of the method [6, 8].
Definition 2
We call an IMEX-RK method implicitly stiffly accurate (ISA) if the correspondingDIRK method is stiffly accurate , namely a si = b i , i = 1 , . . . , s. (28) If in addition the explicit methods satisfies ˜ a si = ˜ b i , i = 1 , . . . , s − the IMEX-RK method is said to be globally stiffly accurate (GSA) . The definitions of ISA follows naturally from the fact that s -stage implicit Runge-Kutta methodsfor which a si = b i for i = 1 , · · · , ν are called stiffly accurate (see [19] for details). A ν -stageexplicit Runge-Kutta method for which ˜ a si = ˜ b i for i = 1 , · · · s − First SameAs Last , see [20] for details). Note that FSAL methods have the advantage that they requireonly s − n coincides with thefirst stage of step n + 1.Therefore, an IMEX-RK scheme is globally stiffly accurate if the implicit scheme is stifflyaccurate and the explicit scheme is FSAL. We observe that this definition states also that thenumerical solution of a GSA IMEX-RK scheme coincides exactly with the last internal stage ofthe scheme. Let us consider system (1) and again for simplicity of notation, we set α = 1. The unifiedIMEX-RK approach that generalizes first order scheme (15) corresponds to compute first theinternal stages U and V U = u n e − ∆ tAV x V = v n e − ∆ tε ˜ A ( p ( U ) x − f ( U )) − ∆ tε AV, (30)9nd then the numerical solution: u n +1 = u n − ∆ tb T V x v n +1 = v n − ∆ tε ˜ b T ( p ( U ) x − f ( U )) − ∆ tε b T V, (31)where f ( U ) and p ( U ) are the vectors with component f ( U ) i = f ( U i ) and p ( U ) i = p ( U i ) respec-tively. Solving the second equation for V in (30) with ζ = ε / ∆ t , one obtains V = ( ζI + A ) − (cid:16) ζv n e − ˜ A ( p ( U ) x − f ( U )) (cid:17) . (32)Substituting this relation in the first equation we have the resulting IMEX-RK scheme may bewritten as U − u n e ∆ t + ζA ( ζI + A ) − ev nx + A ( ζI + A ) − ˜ Af ( U ) x = A ( ζI + A ) − ˜ Ap ( U ) xx V − v n e ∆ t + 1 ε ˜ Ap ( U ) x = − ε (cid:16) AV − ˜ Af ( U ) (cid:17) (33)and u n +1 − u n ∆ t + ζb T ( ζI + A ) − ev nx + b T ( ζI + A ) − ˜ Af ( U ) x = b T ( ζI + A ) − ˜ Ap ( U ) xx v n +1 − v n ∆ t + 1 ε ˜ b T p ( U ) x = − ε (cid:16) b T V − ˜ b T f ( U ) (cid:17) . (34)Then setting B = ( ζI + A ) − for small values of ∆ t , from (34), we get the system u t + ζb T Bev x + b T B ˜ Af ( U ) x = b T B ˜ Ap ( U ) xx + O (∆ t ) v t + 1 ε ˜ b T p ( U ) x = − ε (cid:16) b T V − ˜ b T f ( U ) (cid:17) + O (∆ t ) . (35)Now, we rewrite f ( U ) x = f (cid:48) ( U ) U x and p ( U ) x = p (cid:48) ( U ) U x , where f (cid:48) ( U ) and p (cid:48) ( U ) are diagonalmatricies with elements f (cid:48) ( U ) ii = f (cid:48) ( U i ) and p (cid:48) ( U ) ii = p (cid:48) ( U i ) respectively. Furthermore, forsimplicity we assume f (cid:48) ( U i ) = c and p (cid:48) ( U i ) = 1.From the IMEX-RK stages (33) we have U = u n e − ε ABev nx − ∆ tcAB ˜ AU x + ∆ tAB ˜ AU xx , (36)and using the fact that U = u n e + O (∆ t ) and b T e = 1 (consistency of the IMEX-RK scheme)we can write the hyperbolic part in (33) as u t + ζb T Bev x + cb T B ˜ Aeu x = O (∆ t ) v t + 1 ε u x = O (∆ t ) . (37)By computing the eigenvalues of the hyperbolic part we obtainΛ ± (∆ t, ε ) = 12 (cid:32) cb T B ˜ Ae ± (cid:114)(cid:16) cb T B ˜ Ae (cid:17) + 4 b T Be ∆ t (cid:33) . (38)Next, we show that the above characteristic speeds are limited. Here we need to assume thatthe scheme is GSA. In fact, b T Be = b T ( ζI + A ) − e = (cid:40) , for ζ → ∼ ζ , for ζ → ∞ (39)10here for the first case using the ISA property we get b T A − e = e Ts e = 1 (here e Ts = (0 , · · · , , ∈ R s ), whereas for the second case we note that for the consistency of the scheme b T e = 1. Onthe other hand, b T B ˜ Ae = b T ( ζI + A ) − ˜ Ae = (cid:40) , for ζ → ∼ ζ , for ζ → ∞ (40)where for the first case using the GSA property we get b T A − ˜ Ae = e Ts ˜ Ae = ˜ b T e = 1, whereas forthe second case we note that the quantity b T ˜ Ae is a number and if we assume that the schemeis a second order accurate one, this gives b T ˜ Ae = 1 /
2, by (26).As a consequence we haveΛ ± (∆ t,
0) = 12 (cid:32) c ± (cid:114) c + 4∆ t (cid:33) , Λ ± (0 , ε ) = ± ε , (41)and therefore the CFL condition, in the limit ε →
0, becomes∆ t ≤ (cid:32) | c | (cid:114) c t (cid:33) − ∆ x = 2 √ ∆ t | c |√ ∆ t + √ c ∆ t + 4 ∆ x. (42)Now, if ∆ t (cid:28) √ ∆ t ≤ | c |√ ∆ t + √ c ∆ t + 4 ∆ x, the parabolic time step restriction ∆ t ∼ ∆ x .Finally, we prove that the scheme (33)-(34) in the limit case ε → V = ζA − v n e − A − ( I − ζA − ) ˜ A ( f ( U ) − p ( U ) x ) + O ( ζ ) , and substituting in the numerical solution v n +1 we obtain ζv n +1 = ζ (1 − b T A − e ) v n + ( b T A − ˜ A − ˜ b T )( U x − f ( U )) − ζb T A − ˜ A ( f ( U ) − p ( U ) x ) + O ( ζ ) . Consistency as ζ →
0, implies ( b T A − ˜ A − ˜ b T ) = 0, which is satisfied if the scheme is GSA, becausein this case b T = e Ts A , ˜ b T = e Ts ˜ A , therefore 1 − b T A − e = 0 and b T A − ˜ A − ˜ b T = e Ts ˜ A − ˜ b T = 0.Then we have v n +1 = − b T A − ˜ A ( f ( U ) − p ( U ) x ) = e Ts A − ˜ A ( f ( U ) − p ( U ) x )with V = − A − ˜ A ( f ( U ) − p ( U ) x ) . Similarly from (52) we have for the U internal stages U − u n ∆ t = − ζv nx e + ˜ A ( p ( U ) x − f ( U )) x − ζA − ˜ A ( p ( U ) x − f ( U )) + O ( ζ )11nd for the numerical solution u n +1 − u n ∆ t = − ζb T v nx e + b T A − ˜ A ( p ( U ) x − f ( U )) x − ζb T A − ˜ A ( p ( U ) x − f ( U )) x + O ( ζ ) , therefore, this leads to U = u n e + ∆ t ˜ A ( p ( U ) x − f ( U )) x + O ( ξ ) ,u n +1 = u n + ∆ tb T A − ˜ A ( p ( U ) x − f ( U )) x + O ( ξ ) . (43)Assuming that the IMEX Runge-Kutta scheme is GSA, the term b T A − ˜ A = e Ts ˜ A = ˜ b T and, as ε →
0, the scheme relaxes to the explicit one, i.e., U − u n e ∆ t + ˜ Af ( U ) x = ˜ Ap ( U ) xx u n +1 − u n ∆ t + ˜ b T f ( U ) x = ˜ b T p ( U ) xx . (44) Generalization to the case α ∈ [0 , . In the case α ∈ [0 ,
1) analogous computations show that the characteristic speeds of the hyper-bolic part are given byΛ α ± (∆ t, ε ) = 12 (cid:32) cb T B ˜ Ae ± (cid:114)(cid:16) cb T B ˜ Ae (cid:17) + 4 ε − α ∆ t b T Be (cid:33) (45)where here ζ = ε α / ∆ t . We haveΛ α ± (∆ t,
0) = 12 ( c ± | c | ) , Λ α ± (0 , ε ) = ± ε α . A similar analysis performed in this case shows that U = u n e − ∆ t ˜ Af ( U ) x + ∆ tε (1 − α ) ˜ A ( p ( U ) xx + O ( ξ ) ,u n +1 = u n − ∆ t ˜ b T f ( U ) x + ∆ tε (1 − α ) ˜ b T p ( U ) xx + O ( ξ ) . As ε → U − u n e ∆ t + ˜ Af ( U ) x = 0 u n +1 − u n ∆ t + ˜ b T f ( U ) x = 0 . (46)All the previous results can be stated by the following Theorem 1
If the IMEX-RK scheme (30)-(31) applied to (1) satisfies the GSA property then,as ε → , it becomes the explicit RK scheme characterized by the pair ( ˜ A, ˜ w ) applied to thelimit convection-diffusion equation (6) for α = 1 and to the limit scalar conservation law (3) for α ∈ [0 , . Some remarks are in order.
Remark 1 The advantage of formulating the unified IMEX-RK approach in the form (33)-(34) isthat we can now adopt different space discretizations for the various term appearing inthe scheme. Typically, we use classical hyperbolic-type schemes, like WENO, for thespace derivatives characterizing the hyperbolic part (which now has finite characteristicspeeds) and centered discretization for the second order term characterizing the asymptoticparabolic behavior. • If the IMEX-RK scheme satisfies the GSA property the numerical solution is the same asthe last stage, then V s = v n +1 and, observing that ˜ b s = 0 , from the second equation of(34), we have v n +1 − v n ∆ t + 1 ε s − (cid:88) i =1 ˜ b i p ( U j ) x = − ε (cid:32) s − (cid:88) i =1 ( b i V i − ˜ b i f ( U i )) + b s v n +1 (cid:33) . Solving for v n +1 , after some algebra, the equation can be written as v n +1 − v n ∆ t + 1 ε + b s ∆ t s − (cid:88) i =1 ˜ b i p ( U j ) x = − ε + b s ∆ t (cid:32) s − (cid:88) i =1 ( b i V i − ˜ b i f ( U i )) + b s v n (cid:33) . Assuming p (cid:48) ( u ) = 1 and using the fact that U = u n e + O (∆ t ) we can now write thehyperbolic part in (33) as u t + ζb T Bev x + cb T B ˜ Aeu x = O (∆ t ) v t + 1 ε + a ss ∆ t u x = O (∆ t ) . (47) Therefore, the characteristic speeds read Λ ± (∆ t, ε ) = 12 cb T B ˜ Ae ± (cid:115)(cid:16) cb T B ˜ Ae (cid:17) + 4 ζb T Beε + a ss ∆ t , (48) and now when ε → , as in the first order case discussed in Section 3, we get Λ ± (∆ t,
0) = 12 ( c ± | c | ) . Although the final schemes developed in the previous section will work independently on ε and α , for small values of ε they relax to an explicit RK scheme originating a time step restrictionof the type ∆ t ≈ ∆ x /ε − α , for α ∈ (0 , ε , only the case α = 1 posesstability restriction. For this reason, we shall consider α = 1 in this subsection.The natural idea here is to treat the term p ( u ) x in (1) implicitly and to observe that inthe limit case, i.e. ε →
0, the IMEX-RK scheme relaxes to an IMEX-RK scheme for the limitconvection-diffusion equation (6) where the diffusion term is now evaluated implicitly. Comparedto similar schemes presented in [6, 8, 9] this new approach has the advantage that is not based ona penalization technique and therefore avoids the difficult problem of the optimal determinationof the penalization parameter (see [6]). 13he unified IMEX-RK scheme for system (1), now reads U = u n e − ∆ tAV x V = v n e + ∆ tε ˜ Af ( U ) − ∆ tε A ( V + p ( U ) x ) . (49)and u n +1 = u n − ∆ t b T V x v n +1 = v n + ∆ tε ˜ b T f ( U ) − ∆ tε b T ( V + p ( U ) x ) . (50)Now, solving the second equation in (49) for V , one obtains V = ( ζI + A ) − (cid:16) ζv n e + ˜ Af ( U ) − Ap ( U ) x (cid:17) , (51)where ζ = ε / ∆ t . Using this relation in the first equation of (49) we get for the internal stages U − u n ∆ t + ζA ( ζI + A ) − v nx e = A ( ζI + A ) − ( Ap ( U ) x − ˜ Af ( U )) x ,V − v n ∆ t + 1 ε ˜ Af ( U ) = − ε A ( V + p ( U ) x ) , (52)and similarly for the numerical solution u n +1 − u n ∆ t + ζb T ( ζI + A ) − v nx e = b T ( ζI + A ) − ( Ap ( U ) x − ˜ Af ( U )) x ,v n +1 − v n ∆ t + 1 ε ˜ b T f ( U ) = − ε b T ( V + p ( U ) x ) . (53)Then setting B = ( ζI + A ) − for small values of ∆ t , from (53), we get the system u t + ζb T Bev x + b T B ˜ Af ( U ) x = b T BAp ( U ) xx + O (∆ t ) v t + 1 ε b T p ( U ) x = − ε (cid:16) b T V − ˜ b T f ( U ) (cid:17) + O (∆ t ) , (54)which is similar to Eq. (35), except that now A and b T appear in front of p ( U ) terms in place of˜ A and ˜ b respectively. Therefore, we have the same characteristic speeds as in the unified IMEX-RK approach presented in the previous section and the same conclusions on the hyperbolic CFLcondition holds true.Concerning the AP property, in the limit ζ → V = − p ( U ) x + A − ˜ Af ( U ) , and by the GSA property, i.e. b T A − ˜ A = e Ts ˜ A = ˜ b T , we get from (53) v n +1 = − e Ts ( p ( U ) x − A − ˜ Af ( U )) . Thus, scheme (52)-(53), becomes an IMEX-RK scheme for the convection-diffusion equation(6) U − u n e ∆ t + ˜ Af ( U ) x = Ap ( U ) xx ,u n +1 − u n ∆ t + ˜ b T f ( U ) x = b T p ( U ) xx . (55)14learly, scheme (55) is a consistent approximation of the limit equation (6) where now thediffusion term is evaluated implicitly, therefore the CFL condition of such scheme is uniquelydetermined by the hyperbolic restriction ∆ t ∼ ∆ x .A similar analysis in the case α ∈ [0 ,
1) for scheme (49)-(50) under the GSA assumptionproduces the explicit Runge-Kutta scheme U = u n − ∆ t ˜ Af ( U ) x ,u n +1 = u n − ∆ t ˜ b T ˜ Af ( U ) x . (56)Therefore we can summarize the results in the following Theorem 2
If the IMEX-RK scheme (49)-(50), applied to (1) for α = 1 , satisfies the GSAproperty then, as ε → , it becomes the IMEX-RK method characterized by the pairs ( ˜ A, ˜ b ) and ( A, b ) for the limit convection-diffusion equation (6). Otherwise for α ∈ [0 , the IMEX-RKscheme as ε → yields the explict RK method characterized by the pair ( ˜ A, ˜ b ) for the limit scalarconservation law (4). In this section we present numerical results that confirm the validity of the new approachpresented in section 3.3. In all tests we used the last approach described in section 3.3, sothat we remove the parabolic restriction in the limit of small ε and α = 1. All the numericalexamples presented here refer to the IMEX-RK schemes reported in the Appendix. We shall usethe notation NAME( ν, σ, p ), where the triplet ( ν, σ, p ) where ν, σ and p represent respectivelythe number of explicit function evaluations, the number of implicit function evaluations and theorder of accuracy.In order to avoid spurious numerical oscillations arising near discontinuities of the solutions,we use interpolating non-oscillatory algorithms, like WENO method, [36]. In these numericaltest, as emphasized in Remark 1, we use classical hyperbolic-type schemes, like finite differencediscretization with WENO reconstruction, for the space derivatives characterized the hyperbolicpart while for the second order term p ( u ) xx we used the standard 2-th and 4-th order finitedifference technique. Note that when we consider third-order IMEX R-K scheme we use 4-thorder finite difference to discretize the term p ( u ) xx except at the nearby boundary points wherea 3-rd order formula was implemented, and this guarantees to achieve a global third-order.In all our numerical results we take ∆ t = λ CF L ∆ x in all regimes. The precise choice of λ CF L is reported in the figure captions.As a mathematical model for our numerical experiments we consider the Ruijgrook-Wumodel of the discrete kinetic theory of rarefied gases. The model describes a two-speed gas inone space dimension and corresponds to the system [17, 25, 35]
M ∂ t f + + ∂ x f + = − af + − bf − − cf + f − ) ,M ∂ t f − − ∂ x f − = 1Kn ( af + − bf − − cf + f − ) , (57)where f + and f − denote the particle density distribution at time t , position x and with velocity+1 and − M is the Mach number of the system15nd a , b and c are positive constants which characterize the microscopic interactions. The local(Maxwellian) equilibrium is defined by f + = bf − a − cf − . (58)The macroscopic variables for the model are the density u and momentum v defined by u = f + + f − , v = ( f + − f − ) /M. (59)The nondimensional multiscale problem is obtained taking M = ε α and Kn = ε , the Reynoldsnumber of the system is then defined as usual according to Re = M/Kn = 1 /ε − α . The model,as we will see, has the nice feature to provide nontrivial limit behaviors for several values of α including the corresponding compressible Euler ( α = 0) limit and the incompressible Euler( α ∈ (0 , α = 1) limits. Test 1. Diffusive scaling in the linear case
For this numerical test we consider the case α = 1 with c = 0, a = 1 + Aε and b = 1 − Aε in ther.h.s. of (57).Adding and subtracting the two equations in (57) one obtains the following macroscopicequations for u and v ∂ t u + ∂ x v = 0 ,∂ t v + 1 ε ∂ x u = − ε ( v − Au ) . (60)In the limit ε → v = Au − ∂u∂x and substituting in the fist equation this gives the limiting advection-diffusion equation u t + Au x = u xx . Next we fix A = 1 and observe that the limiting advection-diffusion equation admits the exactsolution u ( x, t ) = e − t sin( x − t ) , v ( x, t ) = e − t (sin( x − t ) − cos( x − t )) , (61)on the domain [ − π, π ] with periodic boundary conditions. We start to consider as initial con-ditions (61) at t = 0, and choose ε = 10 − with final time T = 0 .
1. The numerical results arecompared with (61) at t = T . Relative errors and convergence orders are reported in Tables 1-3.In order to check the temporal order of convergence of these schemes, ∆ x decreases with thetime step ∆ t accordingly to the CFL condition ∆ t = 0 . x . In the Tables we show the orderof convergence as p = log ( || E ∆ t || ∞ / || E ∆ t/ || ∞ ) , with E ∆ t the relative error computed with time step ∆ t . We consider the error obtained with N = 40 time steps up to N = 640 time steps in the interval [0 , T ]. We observe that the classicalorder of the methods is maintained in the limit case for the density u . We omit the convergence16 , −3 −2 −1 0 1 2 3−1.5−1−0.500.511.5 Figure 1: Test 1. Comparison between classical ARS(2,2,2) scheme (left) and the BPR(4,4,2)scheme (right) with N = 40 and ∆ t = ∆ x .results of ARS(2,2,2) scheme (73) since they are essentially the same as the ones obtained withCK(2,2,2) scheme reported in (74).Typically, classical schemes, as ARS(2,2,2), CK(2,2,2) and BPR(3,4,3), which satisfy onlythe GSA property, have few numbers of internal stages but maintain, in the limit ε →
0, theorder in time only for the u -component, while for the v -component they do not guarantee eventhe consistency [8]. As an example, we report the convergence table for the variable v in Table 2for the CK(2,2,2) scheme (74) and we observe that we do not obtain the correct classical orderof accuracy for v in the limit case ε → v , is that this classicalschemes does not satisfy the additional order conditions (71) (see Appendix) that guaranteesthe consistency of the scheme and the order up to 2 for the variable v . As a comparison weconstruct a new scheme, BPR(4,4,2), that satisfies the additional order conditions (71), and inTable 3 we observe the correct convergence rate for both components, u and v . Furthermore,in Fig. 1, we also compare the numerical solutions (star points) for u and v obtained with thesecond order ARS(2,2,2) scheme (73) and BPR(4,4,2) scheme. In the same figure, the exactsolution is plotted by a continuous line. Note again that the new scheme BPR(4,4,2) shows thecorrect behaviour for the v variable.Of course there are some advantage and disadvantages to consider additional order conditions(71). The advantage is that they guarantee the correct order in the limit case ( ε →
0) for bothvariables. On the other hand, to construct schemes that satisfy (71), requires a larger numberof internal stages, due to these extra order conditions.Next, we consider a Riemann problem with initial data (cid:26) u L = 4 . , v L = 0 , − < x < ,u R = 2 . , v R = 0 , < x < , (62)and inflow and outflow boundary conditions. The exact solution for the limit advection-diffusion17able 1: Test 1. Converge rates for the density u with ε = 10 − .Method N L ∞ error OrderARS(1,1,1) 40 6 . e − −− ARS(1,1,1) 80 3 . e −
03 0 . . e −
03 0 . . e −
04 0 . . e −
04 0 . . e − −− CK(2,2,2) 80 3 . e −
05 1 . . e −
05 1 . . e −
06 2 . . e −
07 1 . . e − −− BPR(3,4,3) 80 7 . e −
07 2 . . e −
07 2 . . e −
08 2 . . e −
09 2 . v for the scheme CK(2,2,2) with ε = 10 − .Method N L ∞ v -error Order of v CK(2,2,2) 40 1 . e − −− CK(2,2,2) 80 1 . e − − . . e −
03 3 . . e − − . . e − − . u ( x, t ) = 12 ( u L + u R ) + 12 ( u L − u R )erf (cid:18) t − x √ t (cid:19) , where erf( x ) denotes the error function. We take ∆ x = 0 . T = 3 .
0. In Figure 2 wecompare the numerical solutions computed by schemes SP(1,1,1), BPR(2,4,4) and BPR(3,3,5)for the mass density in the intermediate regime ( ε = 0 .
5) and in diffusive one ( ε = 10 − ) witha reference solution obtained with ∆ x = 0 . t = 0 . x . Test 2. Multiscale limit in the nonlinear case
Here we consider the nonlinear Ruijgrok-Wu model Eq.(57) for α ∈ [0 ,
1] and interaction pa-rameters c = 2 ε and a = b = 1.Adding and subtracting the two equations in (57) one obtains the following macroscopic18able 3: Test 1. Converge rates for u and v for the scheme (75) with ε = 10 − .Method N L ∞ u -error Order u L ∞ v -error Order v BPR(4,4,2) 40 1 . e − −− . e − −− BPR(4,4,2) 80 4 . e −
05 1 . . e −
05 1 . . e −
05 1 . . e −
05 1 . . e −
06 2 . . e −
06 1 . . e −
07 1 . . e −
06 1 . -10 -8 -6 -4 -2 0 2 4 6 8 1022.22.42.62.833.23.43.63.84 reference solutionBPR(4,4,2)ARS(1,1,1)BPR(3,4,3) -10 -8 -6 -4 -2 0 2 4 6 8 1022.22.42.62.833.23.43.63.84 ARS(1,1,1)Reference SolutionBPR(4,4,2)BPR(3,4,3)
Figure 2: Test 1. Solution of problem (60) with initial data (62), ∆ x = 0 . t = 0 . x forthe density u . Left: the rarefied regime ε = 0 .
5. Right: the parabolic regime ε = 10 − .equations for u and v ∂ t u + ∂ x v = 0 ,∂ t v + 1 ε α ∂ x u = 1 ε α (cid:26) − v + 12 (cid:0) u − ε v (cid:1)(cid:27) . (63)For small values of ε the model behaviour can be derived by the Chapman-Enskog expansionand is characterized by the viscous Burgers equation v = 12 u − ε − α ∂ x u + ε α u ∂ x u + O ( ε ) ∂ t u + ∂ x (cid:18) u (cid:19) = ε α ∂ x (cid:20)(cid:18) ε α − u (cid:19) ∂ x u (cid:21) + O ( ε ) . (64)We consider two different initial conditions. The first one is given by two local Maxwelliancharacterized by (cid:26) u L = 1 . , v L = 0 , − < x < ,u R = 2 . , v R = 0 , < x < , (65)with v = [(1 + u ε ) / − /ε . 19e show in Figure 3 the numerical solution for u in the case α = 1 using BPR(4,4,2) andBPR(3,4,3) schemes in the rarefied ( ε = 0 .
4) and parabolic ( ε = 10 − ) regimes with ∆ x = 0 . T = 0 .
2. The solution of both schemes is in very good agreement with thereference solution computed with ∆ x = 0 . -10 -8 -6 -4 -2 0 2 4 6 8 1011.11.21.31.41.51.61.71.81.92 Reference SolutionBPR(3,4,3)BPR(4,4,2) , -10 -8 -6 -4 -2 0 2 4 6 8 1011.11.21.31.41.51.61.71.81.92 BPR(4,4,2)BPR(3,4,3)Reference Solution
Figure 3: Test 2. Numerical solutions of problem (63) for α = 1 with initial conditions (65) at T = 2 . t = 0 . x . Left: the rarefied regime for u with ε = 0 .
4. Right: the parabolicregime for u with ε = 10 − .The last test case that we consider is the propagation of an initial square wave. The initialprofile is specified as u = 1 . , v = 0 . | x | < . , u = v = 0 for | x | > .
125 (66)with reflecting boundary conditions, i.e. v = 0, u x = 0 on the boundary.We integrate the equations over [ − . , .
5] with 200 spatial cells. In Figure 4 we plot thebehavior of the system in the rarefied regime for ε = 0 . α = 0 at time T = 0 . ε = 10 − ) with α = 0 . , .
75. We use BPR(3,4,3) scheme (76). The numericalsolutions for the mass density u and momentum v are computed in the rarefied regime and inthe diffusive regime with ∆ t = 0 . x and are depicted with a reference solution obtained usingfine grids with ∆ x = 0 . α gives the corresponding Reynolds numbers of the problem, i.e., Re = 10000 and Re = 100. The output of the solution in the parabolic regime is given at T = 0 . t = 0 . t = 0 . x ) and a right moving shock is formed. Although ε is underresolved the schemeproposed captures well the correct behavior of the equilibrium equation independently of α . Test 3. Multiscale space varying limit in the nonlinear case
Finally, we consider the multi scale relaxation system of the previous section, in the case wherethe parameter α = α ( x ) ∈ [0 ,
1] depends on the space variable. Therefore, the limit behavior ofthe system may depend on the particular region of the computational domain.Now we apply our schemes to system (63) considering two different cases of the varying α number. The domain is chosen to be x ∈ [ − . , . α increases smoothly from20 , −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5−0.6−0.4−0.200.20.40.60.8 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5−0.200.20.40.60.81 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5−0.100.10.20.30.40.50.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.500.10.20.30.40.50.60.70.8 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5−0.0500.050.10.150.20.250.3 Figure 4: Test 2. Numerical solutions of system (63) with initial conditions (66) for the massdensity u (left) and the momentum v (right) in the rarefied regime (top panels) with ε = 0 . α = 0 and ∆ t = 0 . x = 0 .
005 at T = 0 .
2. In the parabolic regime ε = 10 − , with α = 0 . α = 0 .
75 (bottom panels) with ∆ t = 0 .
004 (i.e. ∆ t = 0 . x ) at time T = 0 . α to O (1), by the formula, α ( x ) = α + 0 . x + 0 . α = 10 − . In the second case we consider the function α which contains a discontinuity (cid:26) α L = 0 . , − . < x < ,α R = 1 . , < x < . . (68)The two different scenarios for α are depicted in the Figure 5. -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.500.10.20.30.40.50.60.70.80.91 parabolic regimehyperbolic regime , -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.500.10.20.30.40.50.60.70.80.91 parabolic regimehyperbolic regime Figure 5: Test 3. Space varying α ( x ). Left: the smooth profile (67). Right: the discontinouosprofile (68)The initial profile is specified by (66) with reflecting boundary conditions and we fix ε = 10 − and ∆ t = 0 . x . Here we use BPR(3,3,5) scheme. The reference solutions are computed usinga fine grid with ∆ x = 0 . T = 0 .
05 (left) and T = 0 . -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.500.10.20.30.40.50.60.70.80.91 reference solution α (x)200 points100 points -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.500.10.20.30.40.50.60.70.80.91 refernce solution α (x)100 points200 points Figure 6: Test 3. For ε = 10 − , the left and right pictures show the behaviour of the component u in the two different cases of α . 22 Conclusions
In this work we have developed a new IMEX-RK approach that is capable correctly to capturethe asymptotic behavior of hyperbolic balance laws with relaxation under different kinds ofscaling. Previous IMEX-RK schemes were designed to deal specifically with one kind of scaling,either hyperbolic [7, 10, 22, 32] or parabolic [6, 21, 25, 29], whereas the present schemes arerobust enough to be able to deal with both scalings. Related approaches were presented in[8, 24, 27, 30].From a physical viewpoint, these scaling limits corresponds to the classical fluid-dynamicscalings in the kinetic theory of rarefied gases that lead from the Boltzmann equation to itscompressible and incompressible limits. Several numerical test for a simple kinetic model whichpossesses different asymptotic limits have confirmed the validity of the present approach. In thenear future, we hope to extend this class of IMEX-RK schemes to the full Boltzmann equationby adopting the penalization techniques developed in [15, 16].
Acknowledgments
The authors are grateful to Prof. M. Lemou and Prof. G. Dimarco for helpful observations.The research that led to the present paper was partially supported by the research grant
Nu-merical methods for uncertainty quantification in hyperbolic and kinetic equations of the groupGNCS of INdAM, by ITN-ETN Marie-Curie Horizon 2020 program ModCompShock,
Modelingand computation of shocks and interfaces , Project ID: 642768 and by the INDAM-GNCS 2017research grant
Numerical methods for hyperbolic and kinetic equation and applications . A Appendix: IMEX-RK schemes
In order to achieve higher then first order accuracy in time for (GSA) IMEX-RK schemes, weneed to observe first that there are no second order GSA IMEX-RK methods of type I with s = 3, i.e. with three stages [6]. Thus, to construct a second order GSA IMEX Runge-Kuttaschemes, we consider the following proposition [6] Proposition 1
The only type of second order GSA IMEX-RK scheme with three levels s = 3 that satisfies the classical second order conditions (26) is the type II where c i = ˜ c i for all i = 2 , ..., s − , with c = 0 and c s = 1 . As usual, order conditions are obtained by matching the Taylor expansion of the exactsolution and the numerical one, up to terms of the prescribed order. In the relaxed case, i.e. ε →
0, system (1) for α = 1 reduces to v = − p ( u ) x + f ( u ) , u t + f ( u ) x = p ( u ) xx . (69)The unified IMEX-RK approach described in Section 3.3 provides a scheme that converges toan explicit-implicit scheme for the limit convection-diffusion equation in (69). Performing thesame analysis as proposed in [8] for system (1) we obtain for the u -component the classical orderconditions, while some additional order conditions for the v -component are required in orderto have consistency and maintain the classical order in the limit case. We recall that the GSAassumption of the method guarantees that IMEX-RK scheme relaxes at the same IMEX-RKone when ε →
0, but in order to maintain the order of accuracy of the scheme in the limit wemust impose some additional order conditions.23elow we list these new additional order conditions that the v -component must satisfy upto second orderconsistency b T A − ˜ Ae = 1 , first order b T A − ˜ Ac = 1 , b T A − ˜ A ˜ c = 1 , second order b T A − ˜ Ac = 1 , b T A − ˜ A ˜ c = 1 , b T A − ˜ A ˜ cc = 1 , b T A − ˜ Ac ˜ c = 1 ,b T A − ˜ AAc = 1 / , b T A − ˜ AA ˜ c = 1 / ,b T A − ˜ A ˜ Ac = 1 / , b T A − ˜ A ˜ A ˜ c = 1 / . (70)Note that in order to reduce the number of the additional order conditions we require that ˜ c = c ,this assumption simplifies a lot the number of coupling conditionsconsistency b T A − ˜ Ae = 1 , first order b T A − ˜ Ac = 1 , second order b T A − ˜ Ac = 1 , b T A − ˜ AAc = 1 / , b T A − ˜ A ˜ Ac = 1 / . (71)Finally, we present the different IMEX-RK schemes, up to order three, used in our numericaltests. Below, these schemes are represented as usual by the double Butcher tableau. On theleft we have the explicit part and on the right the implicit part of the IMEX schemes. All theschemes satisfy the GSA property, but only the new scheme BPR(4 , ,
2) satisfies the additionalorder conditions (71). • First order ARS(1 , ,
1) scheme 0 0 01 1 01 0 0 0 01 0 10 1 (72) • Second order ARS(2 , ,
2) scheme [1]0 0 0 0 γ γ δ − δ δ − δ γ γ
01 0 1 − γ γ − γ γ . (73)where γ = 1 − / √ δ = 1 − / γ , • Second order CK(2 , ,
2) scheme [11]0 0 0 02 / / / / / / / − / √ / − √ / / − √ / − / √ / − √ / / − √ / − / √ / − √ / . (74) • Second order BPR(4 , ,
2) scheme 24 0 0 0 0 01 / / / / − / / / / / / / / / / / / / / /
24 11 /
24 1 / /
24 1 / / /
40 11 /
24 1 / / / . (75) • Third order BPR(3 , ,
3) scheme [6]0 0 0 0 0 01 1 0 0 0 02 / / / / / / / / / / / / / − / / / / / / − / / / / − / / . (76) References [1]
U. M. Ascher, S. J. Ruuth, R. J. Spiteri , Implicit-Explicit Runge-Kutta Methods forTime-Dependent Partial Differential Equations , Appl. Numer. Math, 1997, V. 25, pages151–167.[2]
S. Boscarino , Error analysis of IMEX Runge-Kutta methods derived from differential-algebraic systems , SIAM Journal on Numerical Analysis, Vol. 45, No. 4, pp. 1600–1621,2007.[3]
S. Boscarino , On an accurate third order implicit-explicit RungeKutta method for stiffproblems , Applied Numerical Mathematics, Vol. 59, pp. 1515–1528, 2009.[4]
C. Bardos, C.D. Levermore, F. Golse , Fluid dynamic limit of kinetic equations II:Convergence proofs for the Boltzmann equations , Comm. Pure Appl. Math., 46, (1993),667-753.[5]
S. Boscarino, L. Pareschi , On the asymptotic properties of IMEX Runge–Kutta schemesfor hyperbolic balance laws , J. Comput. Appl. Math. (to appear).[6]
S. Boscarino, L. Pareschi, G. Russo , Implicit-Explicit Runge-Kutta schemes for hy-perbolic systems and kinetic equations in the diffusion limit , SIAM J. Sci. Comput., Vol.35, No. 1, pp. A22A51.[7]
S. Boscarino, G. Russo , On a class of uniformly accurate IMEX Runge-Kutta schemesand applications to hyperbolic systems with relaxation , SIAM J. Sci. Comput. Vol. 31, No.3, pp. 1926-1945.[8]
S. Boscarino, G. Russo , Flux-Explicit IMEX Runge-Kutta schemes for hyperbolic toparabolic relaxation problems , SIAM J. NUMER. ANAL., Vol. 51, No. 1, 163–190.[9]
S. Boscarino, P. Le Floch, G. Russo , High-Order asymptotic-preserving methods forfully non linear relaxation problems , SIAM J. Sci. Comput., Vol. 36, No. 2, pp. A377–A395.2510]
R.E. Caflisch, S. Jin, G. Russo , Uniformly accurate schemes for hyperbolic systemswith relaxation , SIAM J. Num. Anal., 34, 1, (1997), 246-281.[11]
M.H. Carpenter, C.A. Kennedy , Additive Runge-Kutta schemes for convection-diffusion-reaction equations
Appl. Numer. Math. 44 (2003), no. 1-2, 139–181.[12]
C. Cercignani, R. Illner, M. Pulvirenti , The mathematical theory of dilute gases ,Applied Mathematical Sciences, 106, Springer-Verlag, New York, (1994).[13]
G.Q. Chen, C.D. Levermore, T.P. Liu , Hyperbolic conservation laws with stiff relax-ation terms and entropy , Comm. Pure and App. Math., XLVII, (1994), 787-830.[14]
C. Cercignani , The Boltzmann equation and its applications , Springer Varlag New York,(1988).[15]
G. Dimarco, L. Pareschi , Asymptotic-preserving IMEX Runge-Kutta methods for non-linear kinetic equations , SIAM J. Num. Anal. (2013), 1064–1087.[16]
F. Filbet, S. Jin , A class of asymptotic preserving schemes for kinetic equations andrelated problems with stiff sources , J. Comp. Phys., 229 (2010), pp. 7625-7648.[17]
E. Gabetta, B. Perthame , Scaling limits of the Ruijgrook-Wu model of the Boltzmannequation , Proceedings of the International Conference on Nonlinear Equations and Appli-cations, Bangalore 19-23 August 1996. Springer-Verlag.[18]
E. Gabetta, L. Pareschi, G. Toscani , Relaxation schemes for nonlinear kinetic equa-tions , SIAM J. Numer. Anal. 34, (1997), no. 6, 2168–2194.[19]
E. Hairer, G. Wanner , Solving Ordinary Differential Equation II: stiff and DifferentialAlgebraic Problems. Springer Series in Comput. Mathematics, Vol. 14, Springer-Verlag1991, Second revised edition 1996.[20]
E. Hairer, S.P. Norsett, G. Wanner , Solving Ordinary Differential Equation I: Nons-tiff Problems. Springer Series in Comput. Mathematics, Vol. 8, Springer-Verlag 1987, Secondrevised edition 1993.[21]
A.Klar , An asymptotic-induced scheme for nonstationary transport equations in the dif-fusive limit , SIAM J. Numer. Anal. 35, no. 3, (1998), 1073–1094.[22]
S. Jin , Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms ,J. Comput. Phys. 122, (1995), 51–65.[23]
S. Jin, C.D. Levermore , Numerical schemes for hyperbolic conservation laws with stiffrelaxation terms , J. Comput. Phys., 126 (1996), 449–467.[24]
S. Jin, L. Pareschi , Asymptotic-preserving (AP) schemes for multiscale kinetic equa-tions: a unified approach , Hyperbolic Problems: Theory, Numerics, Applications, Vol. 141,International Series of Numerical Mathematics, (2001), 573–582.[25]
S. Jin, L. Pareschi, G. Toscani , Diffusive relaxation schemes for discrete-velocity kineticequations , SIAM J. Numer. Anal. 35, (1998), 2405–2439.[26]
S. Jin, L. Pareschi, G. Toscani , Uniformly accurate diffusive relaxation schemes fortransport equations , SIAM J. Numer. Anal. 38, (2000), 913–936.2627]
M. Lemou, L. Mieussens , A new asymptotic preserving scheme based on micro-macroformulation for linear kinetic equations in the diffusion limit , SIAM J. Sci. Comput., 31(2008), 334–368.[28]
T.P. Liu , Hyperbolic conservation laws with relaxation , Comm. Math. Phys., 108, 153(1987).[29]
G. Naldi, L. Pareschi , Numerical schemes for kinetic equations in diffusive regimes ,Appl. Math. Letters, 11, 2, (1998), 29–35.[30]
G. Naldi, L. Pareschi , Numerical schemes for hyperbolic systems of conservation lawswith stiff diffusive relaxation . SIAM J. Numer. Anal. 37 (2000), no. 4, 1246–1270.[31]
L. Pareschi , Characteristic-based numerical schemes for hyperbolic systems with nonlinearrelaxation , Rend. Circ. Matem. di Palermo, II, 57, (1998), 375–380.[32]
L. Pareschi and G. Russo , Implicit-Explicit Runge-Kutta schemes and applications tohyperbolic systems with relaxations , J. Sci. Comput., 25 (2005), 129–155[33]
R.B. Pember , Numerical methods for hyperbolic conservation laws with stiff relaxation, I.Spurious solutions , SIAM J. Appl. Math., 53, (1993), 1293–1330.[34]
P.L. Roe, M. Arora , Characteristic-based schemes for dispersive waves I. The method ofcharacteristic for smooth solutions , Numerical Meth. for PDE’s, 9, (1993), 459–505.[35]
W. Ruijgrook, T.T. Wu , A completely solvable model of the nonlinear Boltzmann equa-tion , Physica, 113 A, (1982), 401–416.[36]