A unified IMEX Runge-Kutta approach for hyperbolic systems with multiscale relaxation

TL;DR

This paper develops a unified IMEX Runge-Kutta method for hyperbolic systems with multiscale relaxation, ensuring accurate asymptotic behavior across different regimes, validated by numerical examples.

Contribution

A novel IMEX Runge-Kutta scheme that maintains asymptotic preservation for multiscale hyperbolic systems with relaxation.

Findings

The method accurately captures the correct asymptotic limit regardless of scaling.

Numerical tests confirm the theoretical asymptotic preservation.

The approach improves efficiency over standard methods in stiff regimes.

Abstract

In this paper we consider the development of Implicit-Explicit (IMEX) Runge-Kutta schemes for hyperbolic systems with multiscale relaxation. In such systems the scaling depends on an additional parameter which modifies the nature of the asymptotic behavior which can be either hyperbolic or parabolic. Because of the multiple scalings, standard IMEX Runge-Kutta methods for hyperbolic systems with relaxation loose their efficiency and a different approach should be adopted to guarantee asymptotic preservation in stiff regimes. We show that the proposed approach is capable to capture the correct asymptotic limit of the system independently of the scaling used. Several numerical examples confirm our theoretical analysis.

Tables

Table 1. Table 1: Test 1. Converge rates for the density u 𝑢 u with ε = 10 − 6 𝜀 superscript 10 6 \varepsilon=10^{-6} .

Method	N	$L_{\infty }$ error	Order
ARS(1,1,1)	40	$6.4800e-03$	$--$
ARS(1,1,1)	80	$3.5082e-03$	$0.8853$
ARS(1,1,1)	160	$1.9203e-03$	$0.8694$
ARS(1,1,1)	320	$9.6447e-04$	$0.9935$
ARS(1,1,1)	640	$4.8457e-04$	$0.9930$
CK(2,2,2)	40	$1.4911e-04$	$--$
CK(2,2,2)	80	$3.9405e-05$	$1.9199$
CK(2,2,2)	160	$1.1356e-05$	$1.7949$
CK(2,2,2)	320	$2.8331e-06$	$2.0030$
CK(2,2,2)	640	$7.0874e-07$	$1.9991$
BPR(3,4,3)	40	$5.8318e-06$	$--$
BPR(3,4,3)	80	$7.8658e-07$	$2.8903$
BPR(3,4,3)	160	$1.2095e-07$	$2.7012$
BPR(3,4,3)	320	$1.5297e-08$	$2.9831$
BPR(3,4,3)	640	$1.9253e-09$	$2.9901$

Table 2. Table 2: Test 1. Converge rates for v 𝑣 v for the scheme CK(2,2,2) with ε = 10 − 6 𝜀 superscript 10 6 \varepsilon=10^{-6} .

Method	N	$L_{\infty }$ $v$-error	Order of $v$
CK(2,2,2)	40	$1.3156e-02$	$--$
CK(2,2,2)	80	$1.4897e-02$	$-0.1793$
CK(2,2,2)	160	$1.1520e-03$	$3.6928$
CK(2,2,2)	320	$1.2584e-03$	$-0.1274$
CK(2,2,2)	$640$	$1.2877e-03$	$-0.0332$

Table 3. Table 3: Test 1. Converge rates for u 𝑢 u and v 𝑣 v for the scheme ( 134 ) with ε = 10 − 6 𝜀 superscript 10 6 \varepsilon=10^{-6} .

Method	N	$L_{\infty }$ $u$-error	Order $u$	$L_{\infty }$ $v$-error	Order $v$
BPR(4,4,2)	40	$1.9129e-04$	$--$	$2.8704e-04$	$--$
BPR(4,4,2)	80	$4.9963e-05$	$1.9368$	$8.0261e-05$	$1.8385$
BPR(4,4,2)	160	$1.4374e-05$	$1.7974$	$2.0603e-05$	$1.9618$
BPR(4,4,2)	320	$3.5895e-06$	$2.0016$	$5.2702e-06$	$1.9669$
BPR(4,4,2)	640	$9.0120e-07$	$1.9939$	$1.4011e-06$	$1.9113$

Equations

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}u+\partial_{x}v=0,\\ \\ \displaystyle\partial_{t}v+\frac{1}{{\varepsilon}^{2\alpha}}\partial_{x}p(u)=-\frac{1}{{\varepsilon}^{1+\alpha}}\left(v-f(u)\right),\quad\alpha\in[0,1]\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{\tau}u+\partial_{\xi}V=0,\\ \\ \displaystyle\partial_{\tau}V+\partial_{\xi}p(u)=-\frac{1}{{\varepsilon}}\left(V-F(u)\right),\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle v=f(u)-{\varepsilon}^{1-\alpha}\partial_{x}p(u)+{\varepsilon}^{1+\alpha}f^{\prime}(u)^{2}\partial_{x}u+{\cal{O}}({\varepsilon}^{2}),\\ \\ \displaystyle\partial_{t}u+\partial_{x}f(u)={\varepsilon}^{1+\alpha}\partial_{x}\Bigg{[}\left(\frac{p^{\prime}(u)}{{\varepsilon}^{2\alpha}}-f^{\prime}(u)^{2}\right)\partial_{x}u\Bigg{]}+{\cal{O}}({\varepsilon}^{2}).\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle v=f(u),\\ \\ \displaystyle\partial_{t}u+\partial_{x}f(u)=0.\\ \end{array}\right.$$

$$f^{\prime }(u{)}^2<\frac{p^{\prime }(u)}{{\epsilon }^{2\alpha }},$$

$$\left\{\begin{array}[]{l}\displaystyle v=f(u)-\partial_{x}p(u),\\ \\ \displaystyle\partial_{t}u+\partial_{x}f(u)=\partial_{xx}p(u).\\ \end{array}\right.$$

$$\frac{u^{n+1}-u^n}{\Delta t}$$

$${\epsilon }^2\frac{v^{n+1}-v^n}{\Delta t}$$

$$v^{n+1}=\frac{{\epsilon }^2}{{\epsilon }^2+\Delta t}{v}^n-\frac{\Delta t}{{\epsilon }^2+\Delta t}\left(p(u^n{)}_x-f(u^{n+1})\right).$$

$$\frac{u^{n+1}-u^n}{\Delta t}+\frac{{\epsilon }^2}{{\epsilon }^2+\Delta t}{v}_x^{n-1}$$

$$\frac{u^{n+1}-u^n}{\Delta t}$$

$$v^{n+1}=f(u^{n+1})-p(u^n{)}_x$$

$$\frac{u^{n+1}-u^n}{\Delta t}$$

$${\epsilon }^2\frac{v^{n+1}-v^n}{\Delta t}$$

$$v^{n+1}=\frac{{\epsilon }^2}{{\epsilon }^2+\Delta t}{v}^n-\frac{\Delta t}{{\epsilon }^2+\Delta t}\left(p(u^{n+1}{)}_x-f(u^{n+1})\right),$$

$$\frac{u^{n+1}-u^n}{\Delta t}+\frac{{\epsilon }^2}{{\epsilon }^2+\Delta t}{v}_x^{n-1}$$

$$\frac{u^{n+1}-u^n}{\Delta t}$$

$$\frac{(u^{n+1}+\epsilon v^{n+1})-(u^n+\epsilon v^n)}{\Delta t}$$

$$\frac{(u^{n+1}-\epsilon v^{n+1})-(u^n-\epsilon v^n)}{\Delta t}$$

$$\frac{u^{n+1}-u^n}{\Delta t}$$

$${\epsilon }^2\frac{v^{n+1}-v^n}{\Delta t}$$

$$\frac{u^{n+1}-u^n}{\Delta t}$$

$${\epsilon }^2\frac{v^{n+1}-v^n}{\Delta t}$$

$$v^{n+1}=\frac{{\epsilon }^2}{{\epsilon }^2+\Delta t}{v}^n-\frac{\Delta t}{{\epsilon }^2+\Delta t}\left(p(u^n{)}_x-f(u^n)\right).$$

$$\frac{u^{n+1}-u^n}{\Delta t}+\frac{{\epsilon }^2}{{\epsilon }^2+\Delta t}{v}_x^n$$

$$\frac{u^{n+1}-u^n}{\Delta t}+\frac{{\epsilon }^2}{{\epsilon }^2+\Delta t}{v}_x^n+\frac{\Delta t}{{\epsilon }^2+\Delta t}f(u^n{)}_x$$

$$\frac{v^{n+1}-v^n}{\Delta t}+\frac{1}{{\epsilon }^2+\Delta t}p(u^n{)}_x$$

$$u_t+\frac{{\epsilon }^2}{{\epsilon }^2+\Delta t}{v}_x+\frac{\Delta t}{{\epsilon }^2+\Delta t}f(u{)}_x$$

$$v_t+\frac{1}{{\epsilon }^2+\Delta t}p(u{)}_x$$

$$\frac{u^{n+1}-u^n}{\Delta t}=u_t+\mathcal{O}(\Delta t),\frac{v^{n+1}-v^n}{\Delta t}=v_t+\mathcal{O}(\Delta t)$$

$${\lambda }_{\pm }^1(\Delta t,\epsilon )=\frac{{\xi }_1}{2}\left(c\pm \sqrt{c^2+\frac{4{\epsilon }^2}{(\Delta t{)}^2}}\right),$$

$${\lambda }_{\pm }^1(0,\epsilon )=\pm \frac{1}{\epsilon }.$$

$${\lambda }_{\pm }^1(\Delta t,0)=\frac{1}{2}\left(c\pm |c|\right).$$

$$\frac{u^{n+1}-u^n}{\Delta t}+f(u^n{)}_x$$

Full text

Key words. IMEX Runge-Kutta methods, hyperbolic conservation laws with sources, diffusion equations, hydrodynamic limits, stiff systems, asymptotic-preserving schemes.

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

1. 1 Introduction

2. 2 Previous IMEX Runge-Kutta approaches

3. 2.1 Additive approach

4. 2.2 Partitioned approach

5. 2.3 Hybrid additive-partitioned approach

6. 3 A unified IMEX Runge-Kutta approach

7. 3.1 Description of the method

8. 3.2 Extension to general IMEX Runge-Kutta schemes

9. 3.2.1 Notations

10. 3.2.2 The unified IMEX-RK setting

11. 3.3 Removing the parabolic stiffness

12. 4 Numerical applications

13. 5 Conclusions

14. A Appendix: IMEX-RK schemes

1 Introduction

Hyperbolic systems with relaxation often contain multiple space-time scales which may differ of several orders of magnitude due to the various physical parameters characterizing the model. This is the case, for example, of kinetic equations close to the hydrodynamic limit [4, 14, 12]. In such regimes these systems can be more conveniently described in terms of fluid-dynamic equations, when they are considered on a suitable space-time scale.

As a prototype system, that we will use to illustrate the subsequent theory, we consider the following

[TABLE]

where $p^{\prime}(u)>0$. System (1) is hyperbolic with two distinct real characteristics speeds $\pm\sqrt{p^{\prime}(u)}/{\varepsilon}^{\alpha}$.

Note that the scaling introduced in (1) corresponds to the study of the limiting behavior of the solution for the usual hyperbolic system with a singular perturbation source

[TABLE]

under the rescaling $t={\varepsilon}^{\alpha}\tau$, $\xi=x$, $v(x,t)=V(\xi,\tau)/{\varepsilon}^{\alpha}$ and $f(u)=F(u)/{\varepsilon}^{\alpha}$.

For $\alpha=0$, system reduces to the usual hyperbolic scaling (2). However, if $\alpha>0$, we are looking at larger microscopic times.

In particular, for small values of ${\varepsilon}$, using the Chapman-Enskog expansion, the behavior of the solution to (1) is, at least formally, governed by the following nonlinear parabolic equation

[TABLE]

Therefore, as ${\varepsilon}\to 0$ when $\alpha\in[0,1)$ we obtain the scalar conservation law

[TABLE]

Note that, the main stability condition [13, 28] for system (3) corresponds to

[TABLE]

and it is always satisfied in the limit ${\varepsilon}\rightarrow 0$ when $\alpha>0$, whereas for $\alpha=0$ it requires suitable assumptions on the functions $f(u)$ and $p(u)$.

In classical kinetic theory the space-time scaling just discussed leads the so-called hydrodynamical limits of the Boltzmann equation (see [12], Chapter 11). For $\alpha=0$ this corresponds to the compressible Euler limit, whereas for $\alpha\in(0,1)$ the incompressible Euler limit is obtained.

Something special happens when $\alpha=1$. In this case, in fact, to leading order in $\varepsilon$, we obtain the parabolic equation

[TABLE]

In other words, considering larger times than those typical for Euler dynamics ($\alpha=1$ instead of $\alpha\in[0,1)$), dissipative effects become non-negligible. This behavior characterizes the incompressible Navier-Stokes limit in classical kinetic theory.

The development of numerical methods to solve hyperbolic systems with stiff source terms in the case $\alpha=0$ has been an active area of research in the past three decades [10, 18, 22, 18, 31, 33, 34, 32, 8]. Another series of works is concerned with the construction of robust schemes for $\alpha=1$ when a diffusion limit is obtained [21, 23, 25, 26]. However, very few papers have considered the general multiscale problem of type (1) for the various possible values of $\alpha$ [29, 24].

The common goal of this general class methods, often referred to as asymptotic-preserving (AP) schemes, was to obtain the macroscopic behavior described by the equilibrium system by solving the original relaxation system (1) with coarse grids $\Delta t,~{}\Delta x\gg{\varepsilon}$, where $\Delta t$ and $\Delta x$ are respectively the time step and the mesh size.

Note that, since the characteristic speeds of the hyperbolic part of system (1) are of order $1/{\varepsilon}^{\alpha}$, most of the popular methods [10, 22, 32], for the solution to hyperbolic conservation laws with stiff relaxation present several limitations when considering the whole range of $\alpha\in[0,1]$ and fail to capture the right behavior of the limit equilibrium equation unless the small relaxation rate is numerically resolved, leading to a stability condition of the form $\Delta t\sim\varepsilon^{\alpha}\Delta x$. Clearly, this hyperbolic stiffness becomes very restrictive when $\alpha>0$, and for $\alpha=1$ in the parabolic regime ${\varepsilon}\ll\Delta x$ where for an explicit scheme a parabolic time step restriction of the type $\Delta t\sim\Delta x^{2}$ is expected. A special class of IMEX schemes with explicit flux and implicit relaxation is able to deal with the parabolic relaxation ($\alpha=1$), [8]. Such methods, however, converge to an explicit scheme for the limit parabolic equation thus requiring a penalization technique to remove the final parabolic stiffness.

In the present paper, using a reformulation of the problem, we develop high-order IMEX Runge-Kutta schemes for a system like (1) in the stiff regime which work uniformly with respect to the scaling parameter $\alpha$. In the parabolic regime, $\alpha=1$, our approach gives a scheme which is not only consistent with (6) without resolving the small $\varepsilon$ scale, but is also capable to avoid the parabolic stiffness leading to the CFL condition $\Delta t\sim\Delta x^{2}$. In other limiting regimes, that is to say when $\alpha\in[0,1)$, the scheme maintains all the nice properties of the numerical schemes for hyperbolic conservation laws, such as the ability to capture shocks with high resolution.

Here, although the final schemes we develop will work independently on ${\varepsilon}$, we will mainly concentrate on the study of the stiff regime for system (1) that is to say when ${\varepsilon}\ll 1$.

Finally we emphasize that from a physical point of view the problem we consider here is close in spirit to the description of the macroscopic incompressible Navier-Stokes equations of fluid-dynamics by the detailed kinetic equations [4, 14, 12]. Although, for the sake of simplicity, we develop our theory for one-dimensional 2$\times$2, systems the results extend far beyond these models.

The rest of the paper is organized as follows. In Section 2 we recall some of the most popular IMEX Runge-Kutta approaches and emphasize their limited applicability when $\alpha$ ranges on the whole $[0,1]$ interval. Next, in Section 3 we introduce our new approach with the aim to avoid the stiffness induced by the characteristic speeds of system (1). First we present the simple first order scheme and then, using the IMEX Runge-Kutta formalism, we construct high order methods. In the case $\alpha=1$, these methods give rise to an explicit approximation of the limiting parabolic problem. In Section 4 we modify the previous schemes in order to avoid the limiting parabolic stiffness. In these schemes the diffusion term in the limiting equation is integrated implicitly. Finally in Section 5 several numerical examples are presented showing the robustness of the present approach. Some final considerations are contained in the last section and an appendix is also included.

2 Previous IMEX Runge-Kutta approaches

In this section, to motivate the new approach, we recall briefly other ways to tackle stiff problems through an implicit-explicit partitioning of the differential system.

We discretize time first, and then we discretize space on the time discrete scheme. The motivation for not adopting a method of line approach is that we can choose the space discretization which is more suitable for each term. For simplicity of notation, in the sequel we assume that $\alpha=1$, and consider the hyperbolic-to-parabolic relaxation. Similar conclusions are obtained for $\alpha\in(0,1)$. Considerations on the case $\alpha=0$ are reported at the end of each subsection.

2.1 Additive approach

Let us consider the simple implicit-explicit Euler method applied to (1) based on taking the fluxes explicitly and the stiff source implicitly [32]

[TABLE]

Solving the second equation for $v^{n+1}$ one obtains

[TABLE]

Making use of this relation (replacing $n$ by $n-1$) in the first equation we get

[TABLE]

Therefore, in the limit $\varepsilon\to 0$, we a two levels scheme (in time) for problem (6)

[TABLE]

Although Eq. (9) is a consistent time discretization of the limit convection-diffusion equation (6), the presence of the term $u^{n-1}$ degrades the accuracy of the first order scheme.

Furthermore, since as $\varepsilon\to 0$ the equilibrium state of Eq. (8)

[TABLE]

involves two time levels, in general additional conditions on the explicit and implicit schemes are necessary in order to obtain asymptotic preserving high order methods. We refer to [8] for more details. These drawbacks are also present for any value of $\alpha\in(0,1]$. Of course, since the additive approach has been originally designed to deal with the case $\alpha=0$ there are no problems in the regime.

2.2 Partitioned approach

A different way to apply the implicit-explicit Euler method to (1) is based on taking the first equation explicitly and the second implicitly [5, 6, 30]

[TABLE]

Solving for $v^{n+1}$ the second equation we get

[TABLE]

which substituted into the first equation results in the two level scheme

[TABLE]

Of course, since the method has been developed specifically for the case $\alpha=1$, in the limit $\varepsilon\to 0$ it yields a consistent explicit scheme for problem (6)

[TABLE]

It is easy to verify that this approach gives a consistent explicit scheme also for any value of $\alpha\in[0,1)$ in the hyperbolic limit. However, for non vanishingly small values of $\varepsilon$, the discretization of the fluxes at different time levels poses several difficulties. For example, if we consider for simplicity $p^{\prime}(u)=1$, introducing the diagonal variables $u\pm\varepsilon v$ system (10) reads

[TABLE]

and therefore the space derivatives of diagonal variables in (13) are defined as a combination of two time levels and it is then not clear how to discretize the system in characteristic variables. Furthermore, most numerical methods based on conservative variables evaluate the flux at the same time level. We refer to [6] for a discussion on these aspects and extensions to schemes that avoid the parabolic stiffness in the relaxation limit for $\alpha=1$.

2.3 Hybrid additive-partitioned approach

A method which combines the advantages of the previous approaches in the various regimes has been proposed in [25, 26]. The idea is to take a convex combination of the previous schemes in such a way that we have an additive scheme in hyperbolic regimes and a partitioned one in parabolic ones. This is achieved by taking

[TABLE]

where $0\leq\phi(\varepsilon)\leq 1$ is such that $\phi(\varepsilon)\approx 1$ for $\varepsilon=\mathcal{O}(1)$ and $\phi(\varepsilon)\approx 0$ when $\varepsilon\ll 1$. For example $\phi(\varepsilon)=\min\{\varepsilon^{2},1\}$ or the smoother approximation $\phi(\varepsilon)=\tanh(\varepsilon^{2})$ were considered in [23, 25].

In the limit $\varepsilon\to 0$ we have the asymptotic behavior of the partitioned approach, whereas for larger values of $\varepsilon$ we have the usual additive approach. The method can be naturally extended to the general multiscale case and provides a consistent discretization also for any value of $\alpha\in[0,1)$. We refer to [24, 29] for further details. One of the main advantages of this approach is that it results in a convex approximation at different times of the space derivative $p(u)_{x}$, therefore one can use different space discretizations for the derivative appearing at time $n$ and the one at time $n+1$. Typically, at time $n$ one can choose a standard hyperbolic discretization that works for large values of ${\varepsilon}$ whereas at time $n+1$ classical central difference schemes that works in the parabolic limit suffices to avoid a CFL condition of the type $\Delta t=\mathcal{O}(\varepsilon)$. In this sense the function $\phi(\varepsilon)$ can be also interpreted as an interpolation parameter between different fluxes in the evaluation of the space derivatives in (10). The determination of the optimal expression of $\phi(\varepsilon)$ in terms of stability and accuracy is a delicate aspect which is beyond the scope of the present paper.

3 A unified IMEX Runge-Kutta approach

In this section we present a different approach which overcomes some of the drawbacks of the above mentioned methods.

3.1 Description of the method

Again let us consider initially, for simplicity of notation, the case $\alpha=1$.

We now consider the following discretization for system (1)

[TABLE]

Solving the second equation for $v^{n+1}$ one obtains

[TABLE]

Making use of this relation in the first equation we get

[TABLE]

Note that, at variance with all the previous approaches, the first equation now uses only two time levels and the space derivatives appear all at the same time level $n$. Therefore we can rewrite the scheme in the equivalent fully explicit form

[TABLE]

In particular, for small values of $\Delta t$ the scheme (17) corresponds to the system

[TABLE]

where we only used

[TABLE]

leaving all other terms.

Note that the left part of system (18) is hyperbolic with characteristic speeds 111The subscript $1$ in the next expression indicates $\alpha=1$.

[TABLE]

with $\xi_{1}={\Delta t}/({\varepsilon^{2}+\Delta t})\in(0,1)$, ${c}=f^{\prime}(u)$ and for simplicity we have set $p^{\prime}(u)=1$.

Now, if $\Delta t\to 0$ for a fixed ${\varepsilon}$, system (18) converges to the original system (1) for $\alpha=1$ and by (19), the characteristics speeds converge to the usual ones, i.e.

[TABLE]

On the other hand, for a fixed $\Delta t$, the characteristic speeds $\lambda^{1}_{+}$ and $\lambda^{1}_{-}$ in (19) are respectively decreasing and increasing functions of $\varepsilon$ and as ${\varepsilon}\to 0$ they converge to

[TABLE]

Therefore, if we denote by $\Delta x$ the space discretization parameter, we obtain the expected hyperbolic CFL condition $\Delta t\leq\Delta x/|{c}|$.

Concerning the second order derivative on the right hand side of (18) this induces a stability restriction of type $\Delta t\sim(\Delta x)^{2}/\xi_{1}$. In particular, since the discrete system (17) as ${\varepsilon}\to 0$ relaxes towards

[TABLE]

in such a limit we have the natural stability condition which links the time step to the square of the mesh space $\Delta t\sim(\Delta x)^{2}$.

Similarly, in the case $\alpha\in[0,1)$, for small values of $\Delta t$ the scheme applied to (1) corresponds to the system

[TABLE]

The left-hand side in (22) now has characteristic speeds

[TABLE]

with $\xi_{\alpha}={\Delta t}/({\varepsilon^{1+\alpha}+\Delta t})\in[0,1)$. As before, if we fix ${\varepsilon}$ and send $\Delta t\to 0$ we obtain the usual characteristic speeds

[TABLE]

The limit behavior of the discrete system now is given by

[TABLE]

and, similarly to the analysis for $\alpha=1$, we now observe that the characteristic speed do not diverge as $\varepsilon\to 0$ since we have

[TABLE]

Therefore, we obtain the expected hyperbolic CFL condition $\Delta t\leq\Delta x/|{c}|$, coming from the hyperbolic part of the system. As before, the stability restriction coming from the parabolic term, is $\Delta t\sim\Delta x^{2}/\xi_{\alpha}$.

3.2 Extension to general IMEX Runge-Kutta schemes

Now we extend the analysis just performed for the simple first order implicit-explicit Euler scheme to a general IMEX-RK scheme.

3.2.1 Notations

An IMEX-RK scheme can be represented with a double Butcher tableau

[TABLE]

where $\tilde{A}=(\tilde{a}_{ij})$ is an $s\times s$ lower triangular matrix with zero diagonal entries for an explicit scheme, and, since computational efficiency in most cases is of paramount importance, usually the $s\times s$ matrix $A=(a_{ij})$ for an implicit scheme is restricted to the particular class of diagonally implicit Runge-Kutta (DIRK) methods, i.e., ($a_{ij}=0,$ for $j>i$). In fact, the use of a DIRK scheme is enough to ensure that the explicit part of the scheme term is always evaluated explicitly (see [1], [11], [7]).

The coefficients $\tilde{c}$ and $c$ are given by the usual relation $\tilde{c}_{i}=\sum_{j=1}^{i-1}\tilde{a_{ij}}$ , $c_{i}=\sum_{j=1}^{i}a_{ij}$ and the vectors $\tilde{b}=(\tilde{b}_{i})_{i=1\cdots s}$ and $b=(b_{i})_{i=1\cdots s}$ provide the quadrature weights to combine internal stages of the RK method.

The order conditions can be derived as usual matching the Taylor expansion of the exact and numerical solution, we refer to [1, 11, 32] for more details. Let us mention that from a practical viewpoint, coupling conditions becomes rather severe if one is interested in very high order schemes (say higher then third).

Here we recall the order conditions for IMEX-RK schemes up to order $p=2$ and $p=3$, under the assumption that $\tilde{c}=c$.

[TABLE]

where we denote by $e^{T}=(1,\cdots,1)\in\mathbb{R}^{s}$, and from the previous relaxation for $\tilde{c}$ and $c$, we have $\tilde{A}e=\tilde{c}$ and $Ae=c$.

It is useful to characterize the different IMEX schemes we consider in this paper according to the structure of the DIRK method. Following [2] we have

Definition 1

We call an IMEX-RK method of type I or type A (see **[32]**) if the matrix $A\in\mathbb{R}^{s\times s}$ is invertible, or equivalently $a_{ii}\neq 0$, $i=1,\ldots,s$. 2. 2.

We call an IMEX-RK method of type II or type CK (see **[11]**) if the matrix $A$ can be written as

[TABLE]

with $a=(a_{21},\ldots,a_{s1})^{T}\in\mathbb{R}^{(s-1)}$ and the submatrix $\hat{A}\in\mathbb{R}^{(s-1)\ \times\ (s-1)}$ is invertible, or equivalently $a_{ii}\neq 0$, $i=2,\ldots,s$. In the special case $a=0$, $w_{1}=0$ the scheme is said to be of type ARS (see **[1]**) and the DIRK method is reducible to a method using $s-1$ stages.

The following definition will be also useful to characterize the properties of the method [6, 8].

Definition 2

We call an IMEX-RK method implicitly stiffly accurate (ISA) if the corresponding DIRK method is stiffly accurate, namely

[TABLE]

If in addition the explicit methods satisfies

[TABLE]

the IMEX-RK method is said to be globally stiffly accurate (GSA).

The definitions of ISA follows naturally from the fact that $s$-stage implicit Runge-Kutta methods for which $a_{si}=b_{i}$ for $i=1,\cdots,\nu$ are called stiffly accurate (see [19] for details). A $\nu$-stage explicit Runge-Kutta method for which $\tilde{a}_{si}=\tilde{b}_{i}$ for $i=1,\cdots s-1$ is called FSAL (First Same As Last, see [20] for details). Note that FSAL methods have the advantage that they require only $s-1$ function evaluations per time step, because the last stage of step $n$ coincides with the first stage of step $n+1$.

Therefore, an IMEX-RK scheme is globally stiffly accurate if the implicit scheme is stiffly accurate and the explicit scheme is FSAL. We observe that this definition states also that the numerical solution of a GSA IMEX-RK scheme coincides exactly with the last internal stage of the scheme.

3.2.2 The unified IMEX-RK setting

Let us consider system (1) and again for simplicity of notation, we set $\alpha=1$. The unified IMEX-RK approach that generalizes first order scheme (15) corresponds to compute first the internal stages $U$ and $V$

[TABLE]

and then the numerical solution:

[TABLE]

where $f(U)$ and $p(U)$ are the vectors with component $f(U)_{i}=f(U^{i})$ and $p(U)_{i}=p(U^{i})$ respectively. Solving the second equation for $V$ in (41) with $\zeta={\varepsilon}^{2}/\Delta t$, one obtains

[TABLE]

Substituting this relation in the first equation we have the resulting IMEX-RK scheme may be written as

[TABLE]

and

[TABLE]

Then setting $B=(\zeta I+A)^{-1}$ for small values of $\Delta t$, from (45), we get the system

[TABLE]

Now, we rewrite $f(U)_{x}=f^{\prime}(U)U_{x}$ and $p(U)x=p^{\prime}(U)U_{x}$, where $f^{\prime}(U)$ and $p^{\prime}(U)$ are diagonal matricies with elements $f^{\prime}(U)_{ii}=f^{\prime}(U^{i})$ and $p^{\prime}(U)_{ii}=p^{\prime}(U^{i})$ respectively. Furthermore, for simplicity we assume $f^{\prime}(U^{i})=c$ and $p^{\prime}(U^{i})=1$.

From the IMEX-RK stages (44) we have

[TABLE]

and using the fact that $U=u^{n}e+\mathcal{O}(\Delta t)$ and $b^{T}e=1$ (consistency of the IMEX-RK scheme) we can write the hyperbolic part in (44) as

[TABLE]

By computing the eigenvalues of the hyperbolic part we obtain

[TABLE]

Next, we show that the above characteristic speeds are limited. Here we need to assume that the scheme is GSA. In fact,

[TABLE]

where for the first case using the ISA property we get $b^{T}A^{-1}e=e^{T}_{s}e=1$ (here $e^{T}_{s}=(0,\cdots,0,1)\in\mathbb{R}^{s}$), whereas for the second case we note that for the consistency of the scheme $b^{T}e=1$. On the other hand,

[TABLE]

where for the first case using the GSA property we get $b^{T}A^{-1}\tilde{A}e=e^{T}_{s}\tilde{A}e=\tilde{b}^{T}e=1$, whereas for the second case we note that the quantity $b^{T}\tilde{A}e$ is a number and if we assume that the scheme is a second order accurate one, this gives $b^{T}\tilde{A}e=1/2$, by (37).

As a consequence we have

[TABLE]

and therefore the CFL condition, in the limit $\varepsilon\to 0$, becomes

[TABLE]

Now, if $\Delta t\ll 1$ we get from

[TABLE]

the parabolic time step restriction $\Delta t\sim\Delta x^{2}$.

Finally, we prove that the scheme (44)-(45) in the limit case $\varepsilon\to 0$ is a consistent discretization of the limit equation (6), i.e. the scheme is AP.

From (43) we get

[TABLE]

and substituting in the numerical solution $v^{n+1}$ we obtain

[TABLE]

Consistency as $\zeta\to 0$, implies $(b^{T}A^{-1}\tilde{A}-\tilde{b}^{T})=0$, which is satisfied if the scheme is GSA, because in this case $b^{T}=e^{T}_{s}A$, $\tilde{b}^{T}=e^{T}_{s}\tilde{A}$, therefore $1-b^{T}A^{-1}e=0$ and $b^{T}A^{-1}\tilde{A}-\tilde{b}^{T}=e^{T}_{s}\tilde{A}-\tilde{b}^{T}=0$. Then we have

[TABLE]

with

[TABLE]

Similarly from (67) we have for the $U$ internal stages

[TABLE]

and for the numerical solution

[TABLE]

therefore, this leads to

[TABLE]

Assuming that the IMEX Runge-Kutta scheme is GSA, the term $b^{T}A^{-1}\tilde{A}=e^{T}_{s}\tilde{A}=\tilde{b}^{T}$ and, as $\varepsilon\to 0$, the scheme relaxes to the explicit one, i.e.,

[TABLE]

Generalization to the case $\alpha\in[0,1).$

In the case $\alpha\in[0,1)$ analogous computations show that the characteristic speeds of the hyperbolic part are given by

[TABLE]

where here $\zeta=\varepsilon^{1+\alpha}/\Delta t$. We have

[TABLE]

A similar analysis performed in this case shows that

[TABLE]

As ${\varepsilon}\to 0$ we get the explicit Runge-Kutta scheme for system (4)

[TABLE]

All the previous results can be stated by the following

Theorem 1

If the IMEX-RK scheme (41)-(42) applied to (1) satisfies the GSA property then, as $\varepsilon\to 0$, it becomes the explicit RK scheme characterized by the pair $(\tilde{A},\tilde{w})$ applied to the limit convection-diffusion equation (6) for $\alpha=1$ and to the limit scalar conservation law (3) for $\alpha\in[0,1)$.

Some remarks are in order.

Remark 1

• •

The advantage of formulating the unified IMEX-RK approach in the form (44)-(45) is that we can now adopt different space discretizations for the various term appearing in the scheme. Typically, we use classical hyperbolic-type schemes, like WENO, for the space derivatives characterizing the hyperbolic part (which now has finite characteristic speeds) and centered discretization for the second order term characterizing the asymptotic parabolic behavior.

• •

If the IMEX-RK scheme satisfies the GSA property the numerical solution is the same as the last stage, then $V^{s}=v^{n+1}$ and, observing that $\tilde{b}_{s}=0$, from the second equation of (45), we have

[TABLE]

Solving for $v^{n+1}$, after some algebra, the equation can be written as

[TABLE]

Assuming $p^{\prime}(u)=1$ and using the fact that $U=u^{n}e+{\cal O}(\Delta t)$ we can now write the hyperbolic part in (44) as

[TABLE]

Therefore, the characteristic speeds read

[TABLE]

and now when ${\varepsilon}\to 0$, as in the first order case discussed in Section 3, we get

[TABLE]

3.3 Removing the parabolic stiffness

Although the final schemes developed in the previous section will work independently on $\varepsilon$ and $\alpha$, for small values of $\varepsilon$ they relax to an explicit RK scheme originating a time step restriction of the type $\Delta t\approx\Delta x^{2}/{\varepsilon}^{1-\alpha}$, for $\alpha\in(0,1]$. Therefore, for small $\varepsilon$, only the case $\alpha=1$ poses stability restriction. For this reason, we shall consider $\alpha=1$ in this subsection.

The natural idea here is to treat the term $p(u)_{x}$ in (1) implicitly and to observe that in the limit case, i.e. $\varepsilon\to 0$, the IMEX-RK scheme relaxes to an IMEX-RK scheme for the limit convection-diffusion equation (6) where the diffusion term is now evaluated implicitly. Compared to similar schemes presented in [6, 8, 9] this new approach has the advantage that is not based on a penalization technique and therefore avoids the difficult problem of the optimal determination of the penalization parameter (see [6]).

The unified IMEX-RK scheme for system (1), now reads

[TABLE]

and

[TABLE]

Now, solving the second equation in (64) for $V$, one obtains

[TABLE]

where $\zeta={\varepsilon}^{2}/\Delta t$. Using this relation in the first equation of (64) we get for the internal stages

[TABLE]

and similarly for the numerical solution

[TABLE]

Then setting $B=(\zeta I+A)^{-1}$ for small values of $\Delta t$, from (68), we get the system

[TABLE]

which is similar to Eq. (46), except that now $A$ and $b^{T}$ appear in front of $p(U)$ terms in place of $\tilde{A}$ and $\tilde{b}$ respectively. Therefore, we have the same characteristic speeds as in the unified IMEX-RK approach presented in the previous section and the same conclusions on the hyperbolic CFL condition holds true.

Concerning the AP property, in the limit $\zeta\to 0$ one has from (66)

[TABLE]

and by the GSA property, i.e. $b^{T}A^{-1}\tilde{A}=e^{T}_{s}\tilde{A}=\tilde{b}^{T}$, we get from (68)

[TABLE]

Thus, scheme (67)-(68), becomes an IMEX-RK scheme for the convection-diffusion equation (6)

[TABLE]

Clearly, scheme (70) is a consistent approximation of the limit equation (6) where now the diffusion term is evaluated implicitly, therefore the CFL condition of such scheme is uniquely determined by the hyperbolic restriction $\Delta t\sim\Delta x$.

A similar analysis in the case $\alpha\in[0,1)$ for scheme (64)-(65) under the GSA assumption produces the explicit Runge-Kutta scheme

[TABLE]

Therefore we can summarize the results in the following

Theorem 2

If the IMEX-RK scheme (64)-(65), applied to (1) for $\alpha=1$, satisfies the GSA property then, as $\varepsilon\to 0$, it becomes the IMEX-RK method characterized by the pairs $(\tilde{A},\tilde{b})$ and $(A,b)$ for the limit convection-diffusion equation (6). Otherwise for $\alpha\in[0,1)$ the IMEX-RK scheme as $\varepsilon\to 0$ yields the explict RK method characterized by the pair $(\tilde{A},\tilde{b})$ for the limit scalar conservation law (4).

4 Numerical applications

In this section we present numerical results that confirm the validity of the new approach presented in section 3.3. In all tests we used the last approach described in section 3.3, so that we remove the parabolic restriction in the limit of small $\varepsilon$ and $\alpha=1$. All the numerical examples presented here refer to the IMEX-RK schemes reported in the Appendix. We shall use the notation NAME$(\nu,\sigma,p)$, where the triplet $(\nu,\sigma,p)$ where $\nu,\sigma$ and $p$ represent respectively the number of explicit function evaluations, the number of implicit function evaluations and the order of accuracy.

In order to avoid spurious numerical oscillations arising near discontinuities of the solutions, we use interpolating non-oscillatory algorithms, like WENO method, [36]. In these numerical test, as emphasized in Remark 1, we use classical hyperbolic-type schemes, like finite difference discretization with WENO reconstruction, for the space derivatives characterized the hyperbolic part while for the second order term $p(u)_{xx}$ we used the standard $2$-th and $4$-th order finite difference technique. Note that when we consider third-order IMEX R-K scheme we use $4$-th order finite difference to discretize the term $p(u)_{xx}$ except at the nearby boundary points where a $3$-rd order formula was implemented, and this guarantees to achieve a global third-order.

In all our numerical results we take $\Delta t=\lambda_{CFL}\Delta x$ in all regimes. The precise choice of $\lambda_{CFL}$ is reported in the figure captions.

As a mathematical model for our numerical experiments we consider the Ruijgrook-Wu model of the discrete kinetic theory of rarefied gases. The model describes a two-speed gas in one space dimension and corresponds to the system [17, 25, 35]

[TABLE]

where $f^{+}$ and $f^{-}$ denote the particle density distribution at time $t$, position $x$ and with velocity $+1$ and $-1$ respectively. Here ${\rm Kn}$ is the Knudsen number, $M$ is the Mach number of the system and $a$,$b$ and $c$ are positive constants which characterize the microscopic interactions. The local (Maxwellian) equilibrium is defined by

[TABLE]

The macroscopic variables for the model are the density $u$ and momentum $v$ defined by

[TABLE]

The nondimensional multiscale problem is obtained taking $M={\varepsilon}^{\alpha}$ and ${\rm Kn}={\varepsilon}$, the Reynolds number of the system is then defined as usual according to $Re=M/Kn=1/{\varepsilon}^{1-\alpha}$. The model, as we will see, has the nice feature to provide nontrivial limit behaviors for several values of $\alpha$ including the corresponding compressible Euler ($\alpha=0$) limit and the incompressible Euler ($\alpha\in(0,1)$) and Navier-Stokes ($\alpha=1$) limits.

Test 1. Diffusive scaling in the linear case

For this numerical test we consider the case $\alpha=1$ with $c=0$, $a=1+A\varepsilon$ and $b=1-A\varepsilon$ in the r.h.s. of (72).

Adding and subtracting the two equations in (72) one obtains the following macroscopic equations for $u$ and $v$

[TABLE]

In the limit $\varepsilon\to 0$ the second equation relaxes to the local equilibrium

[TABLE]

and substituting in the fist equation this gives the limiting advection-diffusion equation

[TABLE]

Next we fix $A=1$ and observe that the limiting advection-diffusion equation admits the exact solution

[TABLE]

on the domain $[-\pi,\pi]$ with periodic boundary conditions. We start to consider as initial conditions (76) at $t=0$, and choose $\varepsilon=10^{-6}$ with final time $T=0.1$. The numerical results are compared with (76) at $t=T$. Relative errors and convergence orders are reported in Tables 1-3. In order to check the temporal order of convergence of these schemes, $\Delta x$ decreases with the time step $\Delta t$ accordingly to the CFL condition $\Delta t=0.5\,\Delta x$. In the Tables we show the order of convergence as

[TABLE]

with $E_{\Delta t}$ the relative error computed with time step $\Delta t$. We consider the error obtained with $N=40$ time steps up to $N=640$ time steps in the interval $[0,T]$. We observe that the classical order of the methods is maintained in the limit case for the density $u$. We omit the convergence results of ARS(2,2,2) scheme (112) since they are essentially the same as the ones obtained with CK(2,2,2) scheme reported in (121).

Typically, classical schemes, as ARS(2,2,2), CK(2,2,2) and BPR(3,4,3), which satisfy only the GSA property, have few numbers of internal stages but maintain, in the limit $\varepsilon\to 0$, the order in time only for the $u$-component, while for the $v$-component they do not guarantee even the consistency [8]. As an example, we report the convergence table for the variable $v$ in Table 2 for the CK(2,2,2) scheme (121) and we observe that we do not obtain the correct classical order of accuracy for $v$ in the limit case $\varepsilon\to 0$.

The reason for the lack of consistency in the algebraic component $v$, is that this classical schemes does not satisfy the additional order conditions (96) (see Appendix) that guarantees the consistency of the scheme and the order up to $2$ for the variable $v$. As a comparison we construct a new scheme, BPR(4,4,2), that satisfies the additional order conditions (96), and in Table 3 we observe the correct convergence rate for both components, $u$ and $v$. Furthermore, in Fig. 1, we also compare the numerical solutions (star points) for $u$ and $v$ obtained with the second order ARS(2,2,2) scheme (112) and BPR(4,4,2) scheme. In the same figure, the exact solution is plotted by a continuous line. Note again that the new scheme BPR(4,4,2) shows the correct behaviour for the $v$ variable.

Of course there are some advantage and disadvantages to consider additional order conditions (96). The advantage is that they guarantee the correct order in the limit case ($\varepsilon\to 0$) for both variables. On the other hand, to construct schemes that satisfy (96), requires a larger number of internal stages, due to these extra order conditions.

Next, we consider a Riemann problem with initial data

[TABLE]

and inflow and outflow boundary conditions. The exact solution for the limit advection-diffusion equation is

[TABLE]

where erf$(x)$ denotes the error function. We take $\Delta x=0.2$ and final time $T=3.0$. In Figure 2 we compare the numerical solutions computed by schemes SP(1,1,1), BPR(2,4,4) and BPR(3,3,5) for the mass density in the intermediate regime ($\varepsilon=0.5$) and in diffusive one ($\varepsilon=10^{-6}$) with a reference solution obtained with $\Delta x=0.001$. As we can see, the numerical results for the different schemes describe the exact motion of the shock and, in the small relaxation limit, are in excellent agreement with the analytical ones using a hyperbolic time step $\Delta t=0.5\Delta x$.

Test 2. Multiscale limit in the nonlinear case

Here we consider the nonlinear Ruijgrok-Wu model Eq.(72) for $\alpha\in[0,1]$ and interaction parameters $c=2{\varepsilon}$ and $a=b=1$.

Adding and subtracting the two equations in (72) one obtains the following macroscopic equations for $u$ and $v$

[TABLE]

For small values of ${\varepsilon}$ the model behaviour can be derived by the Chapman-Enskog expansion and is characterized by the viscous Burgers equation

[TABLE]

We consider two different initial conditions. The first one is given by two local Maxwellian characterized by

[TABLE]

with $v=[(1+u^{2}\varepsilon^{2})^{1/2}-1]/\varepsilon^{2}$.

We show in Figure 3 the numerical solution for $u$ in the case $\alpha=1$ using BPR(4,4,2) and BPR(3,4,3) schemes in the rarefied ($\varepsilon=0.4$) and parabolic ($\varepsilon=10^{-6}$) regimes with $\Delta x=0.2$ at the final time $T=0.2$. The solution of both schemes is in very good agreement with the reference solution computed with $\Delta x=0.04$.

The last test case that we consider is the propagation of an initial square wave. The initial profile is specified as

[TABLE]

with reflecting boundary conditions, i.e. $v=0$, $u_{x}=0$ on the boundary.

We integrate the equations over $[-0.5,0.5]$ with $200$ spatial cells. In Figure 4 we plot the behavior of the system in the rarefied regime for $\varepsilon=0.7$ with $\alpha=0$ at time $T=0.2$ and in the parabolic regime ($\varepsilon=10^{-8}$) with $\alpha=0.5,0.75$. We use BPR(3,4,3) scheme (147). The numerical solutions for the mass density $u$ and momentum $v$ are computed in the rarefied regime and in the diffusive regime with $\Delta t=0.5\Delta x$ and are depicted with a reference solution obtained using fine grids with $\Delta x=0.001$. We observe that in the rarefied regime no oscillation appears and the scheme describes well the behavior of the system. In the parabolic regime the choice of the parameter $\alpha$ gives the corresponding Reynolds numbers of the problem, i.e., $Re=10000$ and $Re=100$. The output of the solution in the parabolic regime is given at $T=0.5$ with $\Delta t=0.004$ (i.e. $\Delta t=0.8\Delta x$) and a right moving shock is formed. Although $\varepsilon$ is underresolved the scheme proposed captures well the correct behavior of the equilibrium equation independently of $\alpha$.

Test 3. Multiscale space varying limit in the nonlinear case

Finally, we consider the multi scale relaxation system of the previous section, in the case where the parameter $\alpha=\alpha(x)\in[0,1]$ depends on the space variable. Therefore, the limit behavior of the system may depend on the particular region of the computational domain.

Now we apply our schemes to system (78) considering two different cases of the varying $\alpha$ number. The domain is chosen to be $x\in[-0.5,0.5]$. In the first case $\alpha$ increases smoothly from a small value $\alpha_{0}$ to $\mathcal{O}(1)$, by the formula,

[TABLE]

with $\alpha_{0}=10^{-6}$. In the second case we consider the function $\alpha$ which contains a discontinuity

[TABLE]

The two different scenarios for $\alpha$ are depicted in the Figure 5.

The initial profile is specified by (81) with reflecting boundary conditions and we fix $\varepsilon=10^{-8}$ and $\Delta t=0.5\Delta x$. Here we use BPR(3,3,5) scheme. The reference solutions are computed using a fine grid with $\Delta x=0.001$. The results are depicted up to time $T=0.05$ (left) and $T=0.18$ (right). In Figure 6 we observe that the scheme is able to capture the correct behavior of the reference solution even in this test case.

5 Conclusions

In this work we have developed a new IMEX-RK approach that is capable correctly to capture the asymptotic behavior of hyperbolic balance laws with relaxation under different kinds of scaling. Previous IMEX-RK schemes were designed to deal specifically with one kind of scaling, either hyperbolic [7, 10, 22, 32] or parabolic [6, 21, 25, 29], whereas the present schemes are robust enough to be able to deal with both scalings. Related approaches were presented in [8, 24, 27, 30].

From a physical viewpoint, these scaling limits corresponds to the classical fluid-dynamic scalings in the kinetic theory of rarefied gases that lead from the Boltzmann equation to its compressible and incompressible limits. Several numerical test for a simple kinetic model which possesses different asymptotic limits have confirmed the validity of the present approach. In the near future, we hope to extend this class of IMEX-RK schemes to the full Boltzmann equation by adopting the penalization techniques developed in [15, 16].

Acknowledgments

The authors are grateful to Prof. M. Lemou and Prof. G. Dimarco for helpful observations. The research that led to the present paper was partially supported by the research grant Numerical methods for uncertainty quantification in hyperbolic and kinetic equations of the group GNCS of INdAM, by ITN-ETN Marie-Curie Horizon 2020 program ModCompShock, Modeling and computation of shocks and interfaces, Project ID: 642768 and by the INDAM-GNCS 2017 research grant Numerical methods for hyperbolic and kinetic equation and applications.

Appendix A Appendix: IMEX-RK schemes

In order to achieve higher then first order accuracy in time for (GSA) IMEX-RK schemes, we need to observe first that there are no second order GSA IMEX-RK methods of type I with $s=3$, i.e. with three stages [6]. Thus, to construct a second order GSA IMEX Runge-Kutta schemes, we consider the following proposition [6]

Proposition 1

The only type of second order GSA IMEX-RK scheme with three levels $s=3$ that satisfies the classical second order conditions (37) is the type II where $c_{i}=\tilde{c}_{i}$ for all $i=2,...,s-1$, with $c_{1}=0$ and $c_{s}=1$.

As usual, order conditions are obtained by matching the Taylor expansion of the exact solution and the numerical one, up to terms of the prescribed order. In the relaxed case, i.e. $\varepsilon\to 0$, system (1) for $\alpha=1$ reduces to

[TABLE]

The unified IMEX-RK approach described in Section 3.3 provides a scheme that converges to an explicit-implicit scheme for the limit convection-diffusion equation in (84). Performing the same analysis as proposed in [8] for system (1) we obtain for the $u$-component the classical order conditions, while some additional order conditions for the $v$-component are required in order to have consistency and maintain the classical order in the limit case. We recall that the GSA assumption of the method guarantees that IMEX-RK scheme relaxes at the same IMEX-RK one when $\varepsilon\to 0$, but in order to maintain the order of accuracy of the scheme in the limit we must impose some additional order conditions.

Below we list these new additional order conditions that the $v$-component must satisfy up to second order

[TABLE]

Note that in order to reduce the number of the additional order conditions we require that $\tilde{c}=c$, this assumption simplifies a lot the number of coupling conditions

[TABLE]

Finally, we present the different IMEX-RK schemes, up to order three, used in our numerical tests. Below, these schemes are represented as usual by the double Butcher tableau. On the left we have the explicit part and on the right the implicit part of the IMEX schemes. All the schemes satisfy the GSA property, but only the new scheme BPR$(4,4,2)$ satisfies the additional order conditions (96).

• •

First order ARS$(1,1,1)$ scheme

[TABLE]

• •

Second order ARS$(2,2,2)$ scheme [1]

[TABLE]

where $\gamma=1-1/\sqrt{2}$ and $\delta=1-1/2\gamma$,

• •

Second order CK$(2,2,2)$ scheme [11]

[TABLE]

• •

Second order BPR$(4,4,2)$ scheme

[TABLE]

• •

Third order BPR$(3,4,3)$ scheme [6]

[TABLE]