Implicit-Explicit multistep methods for hyperbolic systems with multiscale relaxation
aa r X i v : . [ m a t h . NA ] J a n IMPLICIT-EXPLICIT MULTISTEP METHODS FOR HYPERBOLIC SYSTEMSWITH MULTISCALE RELAXATION
GIACOMO ALBI, GIACOMO DIMARCO AND LORENZO PARESCHI
Abstract.
We consider the development of high order space and time numerical methods based onImplicit-Explicit (IMEX) multistep time integrators for hyperbolic systems with relaxation. More specif-ically, we consider hyperbolic balance laws in which the convection and the source term may have verydifferent time and space scales. As a consequence the nature of the asymptotic limit changes completely,passing from a hyperbolic to a parabolic system. From the computational point of view, standard numer-ical methods designed for the fluid-dynamic scaling of hyperbolic systems with relaxation present severaldrawbacks and typically lose efficiency in describing the parabolic limit regime. In this work, in the con-text of Implicit-Explicit linear multistep methods we construct high order space-time discretizationswhich are able to handle all the different scales and to capture the correct asymptotic behavior, inde-pendently from its nature, without time step restrictions imposed by the fast scales. Several numericalexamples confirm the theoretical analysis.
Key words.
Implicit-explicit methods, linear multistep methods, hyperbolic balance laws, fluid-dynamiclimit, diffusion limit, asymptotic-preserving schemes.
AMS subject classification.
Contents
1. Introduction 22. First order IMEX discretization 33. AP-explicit and AP-implicit IMEX-LM methods 53.1. AP-explicit methods in the diffusive case: α = 1 63.2. AP-explicit methods in the general case: α ∈ [0 , The work of G. Albi has been partially supported by the project Computer Science for Industry 4.0, ’MIUR Depart-ments of Excellence 2018-2022’. The authors acknowledge the partial support of GNCS-INdAM 2019 project Numericalapproximation of hyperbolic problems and applications.Computer Science Department, University of Verona, ( [email protected] ).Mathematics Department, University of Ferrara, Italy ( [email protected] ).Mathematics Department, University of Ferrara, Italy ( [email protected] ). Introduction
The goal of the present work is to develop high order numerical methods based on IMEX linearmultistep (IMEX-LM) schemes for hyperbolic systems with relaxation [14, 38, 42]. These systems oftencontain multiple space-time scales which may differ by several orders of magnitude. In fact, the variousparameters characterizing the models permit to describe different physical situations, like flows which passfrom compressible to incompressible regimes or flows which range from rarefied to dense states. This isthe case, for example, of kinetic equations close to the hydrodynamic limits [5,13,15,49]. In such regimesthese systems can be more conveniently described in terms of macroscopic equations since these reducedsystems permit to describe all the features related to the space-time scale under consideration [5, 49].However, such macroscopic models can not handle all the possible regimes one is frequently interested in.For such reason one has to resort to the full kinetic models. They permit to characterize a richer physicsbut on the other hand they are computationally more expensive and limited by the stiffness induced bythe scaling under consideration [16].The prototype system we will use in the rest of the paper is the following [8, 41](1.1) ∂ t u + ∂ x v = 0 ,∂ t v + 1 ε α ∂ x p ( u ) = − ε α ( v − f ( u )) , α ∈ [0 , ε is the scaling factor and α characterizes the different type of asymptotic limit that can beobtained. The condition p ′ ( u ) > ± p p ′ ( u ) /ε α . Note that, except for the case α = 0, the eigenvalues are unboundedfor small values of ε .System (1.1) is obtained from a classical (2 × p -system with relaxation under the space-time scaling t → t/ε α , x → x/ε and by the change of variables v = ˜ v/ε α , f ( u ) = ˜ f ( u ) /ε α , where ˜ v is the originalunknown and ˜ f ( u ) the original flux associated to the variable u in the non rescaled p -system. For α = 0,system (1.1) reduces to the usual hyperbolic scaling(1.2) ∂ t u + ∂ x v = 0 ,∂ t v + ∂ x p ( u ) = − ε ( v − f ( u )) , whereas for α = 1 yields the so-called diffusive scaling(1.3) ∂ t u + ∂ x v = 0 ,∂ t v + 1 ε ∂ x p ( u ) = − ε ( v − f ( u )) . More in general, thanks to the Chapman-Enskog expansion [15], for small values of ε we get from (1.1)the following nonlinear convection-diffusion equation(1.4) v = f ( u ) − ε − α ∂ x p ( u ) + ε α f ′ ( u ) ∂ x u + O ( ε ) ,∂ t u + ∂ x f ( u ) = ε α ∂ x " (cid:18) p ′ ( u ) ε α − f ′ ( u ) (cid:19) ∂ x u + O ( ε ) . In the limit ε →
0, for α ∈ [0 , v = f ( u ) ,∂ t u + ∂ x f ( u ) = 0 , while, when α = 1, in the asymptotic limit we obtain the following advection-diffusion equation(1.6) v = f ( u ) − ∂ x p ( u ) ,∂ t u + ∂ x f ( u ) = ∂ xx p ( u ) . Note that, the main stability condition for system (1.4) corresponds to(1.7) f ′ ( u ) < p ′ ( u ) ε α , MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 3 and it is always satisfied in the limit ε → α >
0, whereas for α = 0 the function p ( u ) and f ( u )must satisfy the classical sub-characteristic condition [14, 38].The space-time scaling just discussed, in classical kinetic theory, is related to the hydrodynamical limitsof the Boltzmann equation. In particular, for α = 0 it corresponds to the compressible Euler scaling,whereas for α ∈ (0 ,
1) to the incompressible Euler limit. In the case α = 1 dissipative effects becomenon-negligible and we get the incompressible Navier-Stokes scaling. We refer to [13, Chapter 11] and [49]for further details and the mathematical theory behind the hydrodynamical limits of the Boltzmannequation. Moreover, we refer to [36,39] for theoretical results on the diffusion limit of a system like (1.3).The development of numerical methods to solve hyperbolic systems with stiff source terms has attractedmany researches in the recent past [10, 11, 17, 22, 24, 24, 28, 31, 33, 34, 37, 43, 44]. The main computationalchallenge is related to the presence of the different scales that require a special care to avoid loss ofstability and spurious numerical solutions. In particular, in diffusive regimes the schemes should becapable to deal with the very large characteristic speeds of the system by avoiding a CFL condition ofthe type ∆ t = O ( ε α ). A particular successful class of schemes is represented by the so-called asymptotic-preserving (AP) schemes which aims at preserving the correct asymptotic behavior of the system withoutany loss of efficiency due to time step restrictions related to the small scales [18, 20, 30, 32].In the large majority of these works, the authors focused on the specific case α = 0, where a hyperbolicto hyperbolic scaling is studied, or to the case α = 1, where a hyperbolic to parabolic scaling is analysed.Very few papers have addressed the challenging multiscale general problem for the various possible valuesof α ∈ [0 ,
1] (see [8, 22, 32, 40]) and all of them refer to one-step IMEX methods in a Runge-Kuttasetting. We refer to [1, 2, 4, 19, 21, 26, 45, 47] for various IMEX-LM methods developed in the literatureand we mention that comparison between IMEX Runge-Kutta methods and IMEX-LM methods havebeen presented in [26].In the present work, following the approach recently introduced by Boscarino, Russo and Pareschiin [8], we analyze the construction of IMEX-LM for such problems that work uniformly independently ofthe choices of α and ε . By this, we mean that the schemes are designed in such a way as to be stable for alldifferent ranges of the scaling parameters independently of the time step. At the same time, they shouldensure high order in space and time and should be able to accurately describe the various asymptoticlimits. Moreover, whenever possible, the above described properties must be achieved without the needof an iterative solver for non linear equations.In [8], using a convenient partitioning of the original problem, the authors developed IMEX Runge-Kutta schemes for a system like (1.1) which work uniformly with respect to the scaling parameters. Here,we extend these results to IMEX-LM, previous results for IMEX-LM methods refer to the case α = 0(see [19, 26]). Among others, there are two main reasons to consider the development of such schemes.First, in contrast to the IMEX Runge-Kutta case, for IMEX-LM it is relatively easy to construct schemesup to fifth order in time and typically they show a more uniform behavior of the error with respectto the scaling parameters. Second, thanks to the use of BDF methods, it is possible consider only oneevaluation of the source term per time step independently of the scheme order. The latter feature isparticularly significant in term of computational efficiency for kinetic equations where often the sourceterm represents the most expensive part of the computation [18, 19].The rest of the paper is organized as follows. In Section 2, we discuss the discretization of thesemultiscale problems and motivate our partitioning choice of the system by analyzing a simple first or-der IMEX scheme. Next, in Section 3, we introduce the general IMEX-LM methods and discuss theasymptotic-preserving properties of our approach. Two classes of schemes are considered, AP-explicitand AP-implicit accordingly to the way the diffusion term in the limit equation is treated. In Section 4,we perform a linear stability analysis for 2 × First order IMEX discretization
In this part, we discuss a first order IMEX time-discretization of the relaxation system (1.1) and weanalyze its relationship with a reformulated system in which the eigenvalues are bounded for any value ofthe scaling parameter ε . To this aim, following [8], we consider the following implicit-explicit first order GIACOMO ALBI, GIACOMO DIMARCO, AND LORENZO PARESCHI partitioning of system (1.1)(2.1) u n +1 − u n ∆ t = − ∂ x v n +1 ,ε α v n +1 − v n ∆ t = − (cid:0) ε − α ∂ x p ( u n ) + v n +1 − f ( u n ) (cid:1) . One can notice that in system (2.1) besides its implicit form, the second equation can be solved explicitlyby inversion of the linear term v n +1 . This gives(2.2) v n +1 = ε α ε α + ∆ t v n − ∆ tε α + ∆ t (cid:0) ε − α ∂ x p ( u n ) − f ( u n ) (cid:1) . Then, making use of the above relation and inserting it in the first equation, one gets(2.3) u n +1 − u n ∆ t + ε α ε α + ∆ t v nx + ∆ tε α + ∆ t ∂ x f ( u n ) = ∆ t ε − α ε α + ∆ t ∂ xx p ( u n ) , while a simple rewriting of the second equation gives(2.4) v n +1 − v n ∆ t + ε − α ε α + ∆ t ∂ x p ( u n ) = − ε α + ∆ t ( v n − f ( u n )) . Therefore, the IMEX scheme can be recast in an equivalent fully explicit form. Similarly to the continuouscase, depending on the choice of α , as ε →
0, we have different limit behaviors. For α ∈ [0 ,
1) we obtain(2.5) u n +1 − u n ∆ t + ∂ x f ( u n ) = 0 , whereas in the case α = 1 we get(2.6) u n +1 − u n ∆ t + ∂ x f ( u n ) = ∂ xx p ( u n ) . For small values of ∆ t , the scheme (2.3)-(2.4) corresponds up to first order in time to the system(2.7) ∂ t u + ε α ε α + ∆ t ∂ x v + ∆ tε α + ∆ t ∂ x f ( u ) = ∆ t ε − α ε α + ∆ t ∂ xx p ( u ) ,∂ t v + ε − α ε α + ∆ t ∂ x p ( u ) = − ε α + ∆ t ( v − f ( u )) , where the following Taylor expansion has been employed at t = t n = n ∆ tu n +1 − u n ∆ t = ∂ t u (cid:12)(cid:12) t = t n + O (∆ t ) , v n +1 − v n ∆ t = ∂ t v (cid:12)(cid:12) t = t n + O (∆ t ) . The main feature of system (2.7) is that its left-hand side has bounded characteristic speeds. These aregiven by(2.8) λ α ± (∆ t, ε ) = 12 (cid:16) γ (1 − θ α ) ± p γ (1 − θ α ) + 4 ε − α θ α (cid:17) , with θ α (∆ t, ε ) := ε α ε α + ∆ t , and where, for simplicity, we considered f ′ ( u ) = γ , γ ∈ R and p ′ ( u ) = 1 so that(2.9) ∂ x f ( u ) = f ′ ( u ) ∂ x u = γ∂ x u, ∂ x p ( u ) = p ′ ( u ) ∂ x u = ∂ x u. If we fix ε and send ∆ t → λ α ± (0 , ε ) = ± ε α , while for a fixed ∆ t , the characteristic speeds λ α + and λ α − are respectively decreasing and increasingfunctions of ε and, as ε →
0, they converge to(2.10) λ α ± (∆ t,
0) = 12 ( γ ± | γ | ) . In Figure 1, we show the shape of the eigenvalues (2.8) for γ = 1 and different values of the scalingparameter α and the time step ∆ t . We observe that when ε grows the absolute value of the eigenvalues MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 5
Figure 1.
Eigenvalues of the modified system (2.7) as a function of ε for different valuesof the time step ∆ t and choices of α .diminish accordingly and when ε diminishes the eigenvalues grow but remain bounded by the finite timestep.Thus, for a given ∆ t , if we denote by ∆ x the space discretization parameter, from the left hand sideof (2.7) we expect the hyperbolic CFL condition ∆ t ≤ ∆ x/ | γ | in the limit ε →
0. On the other hand,the stability restriction coming from the parabolic term requires ∆ t = O (∆ x ) when α = 1.In the next Section, we will show how to generalize the above arguments to the case of high orderIMEX multistep methods.3. AP-explicit and AP-implicit IMEX-LM methods
In this part, we focus our attention on order s , s -step IMEX-LM methods with s ≥ ε →
0. For clarity of presentation, we separate the discussion of the diffusive case α = 1 fromthe general case α ∈ [0 , GIACOMO ALBI, GIACOMO DIMARCO, AND LORENZO PARESCHI
AP-explicit methods in the diffusive case: α = 1 . In this case, we can write the s − stepIMEX-LM for the original hyperbolic system (1.1) as follows(3.1) u n +1 = − s − X j =0 a j u n − j − ∆ t s − X j = − c j ∂ x v n − j ,v n +1 = − s − X j =0 a j v n − j − ∆ tε s − X j =0 b j ∂ x p ( u n − j ) + s − X j = − c j v n − j − s − X j =0 b j f ( u n − j ) , where we introduced the following coefficients(3.2) Explicit a T = ( a , a , . . . , a s − ) b T = ( b , b , . . . , b s − ) Implicit c − = 0 , c T = ( c , c , . . . , c s − ) . Methods for which c j = 0, j = 0 , . . . , s − a = − a j = 0, j = 1 , . . . , s −
1. We refer to Appendix A for a brief survey of someIMEX multistep methods, and to [1, 2, 4, 19, 21, 26, 45, 47] for further details and additional schemes.In what follows, we rely on the equivalent vector-matrix notation,(3.3) u n +1 = − a T · U − ∆ tc T · ∂ x V − ∆ tc − ∂ x v n +1 ,v n +1 = − a T · V − ∆ tε (cid:0) b T · ∂ x p ( U ) + c T · V + c − v n +1 − b T · f ( U ) (cid:1) , where U = ( u n , . . . , u n − s +1 ) T , V = ( v n , . . . , v n − s +1 ) T , ∂ x p ( U ) = ( ∂ x p ( u n ) , . . . , ∂ x p ( u n − s +1 )) T and f ( U ) = ( f ( u n ) , . . . , f ( u n − s +1 )) T are s -dimensional vectors.Similarly to the one step scheme (2.1), we proceed by rewriting the multistep methods in the fullyexplicit vector form. To that aim, we observe that the second equation in (3.3) can be explicitly solvedin terms of v due to the linearity in the relaxation part. This gives(3.4) v n +1 = − ε ε + ∆ tc − a T · V − ∆ tε + ∆ tc − (cid:0) b T · ∂ x p ( U ) + c T · V − b T · f ( U ) (cid:1) = − a T · V + (cid:18) a T − ε a T ε + ∆ tc − − ∆ tc T ε + ∆ tc − (cid:19) · V − ∆ tb T ε + ∆ tc − · ( ∂ x p ( U ) − f ( U )) , and thus(3.5) v n +1 = − a T · V − ∆ tε + ∆ tc − (cid:0) c T − c − a T (cid:1) · V − ∆ tε + ∆ tc − b T · ( ∂ x p ( U ) − f ( U )) . Substituting equation (3.5) into the first one of (3.3) leads to(3.6) u n +1 = − a T · U − ∆ t (cid:0) c T − c − a T (cid:1) · ∂ x V ++ ∆ t c − ε + ∆ tc − (cid:0) c T − c − a T (cid:1) · ∂ x V + ∆ t c − ε + ∆ tc − b T · ∂ x ( ∂ x p ( U ) − f ( U ))= − a T · U − ∆ t ε ε +∆ tc − (cid:0) c T − c − a T (cid:1) · ∂ x V + ∆ t c − ε + ∆ tc − b T · ∂ x ( ∂ x p ( U ) − f ( U )) . Hence, we obtain the following system(3.7) u n +1 + a T · U ∆ t = − ε (cid:0) c T − c − a T (cid:1) ε + ∆ tc − · ∂ x V + ∆ tc − ε + ∆ tc − b T · ∂ x ( ∂ x p ( U ) − f ( U )) ,v n +1 + a T · V ∆ t = − (cid:0) c T − c − a T (cid:1) ε + ∆ tc − · V − ε + ∆ tc − b T · ( ∂ x p ( U ) − f ( U )) , which is the generalization of system (2.4) to an s -step IMEX scheme.Our aim now is to show that, similarly to the simple first order method analyzed in Section 2, theabove discretization for a small value of ∆ t corresponds, up to first order in time, to a modified hyperbolic MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 7 problem where the characteristic speeds are bounded even in the limit ε →
0. More precisely, see also [3],for a smooth solution in time, by Taylor series expansion about t = t n , we have U = e u (cid:12)(cid:12) t = t n − ∆ tJ∂ t u (cid:12)(cid:12) t = t n + . . . + ( − s ∆ t s s ! J s ∂ st u (cid:12)(cid:12) t = t n + O (∆ t s +1 ) V = e v (cid:12)(cid:12) t = t n − ∆ tJ∂ t v (cid:12)(cid:12) t = t n + . . . + ( − s ∆ t s s ! J s ∂ st v (cid:12)(cid:12) t = t n + O (∆ t s +1 ) , where e is a vector of ones in R s , J = (0 , . . . , s − T , ∂ qt , q = 1 , . . . , s denotes the q -derivative andthe vector powers must be understood component-wise. Similarly, we can expand u n +1 and v n +1 about t = t n . Therefore, we obtain u n +1 + a T · U ∆ t = (1 + a T · e )∆ t u (cid:12)(cid:12) t = t n + (1 − a T · J ) ∂ t u (cid:12)(cid:12) t = t n + O (∆ t ) , (3.8) v n +1 + a T · V ∆ t = (1 + a T · e )∆ t v (cid:12)(cid:12) t = t n + (1 − a T · J ) ∂ t v (cid:12)(cid:12) t = t n + O (∆ t ) . From the order conditions we have1 + a T · e = 0 , − a T · J = b T · e = c − + c T · e =: β, and subsequentely (cid:0) c T − c − a T (cid:1) · e = c T · e + c − = β. Thus, scheme (3.7) for small values of ∆ t can be considered as a first order approximation of the followingmodified system(3.9) ∂ t u + ε ε + ∆ tc − ∂ x v + ∆ tc − ε + ∆ tc − ∂ x f ( u ) = ∆ tc − ε + ∆ tc − ∂ xx p ( u ) ∂ t v + 1 ε + ∆ tc − ∂ x p ( u ) = − ε + ∆ tc − ( v − f ( u )) , where the factor β simplifies in all the terms.Note that, the above system has exactly the same structure as system (2.7) (with c − = 1 and α = 1) which was derived from to the first order time discretization. As a consequence, under thesame simplification assumptions (2.9) on the fluxes f ( u ) and p ( u ), the eigenvalues of the hyperbolic partcorrespond to(3.10) Λ ± (∆ t, ε ) = 12 (cid:18) γ (1 − θ ) ± q γ (1 − θ ) + 4 ε − θ (cid:19) , where(3.11) θ (∆ t, ε ) := ε ε + ∆ tc − . Thus, the bounds for the characteristic velocities are the same as for the first order scheme and we getthe limit cases Λ ± (∆ t,
0) = 12 ( γ ± | γ | ) , Λ ± (0 , ε ) = ± ε . Asymptotic preserving property for α = 1 . Now, we study the capability of the schemes (3.3)to become a consistent discretization of the limit system (1.6). To that aim, letting ε →
0, in thereformulated scheme (3.7), we get from the first equation(3.12) u n +1 + a T · U ∆ t = b T · ∂ x ( ∂ x p ( U ) − f ( U )) , which corresponds to the explicit multistep scheme applied to the limiting convection-diffusion equation(1.6). For this reason, from now on, we refer to this class of IMEX-LM schemes as AP-explicit methods.On the other hand, we have for the second equation c − v n +1 + c T · V = − b T · ( ∂ x p ( U ) − f ( U )) , or equivalently(3.13) v n +1 = − c T c − · V − b T c − · ( ∂ x p ( U ) − f ( U )) . GIACOMO ALBI, GIACOMO DIMARCO, AND LORENZO PARESCHI
Let us observe that in order to have at time t n +1 a consistent projection over the asymptotic limit acondition over the states U = ( u n , .., u n − s +1 ) T and V = ( v n , .., v n − s +1 ) T should be imposed. This canbe resumed by saying that the vector of the initial data should be well prepared to the asymptotic state .This means that for the first variable(3.14) u n − j = ¯ u n − j + ˜ u n − jε , lim ε → ˜ u n − jε = 0 , j = 0 , ..., s − , where ¯ u n − j is a consistent solution of the limit system (1.6) while ˜ u n − jε is a perturbation that disappearsin the limit. An analogous relation should hold true for the second variable v , i.e.(3.15) v n − j = ¯ v n − j + ˜ v n − jε , lim ε → ˜ v n − jε = 0 , j = 0 , ..., s − , where ¯ v n − j = f ( u n − j ) − ∂ x p ( u n − j ), j = 0 , . . . , s − v n − jε is a perturbation that disappears in the limit. Under this assumption, as a consequence of theorder conditions, relation (3.13) is an s -order approximation of the asymptotic limit v = f ( u ) − ∂ x p ( u )and therefore, at subsequent time steps, the numerical solution is guaranteed to satisfy (3.14)-(3.15). Ifsuch conditions are not imposed on the initial values, then the numerical solution may present a spuriousinitial layer and deterioration of accuracy is observed. In particular, for IMEX-BDF methods, expression(3.13) simplifies to(3.16) v n +1 = − b T c − · ( ∂ x p ( U ) − f ( U )) . This shows that, even for non well prepared initial data in v but only in u we obtain an s -order approx-imation of the equilibrium state and subsequently of the numerical solution for all times. This strongerAP property of IMEX-BDF methods is satisfied also in the asymptotic limits analyzed for the variousschemes in the sequel of the manuscript. We refer to [18, 19, 43] for more detailed discussions.3.2. AP-explicit methods in the general case: α ∈ [0 , . The IMEX-LM scheme in vector formreads(3.17) u n +1 = − a T · U − ∆ tc T · ∂ x V − ∆ tc − ∂ x v n +1 ,v n +1 = − a T · V − ∆ tε α b T · ∂ x p ( U ) − ∆ tε α (cid:0) c T · V + c − v n +1 − b T · f ( U ) (cid:1) . Again, we rewrite the second equation by solving it in terms of v n +1 as follows v n +1 = − a T · V − ∆ tε α + ∆ tc − (cid:0) c T − c − a T (cid:1) · V ++ ∆ tε α + ∆ tc − b T · f ( U ) − ∆ tε − α ε α + ∆ tc − b T · ∂ x p ( U ) , and using this solution in the first equation, we obtain the explicit scheme(3.18) u n +1 + a T · U ∆ t = − ε α (cid:0) c T − c − a T (cid:1) ε α + ∆ tc − · ∂ x V − ∆ tc − b T ε α + ∆ tc − · ( ∂ x f ( U ) − ε − α ∂ xx p ( U )) ,v n +1 + a T · V ∆ t = − (cid:0) c T − c − a T (cid:1) ε α + ∆ tc − · V + b T ε α + ∆ tc − · (cid:0) f ( U ) − ε − α ∂ x p ( U ) (cid:1) . Considering, the same Taylor approximations (3.8) the above schemes corresponds up to first order intime to the modified system(3.19) ∂ t u + ε α ε α + ∆ tc − ∂ x v + ∆ tc − ε α + ∆ tc − ∂ x f ( u ) = ∆ t c − ε − α ε α + ∆ tc − ∂ xx p ( u ) ,∂ t v + ε − α ε α + ∆ tc − ∂ x p ( u ) = − ε α + ∆ tc − ( v − f ( u )) . Clearly, system (3.19) has again the same structure as (2.7). As a consequence, under the same simplifi-cation assumptions (2.9), the eigenvalues of the hyperbolic part correspond to(3.20) Λ α ± (∆ t, ε ) = 12 (cid:16) γ (1 − θ α ) ± p γ (1 − θ α ) + 4 ε − α θ α (cid:17) , MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 9 where θ α is defined as follows(3.21) θ α (∆ t, ε ) := ε α ε α + ∆ tc − , and the bounds for the characteristic velocities for ε = 0 and ∆ t = 0 areΛ α ± (∆ t,
0) = 12 ( γ ± | γ | ) , Λ α ± (0 , ε ) = ± ε α . Asymptotic preserving property for α ∈ [0 , . We consider now the analogous asymptotic preserv-ing property proved for the schemes (3.3) in the case of the schemes (3.17). Namely, we want to showthat (3.17) becomes a consistent discretization of the limit system (1.5) when ε →
0. Taking scheme(3.18), we get from the first equation(3.22) u n +1 = − a T · U − ∆ tb T · ∂ x f ( U ) , which is a standard explicit multistep discretization of the asymptotic hyperbolic limit, i.e. of equation(1.5). On the other hand, we have for the second equation c − v n +1 + c T · V = b T · f ( U ) , or equivalentely(3.23) v n +1 = − c T c − · V + b T c − · f ( U ) . As a consequence of the order conditions, equation (3.23) defines an s -order consistent approximationof the asymptotic limit v = f ( u ) provided that the vector of the initial data is well prepared. Theseconditions are the analogous of (3.14) and (3.15) except that now ¯ v n − j = f ( u n − j ), j = 0 , . . . , s − v n +1 = b T c − · f ( U ) , so that, even for non well prepared initial data in v but only in u we obtain an s -order approximation ofthe equilibrium state and subsequently of the numerical solution for all times.3.3. Removing the parabolic stiffness: AP-implicit methods.
Although, the schemes developedin the previous Section overcome the stiffness related to the scaling factor ε , there is another stiffness thatmay appear in the equations close to the asymptotic limit. In fact, as shown previously, all the schemesoriginate a fully explicit scheme in the limit.In diffusive regimes, this typically leads to the time step restriction ∆ t = O (∆ x ) when(3.25) ε − α ∆ tc − / ( ε α + ∆ tc − ) = O (1) , (see the diffusion coefficient in equation (3.19)). Therefore, for small ε and in the case of α ≃
1, the mainstability restriction is due to the second order term of the Chapmann–Enskog expansion (1.4). Note that,beside the case α = 1 and ε → α = 1 as soon as (3.25) holds true.For this reason, we modify the partitioning of the system taking also ∂ x p ( u ) implicit in the secondequation as follows u n +1 = − a T · U − ∆ tc T · ∂ x V − ∆ tc − ∂ x v n +1 , (3.26) v n +1 = − a T · V − ∆ tε α (cid:0) c T · V + c − v n +1 − b T · f ( U ) (cid:1) − ∆ tε α (cid:0) c T · ∂ x p ( U ) + c − ∂ x p ( u n +1 ) (cid:1) . We can still solve the second equation in v , to get v n +1 = − a T · V − ∆ tε α + ∆ tc − (cid:0) c T − c − a T (cid:1) · V ++ ∆ tε α + ∆ tc − b T · f ( U ) − ∆ tε − α ε α + ∆ tc − (cid:0) c T · ∂ x p ( U ) + c − ∂ x p ( u n +1 ) (cid:1) , which, inserted into the first equation of (3.26) yields the IMEX formulation(3.27) u n +1 + a T · U ∆ t = − ε α (cid:0) c T − c − a T (cid:1) ε α + ∆ tc − · ∂ x V − ∆ tc − ε α + ∆ tc − b T · ∂ x f ( U )+ ε − α ∆ tc − ε α + ∆ tc − (cid:0) c T · ∂ xx p ( U ) + c − ∂ xx p ( u n +1 ) (cid:1) ,v n +1 + a T · V ∆ t = − (cid:0) c T − c − a T (cid:1) ε α + ∆ tc − · V + 1 ε α + ∆ tc − b T · f ( U ) − ε − α ε α + ∆ tc − (cid:0) c T · ∂ x p ( U ) + c − ∂ x p ( u n +1 ) (cid:1) . Note that, except for the case in which p ( u ) is linear, in general, the first equation in (3.27) requires theadoption of a suitable solver for nonlinear problems to compute u n +1 . By the same arguments of in theprevious Sections, for small values of ∆ t , the scheme (3.27) corresponds up to first order to the modifiedsystem (3.19). Thus, under the same simplification assumptions (2.9), the eigenvalues of the hyperbolicpart are given by (3.20).3.3.1. Asymptotic preserving property for the AP-implicit methods.
Finally, we conclude our analysis bystudying the asymptotic preserving property. As we will see the main difference is that in the asymptoticlimit the diffusive terms are integrated implicitly. For this reason we refer to this class of IMEX-LMschemes as
AP-implicit methods. We first consider the case α = 1. Taking the reformulated IMEXscheme (3.27) and letting ε → α = 1, gives(3.28) u n +1 + a T · U ∆ t = − b T · ∂ x f ( U ) + c T · ∂ xx p ( U ) + c − ∂ xx p ( u n +1 ) , which correspond to the IMEX multistep scheme applied to the limiting convection diffusion problemwhere the diffusion term is treated implicitly [1, 2]. Note that, for the second equation we have c − v n +1 + c T · V = b T · f ( U ) − c T · ∂ x p ( U ) − c − ∂ x p ( u n +1 ) , or equivalently(3.29) v n +1 = − c T c − · V + b T c − · f ( U ) − c T c − · ∂ x p ( U ) − ∂ x p ( u n +1 ) , which, under the assumption of well prepared initial values, as a consequence of the order conditions,corresponds to a s -order approximation of the equilibrium projection v = f ( u ) − ∂ x p ( u ).On the other hand, in the case α ∈ [0 , u n +1 = − a T · U − ∆ tb T · ∂ x f ( U ) . We will not discuss further this limit system, but we emphasize that (3.30) is obtained as the limit of theimplicit-explicit scheme (3.27) whereas (3.22) is obtained as the limit of the explicit scheme (3.18).4.
Linear stability analysis
Monotonicity properties for IMEX-LM have been previously studied in [4, 19, 21, 26, 27]. Due to thewell-known difficulties in extending the usual stability analysis for linear systems to the implicit-explicitsetting most results are limited to the single scalar equation. In our case, however, the schemes arespecifically designed to deal with systems in the form (1.1), and we need therefore to tackle the stabilityproperties in such a case. Here we show that, in the case of IMEX-BDF methods, we can generalize someof the stability results for the single scalar equations to linear multiscale systems of the form ∂ t u + ∂ x v = 0 , (4.1) ∂ t v + 1 ε α ∂ x u = − ε α ( v − γu ) , where ε > γ > α ∈ [0 , ε , the above system when α = 1reduces to the convection-diffusion equation ∂ t u + γ∂ x u = ∂ xx u , whereas when α = 0, if γ <
1, yields thesimple advection equation ∂ t u + γ∂ x u = 0. To simplify notations, in the sequel we will assume γ = 1 and α >
0. The case α = 0 is rather classical and follows straightforwardly from our analysis. Under these MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 11 assumptions, the Chapman-Enskog expansion for small values of ε gives the limiting convection-diffusionequation(4.2) ∂ t u + ∂ x u = ε − α ∂ xx u + O ( ε α ) . In Fourier variables we get ˆ u ′ = − iξ ˆ v, (4.3) ˆ v ′ = − iξε α ˆ u − ε α (ˆ v − ˆ u ) , where, for example, ξ = sin(2 k ) / ∆ x if we use central differences and k is the frequency of the correspond-ing Fourier mode.The change of variables y = ˆ u , z = ε α ˆ v , λ I = iξ/ε α , λ R = 1 /ε α , λ = λ I + λ R ∈ C transforms thesystem into the problem y ′ = − λ I z, λ ∈ C , (4.4) z ′ = − ( λ I − λ R ε α ) y − λ R z. Let us note that the above problem is equivalent to the second order differential equation(4.5) y ′′ = − λ R y ′ + λ I ( λ I − λ R ε α ) y. AP-explicit methods.
We then apply an IMEX-LM method to system (4.4) as follows y n +1 = − s − X j =0 a j y n − j − λ I ∆ t s − X j = − c j z n − j (4.6) z n +1 = − s − X j =0 a j z n − j − ( λ I − λ R ε α )∆ t s − X j =0 b j y n − j − λ R ∆ t s − X j = − c j z n − j . (4.7)In the case of IMEX-BDF methods the first equation (4.6) permits to write z n +1 = − tc − λ I y n +1 + s − X j =0 a j y n − j and more in general for j = 0 , . . . , s − z n − j = − tc − λ I y n − j + s − X h =0 a h y n − j − h − ! . Thus, by direct substitution into the second equation (4.7) we obtain a discretization to (4.5) in the form y n +1 + s − X j =0 a j y n − j (1 + λ R ∆ tc − ) = − s − X j =0 a j y n − j + s − X h =0 a h y n − j − h − ! + λ I ( λ I − λ R ε α )∆ t c − s − X j =0 b j y n − j . Finally, we can rewrite the resulting scheme as y n +1 = − s − X j =0 a j y n − j −
11 + λ R ∆ tc − s − X j =0 a j y n − j + s − X h =0 a h y n − j − h − ! (4.9) + λ I ( λ I − λ R ε α )∆ t c − λ R ∆ tc − s − X j =0 b j y n − j . Note that, the above IMEX-LM in the limit ε → α = 1 leads to the reduced scheme(4.10) y n +1 = − s − X j =0 a j y n − j − ∆ t ( i + ξ ) ξ s − X j =0 b j y n − j , which corresponds to an explicit LMM for the convection-diffusion equation. Figure 2.
AP-explicit methods. Stability regions of IMEX-BDF methods in terms of z I = iξ ∆ t/ε α and z R = ∆ t/ε α . The different contour lines correspond to differentvalues of the scaling parameter ε α .The characteristic equation for scheme (4.9) reads(4.11) ̺ ( ζ ) + 11 + z R c − σ ( ζ ) − z I ( z I − z R ε α ) c − z R c − σ ( ζ ) = 0 , with z R = λ R ∆ t , z I = λ I ∆ t and ̺ ( ζ ) = ζ s + s − X j =0 a j ζ n − j , σ ( ζ ) = s − X j =0 a j ζ n − j + s − X h =0 a h ζ n − j − h − ! , σ ( ζ ) = s − X j =0 b j ζ n − j . Stability corresponds to the requirement that all roots of (4.11) have modulus less or equal one and thatall multiple roots have modulus less than one. In Figure 2 we plot the stability regions of AP-explicitIMEX-BDF schemes with respect to the variable z R and z I . The contour lines represent different valuesof the scaling parameter ε α . Note that, since we are assuming to use central differences, the stabilityregions are inversely proportional to the value of ε α , since ε − α measures the strength of the diffusive termin agreement with (4.2). As expected, as the order of the methods increase the corresponding stabilityrequirements become stronger and the various stability regions are reduced. MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 13
Figure 3.
AP-implicit methods. Stability regions of IMEX-BDF methods in terms of z I = iξ ∆ t/ε α and z R = ∆ t/ε α . The different contour lines correspond to differentvalues of the scaling parameter ε α .4.2. AP-implicit methods.
Next we apply an IMEX-LM method to (4.4) in the AP-implicit form y n +1 = − s − X j =0 a j y n − j − λ I ∆ t s − X j = − c j z n − j (4.12) z n +1 = − s − X j =0 a j z n − j + λ R ε α ∆ t s − X j =0 b j y n − j − ∆ t s − X j = − c j ( λ R z n − j + λ I y n − j ) . (4.13)Thus, restricting to IMEX-BDF methods, by direct substitution of (4.8) into the second equation (4.13)we get y n +1 + s − X j =0 a j y n − j (1 + λ R ∆ tc − ) = − s − X j =0 a j y n − j + s − X h =0 a h y n − j − h − ! − λ I λ R ∆ t c − ε α s − X j =0 b j y n − j + ( λ I ∆ tc − ) y n +1 , or equivalently y n +1 = − s − X j =0 a j y n − j −
11 + λ R ∆ tc − s − X j =0 a j y n − j + s − X h =0 a h y n − j − h − ! (4.14) − λ I λ R ∆ t c − ε α λ R ∆ tc − s − X j =0 b j y n − j + ( λ I ∆ tc − ) λ R ∆ tc − y n +1 . Now, the above LMM method in the limit ε → α = 1 leads to the scheme(4.15) y n +1 = − s − X j =0 a j y n − j − ∆ t iξ s − X j =0 b j y n − j − ∆ tξ c − y n +1 , which corresponds to an implicit-explicit IMEX-BDF scheme for the convection-diffusion equation.The characteristic equation associated to scheme (4.14) takes the form(4.16) ̺ ( ζ ) + 11 + z R c − σ ( ζ ) + z I z R ε α c − z R c − σ ( ζ ) − z I c − z R c − ζ s = 0 . We report in Figure 3 the stability regions of the AP implicit IMEX-BDF methods with respect to thevariable z R and z I . The contour lines, as for the AP explicit case, represent different values of the scalingparameter ε α . The second order method is uniformly stable when ε α < .
5, all other methods showbetter stability properties compared to the AP-explicit case, in particular for large values of z R . Againthe stability regions diminish for increasing values of ε α , in agreement with the limit problem (4.2), andas the order of the methods increases. 5. Space discretization
In this section we briefly discuss the space discretization adopted in the numerical examples. For thehyperbolic fluxes, we consider a WENO method of order five [48] combined with a Rusanov flux. Westress that the space discretization is not constructed over the original discretized systems, namely (3.1),(3.17) and (3.26). Instead, we introduce the space discretization on the reformulated systems (3.7) forthe AP-explicit with α = 1, on (3.18) for the AP-explicit with α ∈ [0 ,
1) and on (3.27) for the AP-implicit case. In fact, the adopted IMEX partitioning of the system, which guarantees boundedness ofthe eigenvalues, is of paramount importance to avoid instabilities of the fluxes and excessive numericaldissipation typical of diffusive scaling limits [31, 40, 41]. As a consequence, the numerical diffusion ischosen accordingly to (3.20) in the numerical fluxes reported below.Given a generic flux function F ( Q ) of Q ∈ R n , we first reconstruct the unknown values at the interfaces Q − , Q + and successively we employ the numerical Rusanov flux defined as follows(5.1) H ( Q − , Q + ) = 12 (cid:2) F ( Q + ) + F ( Q − ) − Θ( F ′ ( q )) S ( Q + − Q − ) (cid:3) , Θ( F ′ ( Q )) = max Q ∈ [ Q − ,Q + ] {| λ ( F ′ ( Q )) |} where max q ∈ [ Q − ,Q + ] {| λ ( F ′ ( Q )) |} represents the maximum modulus of the eigenvalues of the Jacobianmatrix F ′ ( Q ) and S ∈ R n × n a transformation matrix. Hence, for systems (3.7), (3.18), (3.27) andaccording to (3.10), (3.20), this value will depend on the scaling parameter ε and the choice of thediscretization steps. In particular for the general hyperbolic system (1.1), independently on the scalingfactor α , we have two unknowns Q = ( u, v ) T and three fluxesˆ f i + = 12 h ( f ( u + i + ) + f ( u − i + )) − Θ( u, v )( u + i + − u − i + ) i , (5.2) ˆ v i + = 12 h ( v + i + + v − i + ) − Θ( u, v )( u + i + − u − i + ) i (5.3)and ˆ p i + = 12 h ( p ( u + i + ) + p ( u − i + )) − Θ( u, v )( v + i + − v − i + ) i , (5.4) MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 15 where to comply with (5.1) we consider(5.5) S = , Θ( u, v ) = 12 (cid:16) γ (1 − θ α ) ± p γ (1 − θ α ) + 4 ε − α θ α (cid:17) , γ = f ′ ( u ) . The generic variables w reconstructed at the grid interfaces ( i + ) and ( i − ) respectively on the right w − i + and on the left side w + i − are given by(5.6) w − i + = X r =0 ω r w ( r ) i + , w + i − = X r =0 ω r ˜ w ( r ) i + with weights(5.7) ω r = α r P s =0 α s , α r = d r ( ε + β r ) , ˜ ω r = ˜ α r P s =0 ˜ α s , ˜ α r = ˜ d r ( ε + β r ) , and with standard smoothness indicators β = 1312 ( w i − w i +1 + w i +2 ) + 14 (3 w i − w i +1 + w i +2 ) β = 1312 ( w i − − w i + w i +1 ) + 14 ( w i − + w i +1 ) β = 1312 ( w i − − w i − + w i ) + 14 (3 w i − − w i − + w i ) , where ε = 10 − , d = 3 /
10 = ˜ d , d = 3 / d and d = 1 /
10 = ˜ d . Finally, the values w ( r ) i + and w ( r ) i + represent the third order reconstructions of the pointwise values ¯ w i . These are obtained throughthe formulas(5.8) w ( r ) i + = X j =0 c rj ¯ w i − r + j , w ( r ) − = X j =0 ˜ c rj ¯ w i − r + j , r = 0 , , w i − r + j are the pointwise values of the unknown evaluated at the points S r ( i ) = { x i − r , .., x i − r +2 } r = 0 , ,
2. Since, we use equispaced grid points, the coefficients c rj can be precomputed and their valuesare reported in Table 1. Table 1.
Coefficients c rj for the 5-th order WENO reconstruction on equispaced grid points. r j = 0 j = 1 j = 2-1 11/6 -7/6 1/30 1/3 5/6 -1/61 -1/6 5/6 1/32 1/3 -7/6 11/6 In addition, we have a second order term ∂ xx p ( u ) which, for the AP-explicit case, may be discretizedby two consecutive application of the Rusanov flux with WENO reconstruction of the state variablesand with numerical diffusion Θ( u, v ) fixed equal to zero, or by a specific space discretization which isconsistent with the limit problem. For example, by a sixth order finite difference formula(5.9) ∂ xx p ( u ( x i )) ≃ ap ( u i − ) + bp ( u i − ) + cp ( u i − ) + dp ( u i ) + cp ( u i +1 ) − bp ( u i +2 ) + ap ( u i +3 )(∆ x ) with a = 1 / , b = − / , c = 3 / , d = − /
18. This latter approach has been adopted in the case ofAP-implicit schemes, since the term ∂ xx p ( u ) is implicit and we want to avoid nonlinearities induced bythe WENO reconstructions. Numerical validation & applications
In this Section, we present different numerical tests to validate the analysis performed in the previousSections. In particular, we report results for the IMEX linear multistep from order two up to orderfive both for the AP-explicit and for the AP-implicit formulations. For the details about the IMEX-LMmethods used we refer to Appendix A. We remark that other IMEX-LM methods can be included as wellin the present formulation, see for example [45, 47]. In all test cases, the initial data has been chosenwell prepared and the IMEX-LM methods have been initialized with a third order IMEX Runge-Kuttascheme (see [7]) with a time step which satisfies the accuracy constraints.6.1.
Test 1. Numerical convergence study for a linear problem.
We consider the following linearhyperbolic model with diffusive scaling for α = 1(6.1) ∂ t u + ∂ x v = 0 ,∂ t v + 1 ε ∂ x u = − ε ( v − γu ) , where γ >
0. In the diffusive limit ε → v = γu − ∂ x u, substituting into the first equation this gives the limiting advection-diffusion equation(6.2) ∂ t u + γ∂ x u = ∂ xx u. In particular, we consider the model (6.1) solved on the domain x ∈ [0 , γ = 1 and periodicboundary conditions and with smooth initial data given by(6.3) u ( x,
0) = sin(2 xπ ) , v ( x,
0) = sin(2 xπ ) − cos(2 xπ ) . Note that, the initial data is well prepared, in the sense that v ( x,
0) = u ( x, − ∂ x u ( x, ε = 1 , . , . , .
001 by measuring the space and time L -error of the numerical solutions computed byusing as reference solution the thinner grids. Namely, given a coarser grid ∆ x = 1 /N with N = 128 weconsider(6.4) ∆ x k +1 = ∆ x k / k = 1 , .., . The time step for AP-explicit methods is chosen as ∆ t = λ ∆ x max { ε, ∆ x } , with λ = 0 .
25. Namely, thelargest between the CFL condition imposed by the hyperbolic part and by the limiting parabolic part ofthe equations. Instead, for implicit schemes we choose ∆ t = λ ∆ x max { ε, } since the diffusion term isintegrated implicitly in the limit.The local truncation error is measured for the two components u and v as follows E k ∆ x, ∆ t ( u ) = | u k ∆ x, ∆ t ( · , T ) − u k − x, ∆ t ( · , T ) | , E k ∆ x, ∆ t ( v ) = | v k ∆ x, ∆ t ( · , T ) − v k − x, ∆ t ( · , T ) | , and the order of convergence is estimated by computing the rate between two L errors of distinctnumerical solutions(6.5) p k ( u ) = log ( k E k − x, ∆ t ( u ) k / k E k ∆ x, ∆ t ( u ) k ) , p k ( v ) = log ( k E k − x, ∆ t ( v ) k / k E k ∆ x, ∆ t ( v ) k ) . The analysis is performed for several different IMEX linear multi-step schemes from second to fifth order(see Appendix A). In particular, we focus on the BDF and the TVB classes of IMEX multistep methodsthanks to their favorable stability properties (see [26,27] for details and derivation). We report in Table 2the space-time L errors and the relative rates of convergence for increasing size of the meshes considering N = 2 k points with k = 8 , ..,
11 for the u variable, while the space-time L errors and the relative ratesof convergence for the v variable are shown in Table 3 for the AP-explicit schemes. In Table 4 and 5, wereport the corresponding L errors and rates of convergence for the AP-implicit schemes. In all cases wecan conclude that the expected orders of convergence are achieved by the schemes for the different valuesof the asymptotic parameter ε and that the behavior of the schemes outperforms the corresponding IMEXRunge-Kutta methods for the non conserved quantity v (see Table 2 in [8] for example). In particular,we observe the tendency of 4-th order methods to achieve higher then expected convergence rates on theconserved quantity u for moderate values of the stiffness parameter. The same tendency, on both variables u and v , is observed for 5-th order methods close to the diffusion limit particularly in the AP-implicit case. MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 17
On the contrary, the 4-th order schemes in the AP-explicit implementation produce a slight deteriorationon the non conserved quantity v which is not observed in the AP-implicit setting. The SG(3,2) schemealso suffers of a slight deterioration of accuracy close to the diffusion limit in the AP-explicit form. Itshould be noted that scheme in AP-explicit form close to the fluid limit have a smaller truncation errorwith respect to time thanks to the CFL condition ∆ t = O (∆ x ). This can be observed by comparing the L -errors of the AP-explicit and AP-explicit formulations for ε = 0 . ε for the AP-explicit and AP-implicit methods. These plots emphasize that the convergence rate for allschemes is almost uniform. Table 2. L error and estimated convergence rates for u in the AP-explicit case. ε = 1 ε = 0 . ε = 0 . ε = 0 . N k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate I M E X - S G ( , )
128 9.9402e-05 – 8.2324e-05 – 0.00019196 – 1.0662e-07 –256 3.0565e-05 1.7014 2.5215e-05 1.707 5.7846e-05 1.7305 4.262e-08 1.3229512 8.3499e-06 1.8721 6.8755e-06 1.8748 1.5671e-05 1.8841 1.1871e-08 1.8441 I M E X - B D F
128 4.4737e-05 – 3.6789e-05 – 8.9909e-05 – 9.2432e-08 –256 1.4423e-05 1.6331 1.1759e-05 1.6455 2.8204e-05 1.6726 2.4676e-08 1.9053512 4.0033e-06 1.8491 3.2494e-06 1.8555 7.7212e-06 1.869 6.0694e-09 2.0235 I M E X - T V B ( , )
128 7.3451e-08 – 1.2994e-07 – 5.0813e-06 – 3.1308e-08 –256 1.2953e-08 2.5035 1.8201e-08 2.8358 7.5974e-07 2.7416 3.8147e-09 3.0369512 1.8478e-09 2.8094 2.4055e-09 2.9197 1.0251e-07 2.8898 4.4859e-10 3.0881 I M E X - B D F
128 8.7902e-08 – 1.6714e-07 – 6.6227e-06 – 5.8099e-08 –256 1.5388e-08 2.514 2.4439e-08 2.7738 9.8974e-07 2.7423 6.0486e-09 3.2638512 2.1902e-09 2.8127 3.2815e-09 2.8968 1.3293e-07 2.8964 6.5821e-10 3.2 I M E X - T V B ( , )
128 2.2234e-08 – 6.9585e-09 – 1.2029e-06 – 2.0263e-07 –256 8.4006e-10 4.7262 1.7448e-10 5.3177 9.1633e-08 3.7145 5.8212e-09 5.1214512 4.7195e-11 4.1538 3.7747e-12 5.5305 6.0883e-09 3.9118 2.2316e-10 4.7052 I M E X - B D F
128 2.3209e-08 – 7.4749e-09 – 4.5554e-07 – 2.449e-08 –256 6.8151e-10 5.0898 2.1477e-10 5.1212 3.3431e-08 3.7683 1.1948e-09 4.3573512 2.9394e-11 4.5352 5.499e-12 5.2875 2.2255e-09 3.909 6.3684e-11 4.2297 I M E X - T V B ( , )
128 1.043e-08 – 4.1161e-08 – 9.4337e-09 – 1.9208e-07 –256 3.2924e-10 4.9854 1.2742e-09 5.0137 3.3621e-10 4.8104 2.4637e-09 6.2848512 1.4718e-11 4.4835 4.7629e-11 4.7416 1.0529e-11 4.9969 5.4809e-11 5.4903 I M E X - B D F
128 4.9951e-08 – 3.49e-08 – 4.906e-09 – 1.7525e-08 –256 1.5618e-09 4.9992 1.0894e-09 5.0016 1.5768e-10 4.9594 4.2623e-10 5.3616512 5.6577e-11 4.7869 3.2187e-11 5.0809 4.6748e-12 5.076 1.1341e-11 5.232
Test 2: Riemann problem for the linear model.
Next, we consider a Riemann problem definedon the space interval [0 ,
4] with discontinuous initial data as follows(6.6) u L = 4 . , v L = 0 , ≤ x ≤ ,u R = 2 . , v R = 0 , < x ≤ . Table 3. L error and estimated convergence rates for v in the AP-explicit case. ε = 1 ε = 0 . ε = 0 . ε = 0 . N k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate I M E X - S G ( , )
128 4.1297e-05 – 0.0015868 – 0.30746 – 0.017621 –256 1.3046e-05 1.6624 0.00049005 1.6951 0.092824 1.7278 0.006956 1.341512 3.6092e-06 1.8538 0.00013416 1.869 0.025173 1.8826 0.0019759 1.8158 I M E X - B D F
128 5.3468e-05 – 0.0010541 – 0.14541 – 0.015341 –256 1.7643e-05 1.5996 0.00034308 1.6194 0.045795 1.6669 0.0040821 1.91512 4.9529e-06 1.8328 9.5659e-05 1.8426 0.012564 1.8658 0.001048 1.9616 I M E X - T V B ( , )
128 7.2288e-07 – 8.2789e-06 – 0.0083555 – 0.0044838 –256 1.1051e-07 2.7096 1.2742e-06 2.6999 0.0012579 2.7317 0.00064242 2.8031512 1.5103e-08 2.8712 1.7412e-07 2.8714 0.0001703 2.8849 8.423e-05 2.9311 I M E X - B D F
128 6.3076e-07 – 8.0827e-06 – 0.010784 – 0.0087392 –256 9.8098e-08 2.6848 1.2676e-06 2.6727 0.0016207 2.7342 0.00099636 3.1328512 1.3503e-08 2.861 1.7436e-07 2.862 0.0002183 2.8923 0.00011713 3.0885 I M E X - T V B ( , )
128 8.1836e-08 – 8.5455e-08 – 0.0019902 – 0.033128 –256 3.443e-09 4.571 7.8768e-09 3.4395 0.00015204 3.7104 0.00099626 5.0554512 1.5925e-10 4.4343 6.285e-10 3.6476 1.0128e-05 3.9079 4.1573e-05 4.5828 I M E X - B D F
128 6.8621e-08 – 2.6004e-08 – 0.00076076 – 0.0045931 –256 2.4028e-09 4.8359 2.607e-09 3.3183 5.5621e-05 3.7737 0.00022733 4.3366512 5.3325e-11 5.4938 2.1977e-10 3.5683 3.7004e-06 3.9099 1.2649e-05 4.1676 I M E X - T V B ( , )
128 1.7977e-08 – 2.1366e-07 – 9.0073e-06 – 0.030281 –256 4.6141e-10 5.2839 6.7708e-09 4.9798 3.4995e-07 4.6858 0.00039906 6.2457512 6.7507e-11 2.773 1.182e-10 5.84 1.1336e-08 4.9481 9.4581e-06 5.3989 I M E X - B D F
128 1.2276e-07 – 1.533e-07 – 1.6471e-06 – 0.0022391 –256 3.8413e-09 4.9981 4.7927e-09 4.9993 5.7491e-08 4.8405 5.2309e-05 5.4197512 2.5154e-10 3.9327 1.9215e-10 4.6405 1.4584e-09 5.3009 1.3736e-06 5.251
For the above initial data, the linear hyperbolic system in the form (1.1) with zero-flux boundary condi-tions is solved and comparisons with different values of the relaxation parameter ε are shown. The sameproblem has been studied in [8] using IMEX Runge-Kutta schemes.In the limit ε →
0, the exact solution of the corresponding advection-diffusion equation is known andit reads(6.7) u ( x, t ) = 12 ( u L + u R ) + 12 ( u L − u R )erf (cid:18) t − x + 22 √ t (cid:19) , with erf( x ) the error function. We report in Figure 6 the numerical solution for u at final time T =0 .
25 computed using two different time integration schemes, namely the IMEX-BDF2, and the IMEX-TVB(4,4). We choose N = 80 points in space, and compare the AP-explicit approach (3.18), withthe AP-implicit one (3.27). The reference solution is computed with the IMEX-BDF5 scheme using∆ t ref = ∆ t/
10 and N = 200 space points. In order to preserve the CFL conditions the time steps for thedifferent regimes of ε are selected as follows: in the AP-explicit case ∆ t = λ ∆ x max { ε, ∆ x } and in theAP-implicit case ∆ t = λ ∆ x max { ε, } with λ = 0 .
4. In Figure 6, for the diffusive limit, we observe thatthe different schemes agree well with the exact solution (6.7). In the hyperbolic regime, ε = 0 .
5, the shock
MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 19
Table 4. L error and estimated convergence rates for u in the AP-implicit case. ε = 1 ε = 0 . ε = 0 . ε = 0 . N k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate I M E X - S G ( , )
128 0.0001009 – 8.4701e-05 – 0.00019279 – 0.015388 –256 3.1168e-05 1.6948 2.6145e-05 1.6958 5.8207e-05 1.7278 0.0047058 1.7093512 8.5342e-06 1.8687 7.1564e-06 1.8692 1.5785e-05 1.8826 0.0012828 1.8751 I M E X - B D F
128 4.5983e-05 – 3.8769e-05 – 9.0615e-05 – 0.0074794 –256 1.4986e-05 1.6175 1.2626e-05 1.6185 2.8544e-05 1.6666 0.0023976 1.6413512 4.1817e-06 1.8414 3.5217e-06 1.8421 7.8321e-06 1.8657 0.00066348 1.8535 I M E X - T V B ( , )
128 5.6622e-08 – 1.5426e-07 – 5.1657e-06 – 0.00024864 –256 9.4982e-09 2.5756 2.2969e-08 2.7476 7.7718e-07 2.7326 3.7993e-05 2.7103512 1.3103e-09 2.8577 3.121e-09 2.8796 1.0522e-07 2.8849 5.1683e-06 2.878 I M E X - B D F
128 7.3382e-08 – 1.8956e-07 – 6.7012e-06 – 0.0002871 –256 1.2326e-08 2.5737 2.8971e-08 2.7099 1.0064e-06 2.7352 4.443e-05 2.6919512 1.7158e-09 2.8448 3.9775e-09 2.8647 1.3555e-07 2.8923 6.0575e-06 2.8747 I M E X - T V B ( , )
128 2.9453e-08 – 2.5367e-08 – 1.2204e-06 – 2.9521e-05 –256 1.2448e-09 4.5645 7.248e-10 5.1292 9.3626e-08 3.7044 2.2475e-06 3.7153512 5.0075e-11 4.6356 1.6962e-11 5.4172 6.2417e-09 3.9069 1.5677e-07 3.8416 I M E X - B D F
128 2.582e-08 – 7.3488e-09 – 4.6271e-07 – 1.2261e-05 –256 8.5367e-10 4.9187 2.0595e-10 5.1572 3.4155e-08 3.7599 8.3808e-07 3.8709512 1.6737e-11 5.6726 1.69e-12 6.9291 2.2839e-09 3.9025 5.5639e-08 3.9129 I M E X - T V B ( , )
128 3.6171e-08 – 3.5169e-08 – 9.4689e-09 – 1.1944e-06 –256 1.1333e-09 4.9962 1.1385e-09 4.9491 3.3274e-10 4.8308 7.5898e-09 7.298512 5.1332e-11 4.4646 6.4841e-11 4.1341 1.2796e-11 4.7006 3.7672e-10 4.3325 I M E X - B D F
128 2.0507e-08 – 2.7439e-08 – 4.9118e-09 – 1.3221e-06 –256 6.2502e-10 5.0361 8.5434e-10 5.0053 1.5808e-10 4.9575 2.1912e-08 5.915512 2.5331e-11 4.625 3.7038e-11 4.5277 1.5295e-12 6.6914 1.0639e-09 4.3643 is correctly captured compared to the reference solution by both schemes with a slightly better resolutionin the case of IMEX-TVB(4,4) for the AP-implicit approach. Note that, no spurious oscillations areobserved for IMEX-BDF2 scheme, even if it does not satisfy any specific TVB stability property in thehyperbolic regime.6.3.
Test 3: Barenblatt solution for the porous media equation.
We then consider a nonlineardiffusion limit, by considering the following hyperbolic system with diffusive relaxation(6.8) ∂ t u + ∂ x v = 0 ,∂ t v + 1 ε ∂ x u = − ε α k ( u ) v. This problem has been previously studied in [33, 41]. Note that, when k ( u ) = 1, the system (6.8) is amodel of relaxing heat flow and, as ε →
0, it relaxes towards the heat equation(6.9) ∂ t u = ε α ∂ xx u, v = − ∂ x u. Table 5. L error and estimated convergence rates for v in the AP-implicit case. ε = 1 ε = 0 . ε = 0 . ε = 0 . N k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate k E k ∆ x, ∆ t k Rate I M E X - S G ( , )
128 4.4648e-05 – 0.00089276 – 0.30567 – 0.11134 –256 1.3846e-05 1.6891 0.00027568 1.6953 0.092283 1.7278 0.032643 1.7702512 3.7984e-06 1.866 7.5475e-05 1.8689 0.025026 1.8826 0.0087961 1.8918 I M E X - B D F
128 2.3413e-05 – 0.00042114 – 0.14368 – 0.060252 –256 7.6904e-06 1.6062 0.00013723 1.6177 0.045257 1.6667 0.0174 1.7919512 2.1539e-06 1.8361 3.8288e-05 1.8416 0.012418 1.8657 0.0046647 1.8992 I M E X - T V B ( , )
128 2.0289e-07 – 1.2756e-06 – 0.0081773 – 0.001989 –256 3.279e-08 2.6294 1.9536e-07 2.707 0.0012317 2.731 0.00026696 2.8974512 4.5301e-09 2.8557 2.6669e-08 2.8729 0.00016679 2.8845 3.4663e-05 2.9451 I M E X - B D F
128 2.3141e-07 – 1.3649e-06 – 0.010612 – 0.0023402 –256 3.8064e-08 2.6039 2.1348e-07 2.6766 0.0015951 2.7339 0.00032067 2.8674512 5.2868e-09 2.848 2.9369e-08 2.8617 0.00021487 2.8921 4.2074e-05 2.9301 I M E X - T V B ( , )
128 7.1151e-08 – 9.5652e-08 – 0.0019469 – 0.00028374 –256 2.6963e-09 4.7218 2.8851e-09 5.0511 0.0001488 3.7098 1.6152e-05 4.1348512 1.0634e-10 4.6642 3.3063e-10 3.1253 9.9148e-06 3.9076 1.0431e-06 3.9527 I M E X - B D F
128 6.5438e-08 – 1.7232e-08 – 0.00074557 – 9.082e-05 –256 2.1846e-09 4.9047 4.6192e-10 5.2213 5.4526e-05 3.7733 5.8265e-06 3.9623512 3.0347e-11 6.1697 3.0551e-11 3.9183 3.6281e-06 3.9097 3.7459e-07 3.9592 I M E X - T V B ( , )
128 8.0907e-08 – 1.4931e-07 – 8.7961e-06 – 3.1526e-05 –256 2.3712e-09 5.0926 4.7065e-09 4.9875 3.4175e-07 4.6859 2.7811e-07 6.8247512 1.4198e-10 4.0619 5.5157e-11 6.415 1.1008e-08 4.9563 9.1116e-09 4.9318 I M E X - B D F
128 4.6717e-08 – 1.0372e-07 – 1.6077e-06 – 9.3079e-06 –256 1.4517e-09 5.0081 3.1885e-09 5.0237 5.6036e-08 4.8425 1.6743e-07 5.7968512 8.9081e-12 7.3484 1.0165e-10 4.9711 1.4002e-09 5.3226 7.4625e-09 4.4878
On the other hand, by choosing k ( u ) = (2 u ) − , the limiting equation for this model results in the porousmedia equation(6.10) ∂ t u = ε α ∂ xx u , v = − u∂ x u. For this problem we consider the IMEX-BDF2 and the IMEX-BDF5 schemes for α = 0 by computingthe numerical solution with N = 80 mesh points, using the AP-explicit approach and selecting the timestep as ∆ t = λ ∆ x max { ε, ∆ x } with λ = 0 . ε → u ( x, t ) = r ( t ) " − (cid:18) xr ( t ) (cid:19) , | x | ≤ r ( t ) , , | x | > r ( t ) , where r ( t ) = [12( t + 1)] / , t ≥ x ∈ [ − , v ( x, t = 0) = 0. In Figure 7, this comparison is shown MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 21 -3 -2 -1 AP Explicit - u
IMEX-SG(3,2)IMEX-BDF2IMEX-TVB(3,3)IMEX-BDF3IMEX-TVB(4,4)IMEX-BDF4IMEX-TVB(5,5)IMEX-BDF5 10 -3 -2 -1 AP Explicit - v
IMEX-SG(3,2)IMEX-BDF2IMEX-TVB(3,3)IMEX-BDF3IMEX-TVB(4,4)IMEX-BDF4IMEX-TVB(5,5)IMEX-BDF5
Figure 4.
Order of convergence of AP-explicit methods in the L − norm using (6.5) forthe u variable (left) and the v variable (right). AP Implicit - u AP Implicit - v
Figure 5.
Order of convergence of AP-implicit methods in the L − norm using (6.5) forthe u variable (left) and the v variable (right).using ε = 10 − at time T = 3. The sharp front of the Barenblatt solution is very well captured by bothschemes without observing any relevant different by the increase order of accuracy of the IMEX-BDF5method. This is somewhat expected, since the CFL condition ∆ t = O (∆ x ), in practice, gives to theIMEX-BDF2 scheme the same accuracy of a fourth order method. Therefore, the accuracy barrier hereis represented by the 5-th order space discretization method.6.4. Test 4: Applications to the Ruijgrook–Wu model.
We consider, in this last test case, anapplication of the schemes to the so-called Ruijgrook–Wu model of rarefied gas dynamic [8,22,23,33,41,46](6.12)
M ∂ t f + + ∂ x f + = − Kn ( af + − bf − − cf + f − ) ,M ∂ t f − − ∂ x f − = 1 Kn ( af + − bf − − cf + f − ) , where f + and f − denote the particle density distribution at time t , position x and with velocity +1 and − Kn is the Knudsen number, M is the Mach number of the system and a , b and c arepositive constants which characterize the microscopic interactions . The local (Maxwellian) equilibrium Figure 6.
Test 2. Riemann problem for linear system with diffusive scaling (6.1) atfinal time T = 0 .
25. Hyperbolic regime with ε = 0 . ε = 10 − (right). Top row AP-explicit approach, bottom row AP-implicit approach. -10 -8 -6 -4 -2 0 2 4 6 8 1000.050.10.150.20.250.30.350.40.450.5 3.5 4 4.500.020.04 Figure 7.
Test 3. Comparison between the Barenblatt analytical solution of system(6.10), and the numerical solution obtained from system (6.8) with k ( u ) = (2 u ) − and ε = 10 − at time T = 3. MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 23 f ±∞ is characterized by(6.13) f + ∞ = bf −∞ a − cf −∞ . The macroscopic variables for the model are the density u and momentum v defined by(6.14) u = f + + f − , v = ( f + − f − ) /M. The nondimensional multiscale problem is obtained taking M = ε α and Kn = ε , the Reynolds number ofthe system is then defined as usual according to Re = M/Kn = 1 /ε − α . In macroscopic variables taking a = b = 1 / c = M = ε α this can be written as [22](6.15) ∂ t u + ∂ x v = 0 ,∂ t v + 1 ε α ∂ x u = − ε α (cid:20) v − (cid:0) u − ε α v (cid:1)(cid:21) . The model has the nice feature to provide nontrivial limit behaviors for several values of α including thecorresponding compressible Euler ( α = 0) limit and the incompressible Euler ( α ∈ (0 , α = 1) limits (see [22, 23]). In [41] the above model has been used under a similar but differentscaling. If we denote with the pair ˜ α, ˜ ε the scaling parameters in [41], these correspond to take M = ˜ ε , Kn = ˜ ε α . The two scaling can be made equivalent for α = 0, if we map the pair α, ε into ˜ α, ˜ ε as follows(6.16) ˜ α = 1 − αα , ˜ ε = ε α , α = 0 . Note, that the non-linearity on the source term depends both on u and v , so that we have f = f ( u, v ) in(1.1). Nevertheless, following the same strategy described in the previous Sections, our methods can beapplied in a straightforward way also in this situation.For the Ruijgrook–Wu model (6.15) it can be shown, via Chapman-Enskog expansion, that for α ∈ (1 / ,
1] so that 2 α > − α and small values of ε we have(6.17) v = 12 u − ε − α ∂ x u + O ( ε α ) . Then, the solution behavior is characterized by the viscous Burgers equation(6.18) ∂ t u + ∂ x (cid:18) u (cid:19) = ε − α ∂ xx u + O ( ε α ) . In the sequel, for
Test 4a we apply the AP-explicit approach with ∆ t = λ ∆ x max { ε, min { , ∆ x/ε − α }} , λ = 0 . Test 4b and
Test 4c we use the AP-implicitapproach with ∆ t = λ ∆ x max { ε, } and λ = 0 .
1. Therefore, for ε ≤ ε and α . The reference solution is computed usingthe IMEX-BDF5 scheme with time step ∆ t ref = ∆ t/
10 and N ref = 400 space points. The initial datafor the non conserved quantity v has been taken well prepared accordingly to (6.17). In the sequel wewill restrict to ˜ α ∈ [0 ,
1] as in [41], therefore, from (6.16), we have α ∈ [0 . , α >
0. We refer to [22] for results in the fluid limit α = 0. Test 4a. Riemann problem.
We select the Riemann problem characterized by (6.6) as initial data. Thenumerical parameters are also the same of the previous test case, therefore we fix the final time T = 0 . u ( x ) using two different time AP-explicit integration schemes, BDF2 andTVB(4,4) with N = 100 points in space.In Figure 8 we report in the left column the behavior in the rarefied (non stiff) regime ε = 0 . ε = 10 − , whereas top row depicts the behavior for α = 1 andbottom row α = 2 /
3. Even if the rarefied solutions for the different values of α are similar, we remarkthe difference of the solution profiles in the limit ε →
0. For α = 1 we obtain the classical viscous Burgerequation where the discontinuity is smeared out by the diffusion term, while for α = 2 / Figure 8.
Test 4a . Riemann problem for the non linear Ruijgrook–Wu model (6.15) atfinal time T = 0 .
25 with AP-explicit methods. Rarefied regime with ε = 0 . ε = 10 − (right). Top row α = 1, bottom row α = 2 / Test 4b. Propagation of a square wave.
Next, we consider the Ruijgrook–Wu model, in the space interval[ − . , .
5] with initial data defined as follows u ( x ) = ( | x | ≤ /
80 otherwise , v ( x ) = 0 , (6.19)where we account for reflecting boundary conditions, i.e. v = 0 , ∂ x u = 0 on the boundaries x = ± . α and ε and we solve the modelwith N = 100 space grid points using the AP-implicit IMEX-BDF4 method. We report in Figure 9 theevolution of the solution for the density u and the momentum v , respectively on the top and bottomrows. The first column represents the rarefied regime, for ε = 0 . , α = 1. In this regime, the transportpart dominates and the initial data propagates in the directions of the particles. This behavior is welldescribed by the method without spurious numerical oscillations. In the second column we have thehyperbolic limit for ε = 10 − , α = 2 /
3, corresponding to the inviscid Burger equation with a shockpropagating in the right direction. Even in this case the shock profile is well captured. Finally, in the lastcolumn we report the parabolic limit for ε = 10 − , α = 4 /
5, corresponding to a viscous Burger equation.As expected, the shock profile is regularized by the presence of the diffusive term.
MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 25 -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.80.91 -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.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.100.10.20.30.40.5 -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 9.
Test 4b . Non linear Ruijgrook–Wu model (6.15) with discontinuous initialdata (6.19), top column reports the density, bottom column the momentum. Left row ε = 0 . , α = 1 and final time T = 0 .
2, middle row ε = 10 − , α = 2 / T = 0 . ε = 10 − , α = 4 /
5, with final time T = 0 .
5. The AP-implicit formulation hasbeen used.
Test 4c. Anisotropy of the multiscale parameter α . In the last test case, we solve the model (6.15)considering a multiscale parameter α to be a function of the space x , whereas the relaxation parameter isfixed to ε = 10 − . This test aims at reproducing a realistic situation where the scaling terms may dependon the physical quantities and vary in the different regions of the computational domain. We report theresults obtained using the AP-implicit scheme with IMEX-BDF3 and the same discretization parametersof the previous test case. The initial date is defined in the space interval [ − . , .
5] as follows u ( x ) = ( | x | ≤ / . , v ( x ) = 0 . (6.20)In the left column of Figure 10, we report the value of function α ( x ) as a function of the space,varying between 0 . u ( x, t ) and v ( x, t ) showing the initial and final profiles, a similar test casehas been presented in [8] for IMEX Runge-Kutta methods.In the first row, we account a single variation of the regime from the hyperbolic to the parabolic α ( x ) = 1 − H ( x ) , where H ( x ) = 11 + exp( x/δ ) , δ = 0 . . At final time T = 0 .
05 we observe a rarefaction wave moving to the left (in the hyperbolic regime)and a smooth profile on the right (in the parabolic regime). Note that, the method captures well thecomplicated shock structure even at the interface between the two regions.Second row considers two variations of the the regime from hyperbolic to parabolic, choosing α ( x ) as α ( x ) = 12 −
12 ( H ( x + 0 . − H ( x − . , at final time is T = 0 . u ( x ) is blunted within the parabolicregime, whereas a shock and rarefaction waves emerge in the hyperbolic one, both well described by thenumerical method. -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.50.50.550.60.650.70.750.80.850.90.951 Parabolic regimeHyperbolic regime -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.50.50.550.60.650.70.750.80.850.90.951 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.500.511.5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.50.50.550.60.650.70.750.80.850.90.951
Parabolic regimeHyperbolic regime -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.50.50.550.60.650.70.750.80.850.90.951 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5-2.5-2-1.5-1-0.500.51
Figure 10.
Test 4c . Non linear Ruijgrook–Wu model with space dependent α and ε = 10 − . Left column shows the space variation of the multiscale parameter, whereasthe central and right columns depict the evolution of the density u and momentum v from the initial data u ( x ) , v ( x ) to final time. Top row accounts a single variation fromhyperbolic to parabolic and final time T = 0 .
05, whereas second row has two transitionsand final time T = 0 .
1. 7.
Conclusions
In this work we have developed a unified IMEX multistep approach for hyperbolic balance laws underdifferent scalings. These problems, inspired by the classical hydrodynamical limits of kinetic theory[13], are challenging for numerical methods since the nature of the asymptotic behavior is not knowna-priori and depends on the scaling parameters. Therefore, the schemes should be able to capturecorrectly asymptotic limits characterized by hyperbolic conservation laws as well as diffusive parabolicequations. A major difficulty in the schemes construction is represented by the unbounded growth of thecharacteristic speeds of the system in diffusive regimes. For these problems, we developed two differentkind of approaches, originating a problem reformulation with bounded characteristic speeds, and whichlead respectively to explicit (AP-explicit) or explicit-implict (AP-implicit) time discretizations of theasymptotic limit. Several numerical results for linear and non linear hyperbolic relaxation systems haveconfirmed that the IMEX multistep methods are capable to describe correctly the solution for a widerange of relaxation parameters and for different values of the scaling coefficient α . Compared to theIMEX Runge-Kutta approach developed in [8] the IMEX multistep schemes here constructed presentseveral advantages. In particular, it is possible to achieve easily high order accuracy and in general amore uniform behavior of the error with respect to the scaling parameters is observed. In addition,when dealing with computationally challenging problems such as the case of kinetic equations with stiffcollision terms its is possible to strongly reduce the number of evaluations of the most expensive part ofthe computation represented by the source term. Future research will be in the direction of extendingthe present results to the more difficult case of diffusion limits for non linear kinetic equations and, morein general, to the case of low Mach number limits and all Mach number flows [7, 25, 28, 29, 37]. MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 27
Appendix A. Order conditions for IMEX-LM methods and examples
In this appendix we give the details of the particular IMEX-LM methods used in the manuscript. Letus recall that an order p , s -step, IMEX-LM scheme is obtained provided that(A.1) 1 + s − X j =0 a j = 0 , − s − X j =0 ja j = s − X j =0 b j = s − X j = − c j ,
12 + s − X j =0 j a j = − s − X j =0 jb j = c − − s − X j =0 jc j . ...1 p ! + s − X j =0 ( − j ) p p ! a j = s − X j =0 ( − j ) p − ( p − b j = c − ( p − s − X j =0 ( − j ) p − ( p − c j , Moreover the following theorem holds true [4]
Theorem 1.
For an s -step IMEX scheme we have (1) If p ≤ s the p + 1 constraints of (A.1) are linearly independent, therefore there exist s -stepIMEX-LM schemes of order s . (2) An s -step IMEX-LM scheme has accuracy at most s . (3) The family of s -step IMEX-LM schemes of order s has s parameters. Listed below the IMEX-LM methods analyzed along the paper, for further details and additionalmethods we refer to [1, 2, 19, 26, 45, 47].
Acknowledgments
The authors are grateful to the unknown referees for having reported a series of inaccuracies in thefirst version of the manuscript and for the indications that led to this improved revision.
References [1]
G. Akrivis , Implicit–explicit multistep methods for nonlinear parabolic equations , Math. Comput., 82 (2012) 45–68.[2]
G. Akrivis, M. Crouzeix and C. Makridakis , Implicit–explicit multistep methods for quasilinear parabolic equations ,Numer. Math., 82 (1999) 521–541.[3]
U. M. Ascher, S. J. Ruuth, R. J. Spiteri , Implicit-Explicit Runge-Kutta Methods for Time-Dependent PartialDifferential Equations , Appl. Numer. Math, 1997, V. 25, pages 151–167.[4]
U. M. Ascher, S. J. Ruuth, B. T. R. Wetton , Implicit-Explicit methods for time-dependent partial differentialequations , SIAM J. Num. Anal. 32, (1995), 797–823.[5]
C. Bardos, C.D. Levermore, F. Golse , Fluid dynamic limit of kinetic equations II: Convergence proofs for theBoltzmann equations , Comm. Pure App. Math., 46, (1993), 667-753.[6]
G.I. Barenblatt , Scaling,
Cambridge Texts in Applied Mathematics , Cambridge Univ. Press, Cambridge, 2003.[7]
S. Boscarino, L. Pareschi, G. Russo , Implicit-Explicit Runge-Kutta schemes for hyperbolic systems and kineticequations in the diffusion limit , SIAM J. Sci. Comput., 35 (2011), A22-A51.[8]
S. Boscarino, L. Pareschi, G. Russo , A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscalerelaxation , SIAM J. Num. Anal., 55, No. 4, (2017), 2085-2109.[9]
S. Boscarino, G. Russo , On a class of uniformly accurate IMEX Runge-Kutta schemes and applications to hyperbolicsystems with relaxation , SIAM J. Sci. Comput. Vol. 31, No. 3, 1926-1945.[10]
S. Boscarino, G. Russo , Flux-Explicit IMEX Runge-Kutta schemes for hyperbolic to parabolic relaxation problems ,SIAM J. Num. Anal., Vol. 51, No. 1, 163–190.[11]
R. E. Caflisch, S. Jin, G. Russo , Uniformly accurate schemes for hyperbolic systems with relaxation , SIAM J. Num.Anal., 34, 1, (1997), 246-281.[12]
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.[13]
C. Cercignani, R. Illner, M. Pulvirenti , The mathematical theory of dilute gases,
Applied Mathematical Sciences ,106, Springer-Verlag, New York, (1994).
Table 6.
Examples of IMEX-LM methods. The numbers between brackets denoterespectively: s number of steps, p order of the method. For BDF methods the numberof steps is equivalent to the order of the method. IMEX Explicit ImplicitSG(3,2) a T ( − , , − , , c T (0 , , , , b T ( , , , , c − a T ( − , , , , c T (0 , , , , b T ( , − , , , c − TVB(3,3) a T ( − , , − , , c T ( − , − , , , b T ( , − , , , c − / a T ( − , , − , , c T (0 , , , , b T ( , − , , , c − TVB(4,4) a T ( − , , − , , c T ( − , − , , − , b T ( , − , , − , c − BDF4 a T ( − , , − , , c T (0 , , , , b T ( , − , , − , c − TVB(5,5) a T ( − , , − , , − ) c T ( − , , , − , ) b T ( , − , , − , ) c − BDF5 a T ( − , , − , , − ) c T (0 , , , , b T ( , − , , − , ) c − [14] G. Q. Chen, C. D. Levermore, T. P. Liu , Hyperbolic conservation laws with stiff relaxation terms and entropy ,Comm. Pure and App. Math., XLVII, (1994), 787 –830.[15]
C. Cercignani , The Boltzmann equation and its applications, Springer Varlag, New York, (1988).[16]
G. Dimarco, L. Pareschi , Numerical methods for kinetic equations , Acta Numerica 23, 369–520, 2014.[17]
G. Dimarco, L. Pareschi , Exponential Runge-Kutta methods for stiff kinetic equations , SIAM J. Num. Anal. (2011),2057–2077.[18]
G. Dimarco, L. Pareschi , Asymptotic-preserving IMEX Runge-Kutta methods for nonlinear kinetic equations , SIAMJ. Num. Anal. (2013), 1064–1087.[19]
G. Dimarco, L. Pareschi , Implicit-Explicit linear multistep methods for stiff kinetic equations . SIAM J. Num. Anal.55(2), (2017), 664–690.[20]
F. Filbet, S. Jin , A class of asymptotic preserving schemes for kinetic equations and related problems with stiffsources , J. Comp. Phys., 229 (2010), 7625–7648.[21]
J. Frank, W. Hundsdorfer, J. G. Verwer , On the stability of implicit-explicit linear multistep methods . Appl. Num.Math. 25, (1997), 193–205.[22]
E. Gabetta, L. Pareschi, M. Ronconi , Central schemes for hydrodynamical limits of discrete-velocity kinetic equa-tions , Transp. Theo. Stat. Phys. 29, 3-5, (2000), 465–477.[23]
E. Gabetta, B. Perthame , Scaling limits of the Ruijgrook-Wu model of the Boltzmann equation , Proceedings of theInternational Conference on Nonlinear Equations and Applications, Bangalore 19-23 August 1996. Springer-Verlag.[24]
E. Gabetta, L. Pareschi, G. Toscani , Relaxation schemes for nonlinear kinetic equations , SIAM J. Numer. Anal.34, (1997), no. 6, 2168–2194.[25]
J. Haack, S. Jin, J.-G. Liu , An All-Speed Asymptotic-Preserving Method for the Isentropic Euler and Navier-StokesEquations , Commun. Comput. Phys. 12, (2012), 955–980.[26]
W. Hundsdorfer, S. J. Ruuth , IMEX extensions of linear multistep methods with general monotonicity and bound-edness properties , J. Comp. Phys. 225, (2007), 2016–2042.[27]
W. Hundsdorfer, S. J. Ruuth, R. J. Spiteri , Monotonicity-preserving linear multistep methods , SIAM J. Numer.Anal. 41, (2003) 605–623.[28]
A. Klar , An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit , SIAM J. Numer.Anal. 35(3), (1998), 1073–1094.[29]
A. Klar , An Asymptotic Preserving Numerical Scheme for Kinetic Equations in the Low Mach Number Limit , SIAMJ. Numer. Anal. 36(5), (1999), 1507–1527.
MEX MULTISTEP METHODS FOR HYPERBOLIC SYSTEM WITH RELAXATION 29 [30]
S. Jin , Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations , SIAM J. Sci. Comput.,21(2), (1999), 441–454.[31]
S. Jin, C. D. Levermore , Numerical schemes for hyperbolic conservation laws with stiff relaxation terms , J. Comput.Phys., 126 (1996), 449–467.[32]
S. Jin, L. Pareschi , Asymptotic-preserving (AP) schemes for multiscale kinetic equations: a unified approach , Hy-perbolic Problems: Theory, Numerics, Applications, Vol. 141, International Series of Numerical Mathematics, (2001),573–582.[33]
S. Jin, L. Pareschi, G. Toscani , Diffusive relaxation schemes for discrete-velocity kinetic equations , SIAM J. Numer.Anal. 35, (1998), 2405–2439.[34]
S. Jin, L. Pareschi, G. Toscani , Uniformly accurate diffusive relaxation schemes for transport equations , SIAM J.Numer. Anal. 38, (2000), 913–936.[35]
S. Jin, Z. Xin , The relaxation schemes for systems of conservation laws in arbitrary space dimension , Comm. Pureand Appl. Math. 48, (1995), 235–276.[36]
P.L. Lions, G. Toscani , Diffusive limit for finite velocities Boltzmann kinetic models , Rev. Mat. Iberoamericana 13,no. 3, (1997), 473–513.[37]
M. Lemou, L. Mieussens , A new asymptotic preserving scheme based on micro-macro formulation for linear kineticequations in the diffusion limit , SIAM J. Sci. Comput., 31 (2008), 334–368.[38]
T. P. Liu , Hyperbolic conservation laws with relaxation , Comm. Math. Phys., 108, 153 (1987).[39]
P. Marcati, B. Rubino , Hyperbolic to parabolic relaxation theory for quasi-linear first order systems , Journal ofDifferential Equations 162, (2000), 359–399.[40]
G. Naldi, L. Pareschi , Numerical schemes for kinetic equations in diffusive regimes , Appl. Math. Letters, 11, 2,(1998), 29–35.[41]
G. Naldi, L. Pareschi , Numerical schemes for hyperbolic systems of conservation laws with stiff diffusive relaxation .SIAM J. Numer. Anal. 37 (2000), no. 4, 1246–1270.[42]
R. Natalini , Convergence to equilibrium for the relaxation approximations of conservation laws . Commun. Pure Appl.Math., 49:795–823, 1996.[43]
L. Pareschi, G. Russo , Implicit-Explicit Runge-Kutta schemes and applications to hyperbolic systems with relaxations ,J. Sci. Comput., 25 (2005), 129–155[44]
R. B. Pember , Numerical methods for hyperbolic conservation laws with stiff relaxation, I. Spurious solutions , SIAMJ. Appl. Math., 53, (1993), 1293–1330.[45]
R. R. Rosales, B. Seibold, D. Shirokoff, D. Zhou , Unconditional stability for multistep IMEX schemes: Theory ,SIAM J. Numer. Anal. 55(5), (2017), 2336–2360.[46]
W. Ruijgrook, T. T. Wu , A completely solvable model of the nonlinear Boltzmann equation , Physica, 113 A, (1982),401–416.[47]
B. Seibold, D. Shirokoff, D. Zhou , Unconditional stability for multistep IMEX schemes: Practice , J. Comp. Phys.376, (2019), 295–321.[48]
C. W. Shu , Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservationlaws , Advanced numerical approximation of nonlinear hyperbolic equations, Springer, (1998), 325–432.[49]