Implicit-Explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit

# Implicit-Explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit

S. Boscarino Mathematics and Computer Science Department, University of Catania, Italy (boscarino@dmi.unict.it).    L. Pareschi Mathematics Department, University of Ferrara, Italy (lorenzo.pareschi@unife.it).    G.Russo Mathematics and Computer Science Department, University of Catania, Italy (russo@dmi.unict.it).
###### Abstract

We consider Implicit-Explicit (IMEX) Runge-Kutta (R-K) schemes for hyperbolic systems with stiff relaxation in the so-called diffusion limit. In such regime the system relaxes towards a convection-diffusion equation. The first objective of the paper is to show that traditional partitioned IMEX R-K schemes will relax to an explicit scheme for the limit equation with no need of modification of the original system. Of course the explicit scheme obtained in the limit suffers from the classical parabolic stability restriction on the time step. The main goal of the paper is to present an approach, based on IMEX R-K schemes, that in the diffusion limit relaxes to an IMEX R-K scheme for the convection-diffusion equation, in which the diffusion is treated implicitly. This is achieved by an original reformulation of the problem, and subsequent application of IMEX R-K schemes to it. An analysis on such schemes to the reformulated problem shows that the schemes reduce to IMEX R-K schemes for the limit equation, under the same conditions derived for hyperbolic relaxation [8]. Several numerical examples including neutron transport equations confirm the theoretical analysis.

Key words. IMEX Runge-Kutta methods, hyperbolic conservation laws with sources, transport equations, diffusion equations, stiff systems.

AMS subject classification. 65C20, 65M06, 76D05, 82C40.

## 1 Introduction

The development of numerical methods to solve hyperbolic systems in diffusive regimes has been a very active area of research in the last years (see for example [25, 28, 31, 34]).

Classical fields of applications involve diffusion in neutron transport [9, 16, 19, 26, 32], drift-diffusion limit in semiconductors [27, 30] and incompressible Navier-Stokes limits in rarefied gas dynamic [29]. A strictly related field of research concerns the construction of schemes for the compressible Navier-Stokes limit (see [3] and the references therein). In such physical problems, the scaling parameter (mean free path) may differ in several orders of magnitude from the rarefied regimes to the diffusive regimes, and it is desirable to develop a class of robust numerical schemes that can work uniformly with respect to this parameter.

A prototype hyperbolic system of conservation laws with diffusive relaxation that we will use to illustrate the subsequent theory is the following, [25, 27, 34]

 ut+vx=0,vt+1ε2p(u)x=−1ε2(v−q(u)), \hb@xt@.01(1.1)

where . System (LABEL:I3) is hyperbolic with two distinct real characteristics speeds .

In the small relaxation limit, , the behavior of the solution to (LABEL:I3) is, at least formally, governed by the convection-diffusion equation

 ut+q(u)x=p(u)xx,v=q(u)−p(u)x. \hb@xt@.01(1.2)

The so called subcharacteristic condition [15] for system (LABEL:I3) becomes

 |q′(u)|2

and is naturally satisfied in the limit .

About the boundary conditions for system (LABEL:I3), we have to specify the domain in which we solve the problem. In a finite domain , one can use periodic boundary conditions, or can assign one independent condition at each boundary, since the two characteristic velocities have opposite sign.

In practice, we shall assign two conditions at each boundary, one independent and the other compatible with the equations. Furthermore, we shall choose boundary conditions which are compatible to the limit solution as . For example for system (LABEL:I3), if we set at and , compatibility with system (LABEL:I3) requires

 ε2q′(u)vx=p′(u)ux,   x=a,b.

Such condition becomes in the limit . In the sections on the numerical tests we shall specify the boundary conditions we use in each test.

In general, numerical approaches that work for hyperbolic system with stiff relaxation terms do not apply directly in the diffusive scaling since in these systems we have the presence of multiple time-scales.

In fact, together with the stiff relaxation term we have a stiff convection term that contributes to the asymptotic diffusive behavior. Then special care must be taken to ensure that the schemes possess the correct zero-relaxation limit, in the sense that the numerical scheme applied to system (LABEL:I3) should be a consistent and stable scheme for the limit system (LABEL:I4) as the parameter approaches zero independently of the discretization parameters. A notion usually referred to as asymptotic preservation. For a nice survey on asymptotic preserving scheme for various kinds of systems see, for example, the review paper by Shi Jin ([23]). Furthermore, a different approach to the derivation of asymptotic preserving schemes is described in the review by Pierre Degond [17]. In the case of Boltzmann kinetic equations we also refer to the recent review by two of the authors [37].

IMEX Runge-Kutta (R-K) schemes [1, 5, 6, 11, 36] represent a powerful tool for the time discretization of such stiff systems. Unfortunately, since the characteristic speed of the hyperbolic part is of order , standard IMEX R-K schemes developed for hyperbolic systems with stiff relaxation [36, 8] become useless in such parabolic scaling, because the CFL condition would require . Of course, in the diffusive regime where , this is very restrictive since for an explicit method a parabolic condition should suffice.

Most previous works on asymptotic preserving schemes for hyperbolic systems and kinetic equations with diffusive relaxation focus on schemes which in the limit of infinite stiffness become consistent explicit schemes for the diffusive limit equation [25, 28, 31, 32, 34]. Such schemes have been derived by splitting the stiff hyperbolic part into an explicit (non-stiff) term, and an implicit (stiff) term. Here we show that by applying partitioned IMEX R-K schemes, in which the stiffness is associated with the variable and not with the operator, one obtains IMEX R-K schemes that naturally relax to the explicit scheme applied to the limit convection-diffusion equation. All these explicit schemes clearly suffer from the usual stability restriction .

In this paper we present a general methodology to overcome such stability restriction which applies to a broad class of problems. The idea is to reformulate problem (LABEL:I3) by properly combining the limiting diffusion flux with the convective flux. This allows to construct a class of IMEX R-K schemes that work with high order accuracy in time and that, in the diffusion limit (i.e. when ), originate an IMEX method for the limiting convection-diffusion equation (LABEL:I4). Other reformulations whose goal is to obtain asymptotic-preserving methods have been proposed in [25, 24]. Schemes that avoid such time step restriction and provide fully implicit solvers in the case of transport equations have been analyzed in [9].

Our new approach allows a hyperbolic CFL condition independent of when applied to (LABEL:I3) in all regimes. The aim of this paper is to derive and analyze different types of IMEX R-K schemes when applied to the reformulated problem in the stiff regime ().

The rest of the paper is organized as follows. The next section is devoted to partitioned IMEX R-K schemes. It is shown that they relax to the explicit scheme applied to the convection diffusion limit. In Section LABEL:Sec3 the new approach is introduced and analyzed. In particular, following [8], we prove that under suitable assumptions the IMEX R-K schemes are consistent with the diffusion limit. The analysis is based on a power series expansion in of the exact and numerical solution. It is shown that, at lowest order in , the model system is a set of differential-algebraic equations of index 1, i.e. it can be transformed into a set of ordinary differential equations by one time differentiation. Compatibility between exact and numerical solution at different orders in introduces additional order conditions on the coefficients of the IMEX schemes. To the lowest order such conditions are referred as index 1 order conditions.

After a short section on space discretization obtained by conservative finite difference schemes, in Section LABEL:Sec5 we report several numerical examples and tests. In Section LABEL:Sec6 we consider one-dimensional neutron transport equation and present several numerical results and comparison with schemes available in the literature. Additional technical material is given in the separate Appendices.

## 2 Partitioned IMEX R-K schemes

The first observation is that in system (LABEL:I3) the stiffness is naturally associated to the variable rather then to some operator. The system has the structure of a singular perturbation problem [21], and it can be treated by a partitioned R-K scheme in which the first equation is treated explicitly and the second implicitly

 ut=−vx,(Explicit)% ε2vt=−(p(u)x+v−q(u)).(Implicit) \hb@xt@.01(2.1)

This approach has been used, for example, in [28, 35].

By using a method of lines approach (MOL), we discretize system (LABEL:I33) in space by a uniform mesh and , . We obtain a large system of ODE’s

 Ut=F(V),ε2Vt=G(U)−V, \hb@xt@.01(2.2)

with and , where and . Here and (with a slight abuse of notation) denote the discretization of the convective terms , , while represents the discretization of the term .

As we shall see, if the implicit scheme is -stable [21], in the limit the IMEX R-K scheme will relaxe to the explicit scheme applied to the limit equation

 Ut=:^F(U), \hb@xt@.01(2.3)

where .

For example, implicit-explicit Euler scheme applied to system (LABEL:MLa) gives

 Un+1 = Un+ΔtF(Vn) Vn+1 = ε2Vn+ΔtG(Un+1)ε2+Δt,

where we have discretized the interval of integration by a time mesh and . As , and therefore .

In the case and , the method would relax to explicit Euler scheme applied to the diffusion equation, thus suffering the usual parabolic CFL stability condition . This approach will be denoted as partitioned approach.

### 2.1 Classification of IMEX R-K schemes

IMEX R-K schemes have been widely used in the literature to treat problems that contain both stiff and non stiff terms [1, 11, 36]. The stiff terms are treated implicitly, while the non stiff terms are treated explicitly, thus lowering the computational complexity of the scheme.

Usually such a scheme is characterized by the matrices , and the vectors , , and can be represented by a double in the usual Butcher notation

 ~c~A  ~bT,cA  bT.

The coefficients and are used if the right hand side depends explicitly on time. We assume that they satisfy the usual relation

 ~ci=i−1∑j=1~aij,ci=i∑j=1aij. \hb@xt@.01(2.4)

Matrix is lower triangular with zero diagonal, while matrix is lower triangular, i.e. the implicit scheme is a diagonally implicit Runge-Kutta (DIRK). This choice guarantees that the term in (LABEL:MLa) is always explicitly evaluated.

IMEX R-K schemes presented in the literature can be classified in two main different types characterized by the structure of the matrix of the implicit scheme.

###### Definition 2.1

We call an IMEX R-K method of type A (see [36]) if the matrix is invertible.

###### Definition 2.2

We call an IMEX R-K method of type CK (see [11]) if the matrix can be written as

 A=(00a^A)

with and the submatrix invertible. In the special case the scheme is said to be of type ARS (see [1]).

We note that schemes CK are very attractive because they allow some simplifying assumptions, that make order conditions easier to treat, therefore permitting the construction of higher order IMEX R-K schemes. On the other hand, schemes of type A are more amenable to a theoretical analysis, since the matrix of the implicit scheme is invertible. This is why we start our analysis with the latter schemes.

### 2.2 Analysis of IMEX schemes for the partitioned approach

Now as an example we perform the analysis of type A scheme when applied to system (LABEL:I33). A similar analysis is possible also for CK schemes. We will restrict our analysis to the limit case .

Applying an IMEX R-K scheme to system (LABEL:MLa) we obtain

 Un+1 = Un+Δts∑k=1~bkF(Vk), \hb@xt@.01(2.5) ε2Vn+1 = ε2Vn+Δts∑k=1bk(G(Uk)−Vk),

for the numerical solution and

 Uk = \hb@xt@.01(2.6) ε2Vk = ε2Vn+Δtk∑j=1akj(g(Uj)−Vj),

for the internal stages.

By Definition LABEL:Def:A and invertible we obtain from the second equation in (LABEL:typeAPart)

 Δt(G(Uk)−Vk)=ε2k∑j=1ωkj(Vj−Vn). \hb@xt@.01(2.7)

From now on, are the elements of the inverse matrix . Now inserting (LABEL:gPar) into the numerical solution and setting , we get

 Vn+1 = R(∞)Vn+Δts∑k=1bkωkjVj. \hb@xt@.01(2.8)

Here we denoted by

 R(∞)=1−bTA−11=limz→∞R(z),

where is the stability function of the implicit scheme defined by (see [21], Sect. IV.3)

 R(z)=1+zbT(I−zA)−11, \hb@xt@.01(2.9)

with and .

By (LABEL:gPar) we get when . This yields that , and we obtain

 Un+1 = Un+Δts∑k=1~bk^F(Uk)

with

 Uk = Un+Δtk−1∑j=1~akj^F(Uj)

internal stages and . This represents the explicit scheme of the starting IMEX R-K one of type A applied to the limit equation (LABEL:MLalimit) obtained by (LABEL:MLa) when . As particular case, if and , this is the explicit scheme applied to the limit diffusion equation under the usual parabolic stability restriction ().

## 3 Overcoming parabolic stiffness

In order to overcome such stability restriction, we reformulate system (LABEL:I3) as the equivalent system

 ut=−(v+μp(u)x)x+μp(u)xx,ε2vt=−p(u)x−v+q(u), \hb@xt@.01(3.1)

where the term has been added and subtracted to the first equation in (LABEL:I3). Here is a free parameter such that . The idea is that, since the quantity is close to as , the first term on the right hand side can be treated explicitly in the first equation, while the term will be treated implicitly. This can be done naturally by using an Implicit-Explicit approach, as we will explain later. Let us point out that the choice , as shown in Appendix LABEL:A1 for a first order implicit-explicit scheme, guarantees the largest stability region of the method.

Next we will study the behavior of the different IMEX R-K schemes when applied to system (LABEL:I3b) in the diffusion limit. In particular we will show that such schemes reduce to the same IMEX R-K schemes for the limit equation and no parabolic stability restriction on the time step appears in the diffusive limit.

### 3.1 The new approach

System (LABEL:I3b) can be written in the form

 u′=f1(u,v)+f2(u),ε2v′=g(u,v) \hb@xt@.01(3.2)

where the primes denote the time derivatives and

 f1(u,v)=−(v+μp(u)x)x,f2(u)=μp(u)xx,
 g(u,v)=−p(u)x−v+q(u).

Notice that, throughout this paper, (and therefore ), depends only algebraically on , while it may contain differential operators acting on .

Now we apply an IMEX-RK scheme to system (LABEL:SPP) where is evaluated explicitly and implicitly. Note that if is evaluated explicitly then by cancelation the IMEX-RK scheme will reduce to the typology of asymptotic preserving methods studied in [7, 35].

In the limit from (LABEL:SPP) we obtain a differential algebraic system (DAE)

 u′=f1(u,v)+f2(u),0=g(u,v). \hb@xt@.01(3.3)

In order to guarantee the solvability of system (LABEL:DA1) we assume that the Jacobian matrix is invertible, and then the DAE is said to be of index one 111The index of a DAE is the number of times one has to differentiate the function to obtain a system of ODE’s. For example, differentiating the function , one obtains . If is invertible, system (LABEL:DA1) can be written as , .. Note that if has a bounded inverse in a neighborhood of the exact solution, we can use the inverse function theorem in order to write

 v(t)=G(u(t))

for some which inserted into gives . From now on we always assume that this is the case. Then, as system (LABEL:SPP) reduces to

 u′ = ^f1(u)+f2(u), \hb@xt@.01(3.4)

where and . This approach will be denoted BPR approach.

First order implicit-explicit Euler scheme that uses this approach is reported in Appendix LABEL:A1, where a stability analysis is performed. In particular it is shown that as , the parabolic restriction on time step disappears.

In the sequel we restrict our analysis to the limit case where the main goal is to capture the diffusive limit.

### 3.2 Analysis of TYPE A IMEX schemes

Applying an IMEX R-K scheme of type A to system (LABEL:SPP) we obtain

 un+1 = un+Δts∑k=1~bkf1(Uk,Vk)+Δts∑k=1bkf2(Uk) \hb@xt@.01(3.5) ε2vn+1 = ε2vn+Δts∑k=1bkg(Uk,Vk),

for the numerical solution and

 Uk = un+Δtk−1∑j=1~akjf1(Uj,Vj)+Δtk∑j=1akjf2(Uj) \hb@xt@.01(3.6) ε2Vk = ε2vn+Δtk∑j=1akjg(Uj,Vj),

for the internal stages (notice a slight changes of notation with respect to Section LABEL:Sec2).

Starting from (LABEL:typeA) and (LABEL:stypeA), by Definition (LABEL:Def:A) and invertible, we obtain from the second equation in (LABEL:stypeA)

 Δtg(Uk,Vk)=ε2k∑j=1ωkj(Vj−vn),

Inserting this into the numerical solution we make independent of and setting , we get

 un+1 = un+Δts∑k=1~bk^f1(Uk)+Δts∑k=1bkf2(Uk) vn+1 = R(∞)vn+Δts∑k=1bkωkjVj,

with , and stage values

 Uk = un+Δtk−1∑j=1~akj^f1(Uj)+Δtk∑j=1akjf2(Uj) \hb@xt@.01(3.8) 0 = g(Uk,Vk).

for . The latter equality implies , .

Note that if we require that the implicit part of the scheme is stiffly accurate, i.e. if

 bTA−1=eTs,

where , then by (LABEL:Sfunc)

 R(∞)=1−bTA−11=1−eTs1=1−1=0.

This implies that if the implicit scheme is -stable and stiffly accurate it is also -stable and .

It is interesting to note that, if we consider system (LABEL:I3b) with , when we get a purely diffusive system which means that the term in (LABEL:SPP) disappears. Therefore, by BPR approach, the IMEX R-K scheme of type A in the limit becomes a stiffly accurate DIRK scheme and hence no stability restriction on the time step is required in the diffusive limit, i.e. we got an unconditionally stable method. Another advantage of this new approach is the following. Usually the numerical solution in (LABEL:index1A) in the case will not lie on the manifold since is not necessarily zero. But this approach guarantees that in the limit we obtain a stiffly accurate implicit scheme and hence , implying .

In the general case of systems for which , it is and, by using the BPR approach, in the limit case we obtain an IMEX R-K scheme with a non vanishing explicit term in which the diffusion term is treated implicitly and a classical CFL hyperbolic condition for the time step is required. In general even if all stage values lie on the manifold, (see the second equation in (LABEL:index1bisA)). However, if the explicit scheme has the property that , and the implicit scheme is stiffly accurate, then, in the limit as the numerical solutions are projected on the manifold , because .

From the above discussion it is clear that the property is crucial if we want that the numerical solution is projected to the limit manifold as . We emphasize that there is a class of -stage explicit R-K methods for which ; such methods are called First Same As Last (FSAL), and their coefficients satisfy for and . They have the advantage of requiring function evaluations for each step (see [22] for details). FSAL methods are often used in the contest of embedded methods, such as the popular Dormond-Prince method (DOPRI) [18], on which MATLAB routine ode45 is based on.

From the arguments above, in order to capture the limit as , it is important that the implicit part on an IMEX R-K is stiffly accurate and the explicit part is FSAL. This motivates the following

###### Definition 3.1

We say that a IMEX R-K scheme is globally stiffly accurate if and , with , and , i.e. the numerical solution is identical to the last internal stage value of the scheme.

From (LABEL:typeA) and (LABEL:stypeA) we observe that if an IMEX R-K is globally stiffly accurate, then , , and therefore .

#### General remarks for type A

• It is worth mentioning some important aspects about type A schemes. First of all, in [5] Boscarino emphasized that an important ingredient for the IMEX R-K schemes of type A is for all . Such a choice provides a significant benefit for the differential component , i.e., an order reduction does not appear for this component. On the other hand, conditions

 eTs~A=~bT,eTsA=bT

imply which means that for a stiffly accurate IMEX R-K scheme it is , and therefore we expect to observe order reduction for the differential variable.

• It is impossible to construct a second order stiffly accurate IMEX R-K scheme of type A with internal stages. The proof is given in Appendix LABEL:A2. In practice, in order to satisfy all these order conditions we have to increase the number of the internal stages. In view of such difficulties, for type A schemes, we shall consider second order IMEX R-K schemes with and in order to avoid the order reduction, giving up to the FSAL property of the explicit scheme (and with that the global stiff accuracy of the IMEX scheme). An example is the scheme SSP(3,3,2) in Appendix LABEL:A4. In this case, if , it is as .

• Formulation (LABEL:SPP) in the limit case yields the index-1 DAE. Then using the same technique adopted in [6], we can derive additional order conditions, called algebraic conditions, that guarantee the correct behavior of the numerical solution in the limit and maintain the accuracy in time of the scheme. If the implicit scheme is stiffly accurate, such conditions becomes, to various order of accuracy,

 ~cs=1, (second order)eTs~A~c=1/2, (third order) \hb@xt@.01(3.9)

where . If the IMEX schemes is globally stiffly accurate, then (LABEL:oc_index1) are automatically satisfied, since .

• Finally we observe that, in order to construct an order IMEX R-K of type A and to maintain accuracy we have to increase the number of the classical order conditions too. Usually several simplifying assumptions (see [6], [8], [21] for details) could help to reduce the number of such conditions, but, higher orders type A schemes are more complicated to construct than CK or ARS schemes because of additional order conditions (see [8]) due to the fact that .

### 3.3 Analysis of TYPE CK schemes

Similar considerations about BPR approach, explained for the IMEX R-K scheme of type A in the limit case , can be reproposed here for the type CK when applied to the system (LABEL:DA1), with slightly modifications. Of course, if we consider the general system (LABEL:I3b) we obtain again an IMEX R-K scheme of type CK in the diffusion limit, i.e. , where the diffusion term is treated implicitly and a CFL hyperbolic condition for the time step is required.

Indeed, we consider an IMEX R-K schemes of type CK where, by Definition LABEL:Def:CK, we assume that the submatrix is invertible and . The Butcher tableaux of a CK scheme takes the form

 00  0^ca  ^Ab1  ^bT

with and . In order to simplify the analysis we consider that the implicit part of the scheme is stiffly accurate. Under this circumstance it is easy to prove that

 b1+^bTα=0, \hb@xt@.01(3.10)

where (see [5] for details).

Then, considering a scheme of the type CK, the second equation in (LABEL:stypeA) becomes

 ε2Vk=ε2vn+Δtak1g(un,vn)+Δtk∑j=2akjg(Uj,Vj). \hb@xt@.01(3.11)

with .

Now multiplying by , where are the elements of the inverse of , and summing on , we obtain

 Δtg(Uk,Vk)=ε2s∑j=2^ωkj(Vj−vn)+Δtαkg(un,vn),  for  k=2,…,s

where

 s∑l=2^ωklalj=δkj,−s∑l=2^ωkjaj1=αk.

Inserting the expression into the second equation in (LABEL:typeA) we obtain

 ε2vn+1 = ε2R(∞)vn+ε2Δts∑k=2bkωkjVj+Δt(b1+s∑k=2bkαk)g(un,vn) \hb@xt@.01(3.12)

Then by (LABEL:prorCK) the last term in the second equation in (LABEL:index1) drops and in the limit case for we can write

 vn+1 = R(∞)vn+Δts∑k=2bkωkjVj,

with

 g(Uk,Vk) = αkg(un,vn),  for  k=2,…,s. \hb@xt@.01(3.13)

Note that, for IMEX R-K schemes of type CK, the stability function of the implicit part of the scheme takes the form

 R(z) = 1+z(b1+^bT(I−z^A)−1(1s−1+za)) = (b1−^bT^A−1a)z+(1−^bT^A−11s−1+^bT^A−2a)+O(1z).

We obtained this result, by applying one step of the implicit part of the scheme to the test problem , , with and .

Thus, the only stiffly accurate condition, i.e. is not enough to guarantee that and then an additional condition is required for the implicit part of the scheme, (for details see [8]). This is expressed by the following

###### Proposition 3.2

If

 −^eTs−1^A−1a=∑j≥2^ωsjaj1=0, \hb@xt@.01(3.15)

then , where .

Proof. In fact, assuming invertible, we get and when , from (LABEL:STAB) we obtain , which is zeros if (LABEL:eq:L) is satisfied.
Note that the previous Lemma implies that and by (LABEL:pro2CK) with we obtain , then the last stage values lie on the manifold as . Now we observe that if the IMEX R-K scheme of type CK is globally stiffly accurate, we obtain from (LABEL:index1) and (LABEL:index1A) and and therefore with .

Since an IMEX R-K schemes of type ARS is a particular case of the type CK where the vector , then the same results hold true.

#### General remarks for type CK

• IMEX CK schemes [11] are attractive because of their good properties. The implicit part of this scheme is singly diagonally implicit Runge-Kutta (SDIRK) with for and differs from the classical SDIRK one because . In [11] such implicit schemes are called explicit singly diagonally implicit (ESDIRK). A consequence to set is the possibility to guarantee stage-order higher than the in the case of SDIRK, for which . Moreover here we consider schemes that are stiffly accurate according to Definition LABEL:1def. Such schemes will project the solution on the manyfold in the limit of infinite stiffness. For these schemes , so one of the so-called simplifying conditions cannot be applied [8]. Here we require that for all ; this choice will reduce the number of coupled order conditions.

## 4 IMEX-Finite Difference schemes

When constructing numerical schemes, one has also to take a great care in order to avoid spurious numerical oscillations arising near discontinuities of the solution. This is avoided by a suitable choice of space discretization. To this aim it is necessary to use non-oscillatory interpolating algorithms, in order to prevent the onset of spurious oscillations (like ENO and WENO methods), see [40]. Moreover the choice of the space discretization may be relevant for a correct treatment of the boundary conditions.

In this section we emphasize some requirements about the space discretization of the system (LABEL:I3b). We remark that the dissipative nature of upwind schemes [34, 35] depends essentially on the fact that the characteristic speeds of the hyperbolic part are proportional to . On the other hand central differences schemes avoid excessive dissipation but when is not small or when the limiting equations contain advection terms may lead to unstable discretizations. In order to overcome these well-known facts and to have the correct asymptotic behavior we fix some general requirements for the space discretization.

1. Correct diffusion limit. Let us consider system (LABEL:I3b) with . In the limit case the therm from the second equation. If we want that also in the first equation, we need to use the same space discretization for the term and require that .

2. Compact stencil. Among the advantages of our approach there is the possibility to have a scheme with a compact stencil in the diffusion limit . This property is satisfied if point 1) is satisfied and we use a suitable discretization for the second order derivative that characterize the diffusion limit.

3. Shock capturing. The schemes when should be based on shock capturing high order fluxes for the convection part. This is necessary not only for large values of but also when we consider convection-diffusion type limit equations with small diffusion. The high order fluxes are then necessary for all space derivatives except for the second order term on the right-hand side.

4. Avoid solving nonlinear algebraic equations. In order to achieve this the implicit space derivative in the second equation must be evaluated using only nodal values of which can be obtained from the solution of the first equation.

The above properties are satisfied for example using high accuracy in space obtained by finite difference discretization with Weighted-Essentially Non Oscillatory (WENO) reconstruction, [40].

System (LABEL:I3b) may be written in the form

 ut+(v+μp(u)x)x = μp(u)xx, \hb@xt@.01(4.1) vt = 1ε2(q(u)−(v+p(u)x)).

with introduced in Section 2. The terms on the right-hand side will be treated implicitly. For large value of the explicit flux is just , while for small values of it is . Here we describe a finite difference WENO scheme for a system of the form

 Ut+F(U)x=R(U),

and apply it to the system (LABEL:hyper) with

 F(U)=(v+μp(u)x,0)T, R(U)=(μp(u)xx,1ε2(q(u)−(v+p(u)x))).

As and , the system becomes

 ut+vx = 0, vt = 0

and the characteristic speed of the system is (twice). As and , and the system relaxes to the equation

 ut+q(u)x=p(u)xx

and the characteristic speed of the left hand side is given by .

Conservative finite difference for system (LABEL:hyper) are written as follows, [36]

 dUjdt = −^Fj+12−^Fj−12Δx+G(Uj)

where is an approximation of the pointwise value of at grid nodes, and the numerical flux at cell edge is computed as follows

 ^Fj+12=^F+j(xj+12)+^F−j+1(xj+12).

The function and are suitable reconstructions defined, respectively, in cell and in cell . They are obtained as follows. First, we assume that the flux can be split into a positive and negative component

 F(U)=F+(U)+F−(U),

with , . The quantity are computed at cell center. Then are reconstructed from using high order essentially non oscillatory reconstruction, such as ENO or WENO, that allows pointwise reconstruction of a function from its cell averages, (see, e.g. [40] for details).

The flux may contain derivatives. For example the first equation in system (LABEL:hyper) contains . Such terms are computed by point-wise WENO reconstruction.

In all our examples we used the simple local Lax-Friedrix flux decomposition, i.e. , , , , where denotes the spectrum radius of matrix , and the max defining is taken for varying in a suitable range in a neighborhood of each cell. In our test case we chose , since as , and in our numerical test is either or , with ranging in , therefore .

We remark here that the choice of is based on the characteristic speeds of the limit convection-diffusion equation, while a more detailed analysis is needed to justify its use in intermediate regions, for which the characteristic speeds can be much higher, and the stabilization that compensates for the apparent violation of the hyperbolic CFL condition comes from the implicit treatment of the diffusion term.

Furthermore for large value of , (e.g., ), we want to avoid adding and subtracting terms which may cause loss of accuracy. For a semidiscrete scheme the function will depend also on the grid space . A simple choice for is given by

 μ=exp(−ε2/Δx)

which is what we used in all our numerical tests.

For the diffusion term we used the standard 2-nd order finite difference technique for second order time discretization, and the standard 4-th order finite difference technique where 3-rd order time discretization are used.

## 5 Numerical examples

In this section we test several second and third order IMEX R-K schemes presented in the literature that satisfy the algebraic order conditions (LABEL:oc_index1) and conditions in Definition 3. Usually, IMEX time discretization are identified by an acronym (e.g. the initials of the authors), and three numbers (, , ) denoting, respectively, the effective number of stages (in practice the number of function evaluations) of the explicit and implicit scheme and the classical order of accuracy.

Below we list the IMEX R-K schemes used in the numerical tests.

• SSP(3,3,2): derived by Pareschi, Russo [36], it is a second order IMER R-K of type A, the explicit part is strongly stability preserving, while the implicit part is stiffly accurate. In accordance with the proposition LABEL:prop in the Appendix LABEL:A2, this scheme is not globally stiffly accurate according to Definition LABEL:1def.

• ARS(2,2,2): derived by Asher, Ruuth, Spiteri [1], it is a second order scheme, the double Butcher tableau of this scheme is reproduced in Appendix LABEL:A4. Note that this scheme is globally stiffly accurate according to the Definition LABEL:1def and satisfies the additional order conditions (LABEL:oc_index1).

• ARS(4,4,3): derived by Asher, Ruuth, Spiteri [1], it is a third order scheme, the double Butcher tableau of this scheme is reproduced in Appendix LABEL:A4. Similarly to ARS(2,2,2) this scheme is globally stiffly accurate according to Definition LABEL:1def and satisfies the additional order conditions (LABEL:oc_index1).

• BPR(3,5,3) introduced in this paper is a third order IMEX R-K scheme of type CK and globally stiffly accurate according to Definition LABEL:1def. This scheme has internal stages, explicit stages and implicit stages. The additional order conditions (LABEL:oc_index1) are satisfied. This scheme is more efficient than ARS(4,4,3) for the explicit part, but less efficient for the implicit one. In many cases the computation of the explicit term is more expensive than the solution of the implicit step, thus resulting in an overall improvement in efficiency per time step. Furthermore, the coefficients of the explicit scheme are all non negative, which is an advantage for the stability properties of the scheme. We reproduced the coefficients of this scheme in Appendix LABEL:A4.

In all the computations presented in this paper we denote each scheme with an acronym indicating the IMEX scheme and the type of space discretization.

Space discretization is identified by a short name containing the order of accuracy in space; for example, WENO53 (or WENO32), see for details [40], means a fifth (or a third) order reconstruction which reduced to third (or second) order near singularities and CdS2 stands for second order central discretization scheme.

We remark here that all the analysis performed in the paper and the numerical tests are performed under the assumption that the initial data is well-prepared, which means that the initial condition lies in the limit manyfold as . If this condition is not satisfied, then a loss of accuracy is observed, unless some initial layer fix is adopted. Schemes of type A are more robust against this problem, as is described in [36].

### 5.1 Convergence test

In this section we investigate numerically the convergence rate of the second and third IMEX R-K schemes introduced before for a wide range of the parameter . To this aim we apply these schemes to simple prototype hyperbolic system (LABEL:I3b), with initial conditions chosen in such a way that the exact solutions is smooth and does not present a rapidly varying transient for small values of . This is achieved in practice by imposing that the initial condition satisfies the limit relation between and as .

Numerical convergence rate is calculated by the formula

 p=log2(EΔt1/EΔt2),

where and are the global errors computed with step , and . In the following tests we put and we choose .

For the first test we set and . Then we get

 ut=−vx−μuxx+μuxx,ε2vt=−ux−v, \hb@xt@.01(5.1)

that in the limit case, and leads to the linear diffusive problem

 ut=uxx,   u(x,0)=u0(x). \hb@xt@.01(5.2)

We use periodic boundary conditions with , and , so that is an exact solution of (LABEL:Limdiff). The final time is and .

The results are reported in Table 5.1 and 5.2 showing that the expected convergence rates are reached for the -component.

Next we set and consider the following system

 ut+vx = μuxx−μuxx ε2vt+ux = −(v−u),

where the limiting behavior is given by an advection-diffusion equation. We use periodic boundary conditions with the initial data , with and , on the spatial interval , at the final time anfd . As reference solution we use the truncated Fourier representation of the exact solution

 Uexa(x,t)=+∞∑k=−∞Uk(t)eikx,  Vexa(x,t)=+∞∑k=−∞Vk(t)eikx

with and satisfying

 ˙Uk=−ikVk,ε2˙Vk=−ikUk+Uk−Vk. \hb@xt@.01(5.3)

For each , system (LABEL:FC1) can be written as a constant coefficient homogeneous system which can be solved exactly. The results are given in Table 5.3 showing that again the expected convergence rates are reached for the -component by all schemes.