An asymptotic preserving scheme for front propagation in a kinetic reaction-transport equation

TL;DR

This paper introduces an asymptotic preserving numerical scheme for a nonlinear kinetic reaction-transport equation, effectively handling sharp interface propagation and stiff regimes by leveraging a micro-macro decomposition aligned with the Hopf-Cole transform.

Contribution

The proposed scheme is specifically designed for nonlinear kinetic equations, maintaining stability and accuracy in the singular limit with a maximum principle and viscosity solution approach.

Findings

Successfully captures sharp interfaces in simulations

Maintains stability across kinetic to macroscopic transition

Efficiently handles stiff regimes with reduced computational cost

Abstract

In this work, we propose an asymptotic preserving scheme for a non-linear kinetic reaction-transport equation, in the regime of sharp interface. With a non-linear reaction term of KPP-type, a phenomenon of front propagation has been proved in [9]. This behavior can be highlighted by considering a suitable hyperbolic limit of the kinetic equation, using a Hopf-Cole transform. It has been proved in [6, 8, 11] that the logarithm of the distribution function then converges to the viscosity solution of a constrained Hamilton-Jacobi equation. The hyperbolic scaling and the Hopf-Cole transform make the kinetic equation stiff. Thus, the numerical resolution of the problem is challenging, since the standard numerical methods usually lead to high computational costs in these regimes. The Asymptotic Preserving (AP) schemes have been typically introduced to deal with this difficulty, since they are designed to be stable along the transition to the macroscopic regime. The scheme we propose is adapted to the non-linearity of the problem, enjoys a discrete maximum principle and solves the limit equation in the sense of viscosity. It is based on a dedicated micro-macro decomposition, attached to the Hopf-Cole transform. As it is well adapted to the singular limit, our scheme is able to cope with singular behaviors in space (sharp interface), and possibly in velocity (concentration in the velocity distribution). Various numerical tests are proposed, to illustrate the properties and the efficiency of our scheme.

Figures

Equations

$${\partial }_{t}f(t,x,v)+v\cdot {\nabla }_{x}f(t,x,v)=\rho (t,x)M(v)-f(t,x,v)+r\rho (t,x)(M(v)-f(t,x,v)),$$

$$\rho (t,x)=⟨f⟩:={\int }_{v\in V}f(t,x,v)\mathrm{d}v,t\ge 0,x\in {\mathbb{R}}^d.$$

$$⟨M⟩=1,⟨vM⟩=0,$$

$$\underset{v\in V}{\inf }M(v)>0.$$

$${\partial }_t{f}^{\epsilon }(t,x,v)+v\cdot {\nabla }_x{f}^{\epsilon }(t,x,v)=\frac{1}{\epsilon }\left({\rho }^{\epsilon }(t,x)M(v)-f^{\epsilon }(t,x,v)+r{\rho }^{\epsilon }(t,x,v)(M(v)-f^{\epsilon }(t,x,v))\right).$$

$$f^{\epsilon }=M{\mathrm{e}}^{-{\psi }^{\epsilon }/\epsilon }.$$

$$0\le f^{\epsilon }(0,\cdot ,\cdot )\le M,$$

$$\forall t\ge 0,0\le {\psi }^{\epsilon }(t,\cdot ,\cdot )\le \Vert {\psi }_{\mathrm{in}}{\Vert }_{\infty }.$$

$${\rho }^{0}M-f^0+r{\rho }^0(M-f^0)=0,$$

$$r{\rho }^0(1-{\rho }^0)=0.$$

$$\psi (0,x,v)={\psi }_{\mathrm{in}}(x),i.e{f}^{\epsilon }(0,x,v)={\mathrm{e}}^{-{\psi }_{\mathrm{in}}(x)/\epsilon }M(v),\text{and}{\rho }^{\epsilon }(0,x)={\mathrm{e}}^{-{\psi }_{\mathrm{in}}(x)/\epsilon },$$

$$\left\{\begin{array}[]{l l}\partial_{t}\psi^{0}+H\left(\partial_{x}\psi^{0}\right)=0,&\text{\;on\;}\mathbb{R}^{*}_{+}\times\mathbb{R},\\ \psi^{0}(0,\cdot)=\psi_{\mathrm{in}}(\cdot).&\end{array}\right.$$

$$\left\{\begin{array}[]{l l}\min\left\{\psi^{0},\partial_{t}\psi^{0}+H\left(\nabla_{x}\psi^{0}\right)+r\right\}=0,&\text{\;on\;}\mathbb{R}^{*}_{+}\times\mathbb{R},\\ \psi^{0}(0,\cdot)=\psi_{\mathrm{in}}(\cdot).&\end{array}\right.$$

$$⟨\frac{M}{1+r+H(p)-vp}⟩=\frac{1}{1+r}.$$

$$\begin{array}[]{c c c}P_{\varepsilon}&\xrightarrow[]{\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \varepsilon\to 0\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ }&P_{0}\vspace{4pt}\\ \rotatebox{90.0}{$\xrightarrow[]{\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \Delta t\to 0\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ }$}&&\rotatebox{90.0}{$\xrightarrow[]{\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \Delta t\to 0\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ }$}\\ S_{\varepsilon}^{\Delta t}&\xrightarrow[]{\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \varepsilon\to 0\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ }&S_{0}^{\Delta t}\end{array}$$

$${\partial }_t{\psi }^{\epsilon }+v{\partial }_x{\psi }^{\epsilon }=1+r⟨M{\mathrm{e}}^{-{\psi }^{\epsilon }/\epsilon }⟩-(1+r)⟨M{\mathrm{e}}^{-{\psi }^{\epsilon }/\epsilon }⟩{\mathrm{e}}^{{\psi }^{\epsilon }/\epsilon },$$

$$f^{\epsilon }={\rho }^{\epsilon }M+g^{\epsilon },$$

$${\psi }^{\epsilon }(t,x,v)={\phi }^{\epsilon }(t,x)+{\eta }^{\epsilon }(t,x,v),(t,x,v)\in {\mathbb{R}}_{+}^{\ast }\times \mathbb{R}\times V,$$

$${\rho }^{\epsilon }(t,x)={\mathrm{e}}^{-{\phi }^{\epsilon }(t,x)/\epsilon },\text{that is}{\phi }^{\epsilon }=-\epsilon \ln ({\rho }^{\epsilon }),$$

$$\frac{f^{\epsilon }}{{\rho }^{\epsilon }M}={\mathrm{e}}^{-{\eta }^{\epsilon }/\epsilon },$$

$$⟨M{\mathrm{e}}^{-{\eta }^{\epsilon }/\epsilon }⟩=1.$$

$$1-{\partial }_t({\phi }^{\epsilon }+{\eta }^{\epsilon })-v{\partial }_x({\phi }^{\epsilon }+{\eta }^{\epsilon })+r{\mathrm{e}}^{-{\phi }^{\epsilon }/\epsilon }=(1+r){\mathrm{e}}^{{\eta }^{\epsilon }/\epsilon }.$$

$$\underset{\epsilon \to 0}{\lim }{\phi }^{\epsilon }={\psi }^0,\underset{\epsilon \to 0}{\lim }{\eta }^{\epsilon }=0,$$

$$\frac{M}{1-{\partial }_t({\phi }^{\epsilon }+{\eta }^{\epsilon })-v{\partial }_x\left({\phi }^{\epsilon }+{\eta }^{\epsilon }\right)+r{\mathrm{e}}^{-{\phi }^{\epsilon }/\epsilon }}=\frac{M{\mathrm{e}}^{-{\eta }^{\epsilon }/\epsilon }}{1+r},$$

$$⟨\frac{M}{1-{\partial }_t({\phi }^{\epsilon }+{\eta }^{\epsilon })-v{\partial }_x\left({\phi }^{\epsilon }+{\eta }^{\epsilon }\right)+r{\mathrm{e}}^{-{\phi }^{\epsilon }/\epsilon }}⟩=\frac{1}{1+r}.$$

$${\partial }_t{\phi }^{\epsilon }+H^{\epsilon }+r=0,$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}\varphi^{\varepsilon}+H^{\varepsilon}+r=0\vspace{4pt}\\ \displaystyle 1-\partial_{t}(\varphi^{\varepsilon}+\eta^{\varepsilon})-v\partial_{x}(\varphi^{\varepsilon}+\eta^{\varepsilon})+r\mathrm{e}^{-\varphi^{\varepsilon}/\varepsilon}=(1+r)\mathrm{e}^{\eta^{\varepsilon}/\varepsilon}\vspace{4pt}\\ \displaystyle\left\langle\frac{M}{1-\partial_{t}\eta^{\varepsilon}+H^{\varepsilon}+r-v\partial_{x}\left(\varphi^{\varepsilon}+\eta^{\varepsilon}\right)+r\mathrm{e}^{-\varphi^{\varepsilon}/\varepsilon}}\right\rangle=\frac{1}{1+r}.\end{array}\right.$$

$$x_i=-x_{\text{max}}+\frac{\Delta x}{2}+(i-1)\Delta x,i=1,\dots ,N_x=2N_x^{\prime },$$

$$v_j=-v_{\text{max}}+\frac{\Delta v}{2}+(j-1)\Delta v,j=1,\dots ,N_v=2N_v^{\prime },$$

$${⟨f⟩}_{N_v}=\Delta v\sum _{j=1}^{N_v}f(v_j).$$

Full text

The hyperbolic scaling and the Hopf-Cole transform make the kinetic equation stiff. Thus, the numerical resolution of the problem is challenging, since the standard numerical methods usually lead to high computational costs in these regimes. The Asymptotic Preserving (AP) schemes have been typically introduced to deal with this difficulty, since they are designed to be stable along the transition to the macroscopic regime. The scheme we propose is adapted to the non-linearity of the problem, enjoys a discrete maximum principle and solves the limit equation in the sense of viscosity. It is based on a dedicated micro-macro decomposition, attached to the Hopf-Cole transform. As it is well adapted to the singular limit, our scheme is able to cope with singular behaviors in space (sharp interface), and possibly in velocity (concentration in the velocity distribution). Various numerical tests are proposed, to illustrate the properties and the efficiency of our scheme.

1 Introduction

We are interested in designing a numerical scheme for a non-linear kinetic equation in the asymptotic regime. The model we consider is a non-linear transport-reaction equation

[TABLE]

with $r\geq 0$, supplemented with an initial data $f(0,x,v)=\rho_{\mathrm{in}}(x)M(v)$. Such models have been introduced in [34, 20, 15]. The asymptotic regime of (1) have been studied in [8, 6, 11], both in the linear case $r=0$, and in the non-linear case $r>0$. In (1), the distribution function $f$, which depends on $t>0$, $x\in\mathbb{R}^{d}$, and $v\in V$, where $V$ is a bounded symmetric set of $\mathbb{R}^{d}$, represents the density of particles at time $t$, at the position $x$, and with velocity $v$. The macroscopic density of particles is defined by,

[TABLE]

Note that the brackets $\left\langle\cdot\right\rangle$ denote velocity average throughout the paper. For $r=0$, equation (1) describes the evolution of the density of particles moving according to a velocity-jump process. Indeed, the motion of a particle is composed of phases of free transport, the run phases, with a velocity $v$, and of tumble phases in which the particle changes velocity instantaneously. The post-tumbling velocity is chosen randomly, according to a given probability density $M$. We assume that $M$ is even, nonnegative, and continuous. Moreover, it satisfies

[TABLE]

and we will suppose that

[TABLE]

Equation (1) is complemented with a reaction term in the case $r>0$. It takes into account creation of new particles at rate $r$, and local quadratic saturation. Initial velocity of new particles is drawn randomly from $M$. Averaging with respect to velocity leads to the classical logistic growth $r\rho(1-\rho)$.

We consider the kinetic equation (1) under an hyperbolic scaling $(t,x,v)\longmapsto(t/\varepsilon,x/\varepsilon,v)$, such that the time and space variables are of same magnitude (see [9]). The kinetic equation (1) then reads

[TABLE]

The propogation of fronts for (1) has been studied in [9]. To study the asymptotic behavior of (4) when $\varepsilon$ goes to [math], an analogy is made in [8, 6] with the sharp front limit of the Fisher-KPP equation. A WKB ansatz is introduced, leading to the so-called approximation of geometric optics (see [17, 19]). It consists in rewriting the distribution function $f^{\varepsilon}$ as

[TABLE]

The equation satisfied by $\psi^{\varepsilon}$ in the limit $\varepsilon\to 0$ is then studied. In the case of the kinetic equation (4), if

[TABLE]

a maximum principle ensures that $\psi^{\varepsilon}$ is well defined and remains nonnegative for all $t\geq 0$, see [6]:

Proposition 1**.**

Let $r\geq 0$ and $\psi_{\mathrm{in}}\in\mathrm{Lip}(\mathbb{R}^{d}\times V)$ and bounded. Let $f^{\varepsilon}=M\mathrm{e}^{-\psi^{\varepsilon}/\varepsilon}$ a solution of (4). Then the phase $\psi^{\varepsilon}$ is uniformly locally Lipschitz, and the following a priori bound holds

[TABLE]

In the case of the Fisher-KPP equation, it has been proved that the function $\psi^{\varepsilon}$ converges to a limit function $\psi^{0}$, which is the viscosity solution of a Hamilton-Jacobi equation, see [17, 2, 3, 35, 13]. Moreover, in the asymptotic regime, the settled population $\rho\sim 1$ is contained in the nullspace of $\psi^{0}$, see [16, 4, 18].

The analysis of propagation phenomena at the mesoscopic scale is motivated by concentration waves of chemotactic bacteria, as observed experimentally [33]. Here, the model under investigation does not contain any chemotactic effect, but takes into account cell division. It satisfies the maximum principle, hence it is more amenable for mathematical analysis, following the seminal works by Kolmogorov, Petrovsky, Piskunov [26], and Aronson, Weinberger [1]. The first analytical works, where travelling waves are construced, are [34] and [15]. Note that the latter develops a micro-macro decomposition to handle the construction of travelling waves near the diffusive regime. We also refer to [20], and references therein, for a more general presentation of reaction transport equations in biology.

The asymptotic behavior of (4) in the limit $\varepsilon\to 0$ has been established rigorously in [8, 6]. Before stating the main theorem, let us highlight that the formation of fronts if $r>0$ can be understood with very formal considerations on (4). Indeed, when $\varepsilon$ goes to [math], supposing that the distribution function $f^{\varepsilon}$ and its density $\rho^{\varepsilon}$ converge respectively to $f^{0}$ and $\rho^{0}$, we have formally at order [math] in $\varepsilon$

[TABLE]

which, once integrated in $v$, gives formally an equation for $\rho^{0}$

[TABLE]

Thus, the limit density $\rho^{0}$ is equal either to [math] or $1$. Going back to (7), in both cases, we obtain that $f^{0}=\rho^{0}M$. Formally, when $\varepsilon$ goes to [math], two areas appear where $f^{0}$ is equal either to [math] or $M$. However, the equation for the dynamic of the propagation of the front is lacking. The WKB ansatz is accurate to catch this phenomenon, since the nullset of $\psi^{0}$ corresponds to the set where $f^{0}=M$, and $f^{0}=0$ where $\psi^{0}$ has positive values. In the one-dimensional case, the equation of the dynamic of the interface, and the asymptotic behavior of $\psi^{\varepsilon}$ is given by the following theorem, proved in [8, 6]:

Theorem 1**.**

Suppose that $V$ is a symmetric and bounded set of $\mathbb{R}$, that $r\geq 0$ and that $M$ is a continuous function of $v$ which satisfies (2)-(3). Let $f^{\varepsilon}$ be the solution of (4), and suppose that the initial data is well prepared

[TABLE]

then $(\psi^{\varepsilon})_{\varepsilon}$ converges locally uniformly towards $\psi^{0}$, where $\psi^{0}$ does not depend on $v$. Moreover, $\psi^{0}$ is the unique viscosity solution of one of the following Hamilton-Jacobi equations:

If $r=0$, then $\psi^{0}$ solves the standard Hamilton-Jacobi problem

[TABLE] 2. 2.

If $r>0$, then the limit equation is the following constrained Hamilton-Jacobi equation

[TABLE]

where, in both cases, the Hamiltonian $H(p)$ is defined implicitly by

[TABLE]

If $d\geq 2$ or if (3) is not satisfied, it may happen that the implicit definition for the Hamiltonian $H(p)$ (10) is inconsistent for large values of $p$. The definition of the Hamiltonian has been extended in [11] to deal with large $p$ in full generality. We refer to Section 6.4 for the details. Numerical results are presented in order to illustrate the singular behaviours which arise in this case (concentration in the velocity variable). We also refer to the recent work in [7] for propagation of front in higher dimension.

The goal of this paper is to construct a numerical scheme able to compute the solution of (4) in all regimes of $\varepsilon$, in the one-dimensional case. Indeed, the fast relaxation towards the equilibrium distribution function $M$ and the outbreak of a sharp interface in the space variable, make the rescaled kinetic equation (4) stiff for small $\varepsilon$. Standard numerical methods for partial differential equations require in general the use of refined grids, when they do not take into account the special structure of the problem. Asymptotic Preserving (AP) schemes have been introduced to avoid the problems arising with standard numerical methods when solving stiff asymptotic problems, see [21, 24, 25]. Indeed, they are constructed exactly with the purpose of being stable along the transition from the mesoscopic regime ($\varepsilon\sim 1$) to the macroscopic regime ($\varepsilon\ll 1$). Their properties are often summarized by the following diagram

[TABLE]

It can be understood as follows: considering an $\varepsilon$-dependent problem $P_{\varepsilon}$ which converges to a limit problem $P_{0}$ when $\varepsilon$ goes to [math] (see Th. 1), the AP scheme $S_{\varepsilon}^{\Delta t}$ must be consistent with $P_{\varepsilon}$ when $\varepsilon$ is fixed. In addition, when the discretization parameters $\Delta t$ are fixed, it has to converge when $\varepsilon$ goes to [math] to a limit scheme $S_{0}^{\Delta t}$, which is required to be consistent with $P_{0}$. An AP scheme can also enjoy the stronger property of being Uniformly Accurate (UA), which means that its accuracy does not depend on $\varepsilon$.

The development of AP schemes for stiff kinetic equations is necessary to ensure their numerical resolution in all regimes of $\varepsilon$. Many AP strategies have been proposed for various asymptotics of linear kinetic equation, see for instance [12, 22, 23, 10, 32, 29, 28, 5].

For the particular asymptotics we are considering, an asymptotic scheme has been proposed in [30], for the linear case $r=0$ in (4), complemented with an efficient scheme for the limit equation (8) in [31]. Asymptotic schemes are based on expansions of the solution of the equation in formal power series of $\varepsilon$. Contrary to AP schemes, they are designed to approach the asymptotic behavior of the equation, but are not relevant for the mesoscopic regimes, where $\varepsilon\sim 1$. Indeed, the formal power series in $\varepsilon$ are truncated, which restrains accuracy of the scheme to the small $\varepsilon$ regimes. Moreover, the scheme proposed in [30] is not adapted to the non-linear case $r>0$ in (4), and it is not clear that the approach they use can be easily adapted to the non-linear case. Conversely, AP schemes are designed to be efficient in all regimes of $\varepsilon$. To do so, it is ensured that no approximation in $\varepsilon$ is done when constructing the scheme. As a consequence, the scheme is consistent with the equation when any $\varepsilon>0$ is fixed. In addition, since they are constructed to be efficient in the asymptotic regime, they catch perfectly the asymptotic behavior of the equation. Note also that, contrary to the asymptotic schemes, the UA property provides in addition the correct behavior of the numerical solutions in the intermediate regime, with no constraint on the parameter $\varepsilon$. As the scheme proposed in this paper enjoys the AP property even in the non-linear case $r>0$ in (4), and since the numerical tests suggest that it also enjoys the UA property, it is accurate in all regimes of $\varepsilon$.

We propose here a scheme for (4) in the one-dimensional case, which enjoys the AP property for the limit equations (8)-(9). It is based on the following reformulation of (4)

[TABLE]

where $\psi^{\varepsilon}$ is the Hopf-Cole transform of $f^{\varepsilon}$, defined in (5). This problem is still stiff when $\varepsilon$ goes to [math], and its non-linear character must be carefully dealt with. Indeed, it is necessary to ensure that the computational cost and the accuracy of the resolution of the non-linear system are constant with respect to $\varepsilon$. Another difficulty we have to consider, is the fact that the solution of (11) converges when $\varepsilon$ goes to [math] to the viscosity solution of one of the Hamilton-Jacobi equations (8)-(9). The discretization of the kinetic equation (4) has to be carefully discussed, since it yields the discretization of the limit equation, and that it is required to catch the viscosity solution of the Hamilton-Jacobi equation. Moreover, it must be robust enough to deal with the possible lack of regularity of the solutions of (4) in the asymptotic regime. Indeed, the solution is expected to be only locally Lipschitz regular, and may present some $\mathcal{C}^{1}$ discontinuities. Eventually, when $r>0$ the limit equation (9) is a constrained problem, thus the AP scheme will be designed to respect this constraint.

The scheme we proposed is based on an adaptation of micro-macro schemes for linear kinetic equations (see [28, 29, 27]). The micro-macro decomposition is modified to be compatible with the nonlinearity of (11). In addition, the resolution of the nonlinear system, that is needed to compute the numerical solution of (11), is carefully discussed to ensure the computational cost of the method does not depend on $\varepsilon$. The scheme enjoys the AP property, and the limit scheme catches the viscosity solution of the limit equation. It also satisfies a maximum principle, which is a discrete equivalent of (6).

The paper is organized as follows: the scheme and a formal derivation of the asymptotic behavior of (11), on which the construction of the scheme is based, are presented in the next section. Then, the numerical resolution of the non-linear system is discussed in Section 3, and the discrete maximum principle is proved in Section 4. The proof of the AP character of the scheme is completed in Section 5. Eventually, various numerical tests are proposed in Section 6, to highlight the properties of the scheme.

Acknowledgements. The author wishes to thank Vincent Calvez for having suggested this problem to her, and for very interesting discussions about it. She would also like to thank Emeric Bouin for his kind help.

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon $2020$ research and innovation programme (ERC starting grant MESOPROBIO no$639638$).

2 The numerical scheme

In this section, we construct an AP scheme for (4) in the one-dimensional case, using its reformulation with the Hopf-Cole transform of $f^{\varepsilon}$ (11). We focus here on the relevant one-dimensional case, as usually done when investigating front propagation. We refer to [7] for the analysis of front propagation in the higher dimensional case. This scheme is based on an adaptation of the micro-macro decomposition proposed in [29, 28, 27] for the construction of AP schemes near the diffusive regime of linear kinetic equations Indeed, the usual micro-macro decomposition consists in decomposing the distribution function $f^{\varepsilon}$ as the sum of its part at equilibrium and of a remainder,

[TABLE]

where $\left\langle f^{\varepsilon}\right\rangle=\rho^{\varepsilon}$, and $\left\langle g^{\varepsilon}\right\rangle=0$. Once injected in a linear kinetic equation, an integration with respect to velocity yields an equation for $\rho^{\varepsilon}$. It is then substracted from the original kinetic equation, to get a coupled system of partial differential equations for $\rho^{\varepsilon}$ and $g^{\varepsilon}$. It is equivalent to the original kinetic equation, and the construction of a an AP scheme for the kinetic equation is more explicit with this version. In the case (11) we are considering, such an approach is not adapted, since the non-linearity of the collision operator prevents the simplifications arising in the linear case. Instead, we propose a new micro-macro decomposition of $f^{\varepsilon}$, in which the density appears multiplied by a remainder. Thanks to the logarithmic transform applied to $f^{\varepsilon}$, this becomes a linear decomposition for $\psi^{\varepsilon}$. In the linear case, the micro-macro decomposition is based on a Chapman-Enskog expansion of the distribution function. Here, it makes also appear the limit equation, and it yields in addition an equation for the corrector. Namely, we decompose $\psi^{\varepsilon}$ defined in (5) as follows

[TABLE]

where $\varphi^{\varepsilon}$ stands for the Hopf-Cole transform of $\rho^{\varepsilon}$

[TABLE]

and $\eta^{\varepsilon}$ denotes what remains to reconstruct $f^{\varepsilon}$

[TABLE]

such that

[TABLE]

Once injected in (11), this decomposition and the conserved quantity (14) yield

[TABLE]

With this reformulation of (11), it is easy to highlight formally the asymptotic behavior of the model, presented in Theorem 1. Moreover, the formal computations we present here are also the basis of the construction of the scheme we propose. Theorem 1 states that when $\varepsilon$ goes to [math], the distribution function $f^{\varepsilon}$ converges to a distribution at equilibrium, since $\psi^{\varepsilon}$ converges to a phase $\psi^{0}$ which does not depend on $v$, and which solves one of the limit equations (8)-(9). The asymptotic behavior of $\varphi^{\varepsilon}$ and $\eta^{\varepsilon}$ is formally,

[TABLE]

since $\varphi^{\varepsilon}$ is the Hopf-Cole transform of $\rho^{\varepsilon}$. Moreover, the equation safisfied by $\psi^{0}$ appears naturally from (15). Indeed, it can be reformulated as follows

[TABLE]

and it becomes, after an integration with respect to velocity

[TABLE]

If we replace $\partial_{t}\varphi^{\varepsilon}$ by $-H^{\varepsilon}-r$, the structure of the limit Hamilton-Jacobi equations (8)-(9) appears. Indeed, we define $H^{\varepsilon}$ by

[TABLE]

and it can be injected in (17). Thanks to (16), if $r=0$ the previous expression formally goes to (8) when $\varepsilon$ goes to [math], and if $\varphi^{0}>0$, we obtain (9) when $\varepsilon$ goes to [math]. Note that (15) ensures that the previous equation (17) is equivalent to the conservation of the quantity (14).

This formal asymptotic analysis paves the way for designing an AP scheme for (4). Indeed, as for the usual micro-macro decomposition, the expressions (15) and (17) provide a coupled system of equations for $\varphi^{\varepsilon}$ and $\eta^{\varepsilon}$, on which the asymptotic equation appears clearly

[TABLE]

To simplify the expression of the scheme, we will rather remplace the last line by (14).

The scheme for (18) is written using a symmetric grid for the space variable

[TABLE]

where $\Delta x=x_{\text{max}}/N_{x}^{\prime}$, and for the velocity variable

[TABLE]

with $\Delta v=v_{\text{max}}/N_{v}^{\prime}$. The integrations in $v$ will be denoted by $\left\langle\cdot\right\rangle_{N_{v}}$ and performed with a simple quadrature rule

[TABLE]

In order to write an upwind scheme for the transport in the $x$ variable, we define

[TABLE]

For the time variable, we will denote by $T$ the final time, $\Delta t=T/N_{t}$ the time step, and

[TABLE]

As in [6, 8], the initial data for $f^{\varepsilon}$ will be considered at equilibrium $f(0,\cdot,\cdot)=M\rho_{\text{in}}^{\varepsilon}=M\mathrm{e}^{-{\varphi_{\text{in}}}/{\varepsilon}}$, such that $\varphi_{\text{in}}$ is a given function and $\eta^{\varepsilon}(0,\cdot,\cdot)=0$. Denoting $\eta^{n}_{i,j}$ (resp. $\varphi^{n}_{i}$) an approximation of $\eta^{\varepsilon}(t_{n},x_{i},v_{j})$ (resp. $\varphi^{\varepsilon}(t_{n},x_{i})$) for all $n\geq 1,i\in[\![1,N_{x}]\!],j\in[\![1,N_{v}]\!]$, we propose the following scheme for $\varphi^{\varepsilon}$ and $\eta^{\varepsilon}$

[TABLE]

with the initialization $\varphi^{0}_{i}=\varphi_{in}(x_{i}),\;\eta^{0}_{i,j}=0$, and where the transport part is computed with an upwind scheme

[TABLE]

Note that the dependence in $\varepsilon$ of $\varphi^{n}_{i}$ and $\eta^{n}_{i,j}$ has been omitted to simplify the notations. To write this scheme, we implicited the stiff terms of the equation, to avoid the problems that may arise in the numerical resolution. From the maximum principle, $\varphi^{\varepsilon}$ is nonnegative in (18). Altough $\mathrm{e}^{-\varphi^{\varepsilon}/\varepsilon}$ is then not stiff, we decided to implicit it in the scheme. Indeed, with this formulation, a discrete maximum principle can be proved for the scheme (see Section 4), and we observed that it is not satisfied in the numerical tests if $\mathrm{e}^{-\varphi^{\varepsilon}/\varepsilon}$ is taken explicit. Eventually, the use of an upwind scheme (24) for the transport part is important, since it provides the monotonicity properties needed to prove that the limit scheme catches the viscosity solution of the asymptotic model (8)-(9).

In practise, to compute $\left(\varphi^{n+1}_{i},\eta^{n+1}_{i,j}\right)_{(i,j)\in[\![1,N_{x}]\!]\times[\![1,N_{v}]\!]}$, one must solve the non-linear system

[TABLE]

in which the unknowns are $\left(H^{n+1}_{i},\eta^{n+1}_{i,j}\right)_{(i,j)\in[\![1,N_{x}]\!]\times[\![1,N_{v}]\!]}$. To ensure the AP property of the scheme, the computational cost of the numerical resolution of this system must not depend on $\varepsilon$. A Newton’s method is used for its resolution, it is discussed in the next section. Eventually, we have the following proposition:

Proposition 2**.**

Consider the scheme (23) defined for all $i\in[\![1,N_{x}]\!],j\in[\![1,N_{v}]\!],n\in[\![0,N_{t}-1]\!]$, and suppose that the accuracy and computational cost of the resolution of the non-linear problem (25) are independent of $\varepsilon$. This scheme has the following properties:

The scheme is of order $1$ for any fixed $\varepsilon>0$. 2. 2.

The scheme enjoys a discrete maximum principle : Let $m>0$ such that $0\leq\varphi^{0}_{i}\leq m$, and suppose that $\eta^{0}_{i,j}=0$ for all $(i,j)\in[\![1,N_{x}]\!]\times[\![1,N_{v}]\!]$. Then, the following bounds hold: $\forall n\in[\![0,N_{t}]\!],\forall(i,j)\in[\![1,N_{x}]\!]\times[\![1,N_{v}]\!],$

[TABLE] 3. 3.

The scheme is AP: when the discretization parameters are fixed, the scheme tends when $\varepsilon$ goes to [math] to one of the following scheme:

[TABLE]

with $\varphi^{n}_{i}\geq 0$ for all $(n,i)\in[\![1,N_{t}]\!]\times[\![1,N_{x}]\!]$, and where in both cases $H^{n+1}_{i}$ is defined by

[TABLE]

These two schemes catch respectively the viscosity solutions of the asymptotic problems (8)-(9), when the discretization parameters tend to zero.

Remark*.*

The numerical tests we propose in Section 6.3 suggest that the scheme enjoys the UA property.

The consistency and the order of the scheme for fixed $\varepsilon$ are clear, since we only used finite differences to write it. The second point of this proposition is proved in Section 4 and the AP character of the scheme is proved in Section 5. The properties claimed in this proposition are illustrated by various numerical tests in Section 6.

3 Resolution of the non-linear system

The scheme (23) proposed in the previous section is both non-linear and implicit. In practice, being given the values $\varphi^{n}_{i},\;\eta^{n}_{i,j}$, the numerical resolution of the non-linear system (25) is needed to compute $\varphi^{n+1}_{i},\;\eta^{n+1}_{i,j}$. Furthermore, since the system (25) depends on $\varepsilon$, and contains some stiff terms, its numerical resolution may become costly when $\varepsilon$ goes to [math]. This would break the AP property, as the computational cost is expected not to increase when $\varepsilon$ decreases. In this section, we explain how the Newton method can be tuned to resolve this issue, and how it can be implemented to make the computational cost independent on the stiffness of the problem.

The system (25) contains $N_{x}(N_{v}+1)$ unknowns, for the same number of equations. We propose to consider it as $N_{x}$ systems of $N_{v}+1$ unknowns, each of them being numerically solved with Newton’s method. Indeed, if $i_{0}\in[\![1,N_{x}]\!]$ is fixed, the system

[TABLE]

where the unknowns are $\left(\left(\eta^{n+1}_{i_{0},j}\right)_{j\in[\![1,N_{v}]\!]},H^{n+1}_{i_{0}}\right)$, is a closed system of $N_{v}+1$ unknowns and $N_{v}+1$ equations, which can be reformulated as follows

[TABLE]

with $F:\mathbb{R}^{N_{v}+1}\longrightarrow\mathbb{R}^{N_{v}+1}$. The Jacobian matrix of $F$ can be computed explicitly to apply Newton’s method. It writes

[TABLE]

where the coefficients $\alpha_{j},\gamma_{j},\delta_{j}$ are defined for $j\in[\![1,N_{v}]\!]$ by

[TABLE]

Note that their dependence in $i_{0}$ is omitted to simplify the notations. Since the coefficients of (30) are stiff when $\varepsilon$ is small, it is more convenient to compute analytically the inverse of this matrix, instead of doing it numerically. Indeed, this can be computed explicitly, and it reads

[TABLE]

where

[TABLE]

Having an explicit formula for $\left(DF\left(\left(G^{n+1}_{i_{0},j}\right)_{j\in[\![1,N_{v}]\!]},H^{n+1}_{i_{0}}\right)\right)^{-1}$ presents two advantages. First of all, it reduces the computational cost by avoiding the numerical inversion of (30), which may be of large dimension. One can also notice that the construction of the matrix (31) can be done with a moderate computational cost, since it is mostly made of an assembly of the vectors $\left(\delta_{i}/\alpha_{i}\right)_{i}$, $\left(1/\alpha_{i}\right)_{i}$, and $\left(\gamma_{i}/\alpha_{i}\right)_{i}$. But the main advantage is that the coefficients of (31), which contain

[TABLE]

and

[TABLE]

do not present any singularity as $\varepsilon$ goes to [math]. Indeed, thanks to the discrete maximum principle that is proved in Section 4, $\mathrm{e}^{(r\Delta t-\varphi^{n}_{i}+H^{n+1}_{i}\Delta t)/\varepsilon}$ is bounded with respect to $\varepsilon$. In fact, the discrete maximum principle states that $\varphi^{n+1}_{i}$ defined by the first line of (23) as

[TABLE]

remains nonnegative. The proof of the AP character of the scheme in Section 5 also states as a preliminary result that $\mathrm{e}^{{\eta^{n+1}_{i_{0},j}}/{\varepsilon}}$ and $\mathrm{e}^{-{\eta^{n+1}_{i_{0},j}}/{\varepsilon}}$ are bounded with respect to $\varepsilon$. Using (31) in the numerical computations will avoid singularities in the numerical computation of the inverse of (30), that may appear for the small values of $\varepsilon$.

As a conclusion, we propose to compute the solution of (29) with the Newton method initialized by

[TABLE]

and with the iterations

[TABLE]

Thanks to the discussion above, the accuracy and the computational cost of the resolution of the system is then independant of $\varepsilon$, which is a necessary condition for the scheme (23) to enjoy the AP property.

4 Discrete maximum principle

The scheme (23) has been designed including an upwind choice to deal with the transport terms. This ensures the consistency of the scheme for a given $\varepsilon$, as well as its stability under a CFL condition

[TABLE]

The stiff terms in $\varepsilon$ have been treated implicitly to ensure the stability of the scheme when $\varepsilon$ decreases. Notice that this implies the AP property as well, as we shall see in the next section. Moreover, the scheme enjoys a discrete maximum principle, which is the discrete version of the maximum principle (6) proved in [6, 8] in the continuous case. In this section, we prove the the following proposition, which contains the second point of Proposition 2.

Proposition 3**.**

Let $m>0$. Assume that for all $i\in[\![1,\;N_{x}]\!],j\in[\![1,N_{v}]\!],0\leq\varphi^{0}_{i}\leq m$ and $\eta^{0}_{i,j}=0$, then for all $n\in[\![0,N_{t}-1]\!],$ and for all $(i,j)\in[\![1,N_{x}]\!]\times[\![1,N_{v}]\!]$, the following properties hold :

$0\leq\varphi^{n+1}_{i}\leq m$, 2. 2.

$0\leq\eta^{n+1}_{i,j}+\varphi^{n+1}_{j}\leq m$, 3. 3.

$0\leq\varphi^{n+1}_{i}+\eta^{n+1}_{i,j}-\Delta t\left(1-\mathrm{e}^{\eta^{n+1}_{i,j}/\varepsilon}\right)+r\Delta t\mathrm{e}^{\eta^{n+1}_{i,j}/\varepsilon}\left(1-\mathrm{e}^{-\left(\varphi^{n+1}_{i}+\eta^{n+1}_{i,j}\right)/\varepsilon}\right)\leq m$.

Proof.

As a preliminary result, let us remark that for a fixed $i\in[\![1,N_{x}]\!]$, there exists at least one index $j_{-}$ such that $\eta^{n+1}_{i,j_{-}}\leq 0$ and one index $j_{+}$ such that $\eta^{n+1}_{i,j_{+}}\geq 0$. Otherwise, the equality

[TABLE]

can not be fulfilled, since $\left\langle M\right\rangle=1$.

The proof of this proposition is done by induction. Let us fix $n\in[\![0,N_{t}-1]\!]$ and suppose that for all $(i,j)\in[\![1,N_{x}]\!]\times[\![1,N_{v}]\!]$, $0\leq\varphi^{n}_{i}+\eta^{n}_{i,j}\leq m$. The third point of the proposition comes by fixing $j\in[\![1,N_{v}]\!]$ such that $v_{j}\geq 0$ (the demonstration is similar for $v_{j}<0$). Indeed, the second line of (23) can be recast as follows

[TABLE]

or, equivalently

[TABLE]

If $0\leq\varphi^{n}_{i}+\eta^{n}_{i,j}\leq m$ for all $(i,j)\in[\![1,N_{x}]\!]\times[\![1,N_{v}]\!]$, and under the CFL condition (32), the right-hand side of (33) is a convex combination of quantities lying in the interval $[0,m]$. This yields the inequality

[TABLE]

Let us suppose in addition that for all $i\in[\![1,N_{x}]\!]$, $\varphi^{n}_{i}\geq 0$. The first point of the proposition is a consequence of the previous inequality. We argue by contradiction. Suppose that $\varphi^{n+1}_{i_{0}}<0$ for some $i_{0}$. In addition, let $j_{0}\in[\![1,N_{v}]\!]$ such that $\eta^{n+1}_{i_{0},j_{0}}\leq 0$. Then the inequality

[TABLE]

cannot be fulfilled, since the right-hand side is a sum of negative terms. Hence, for all $j\in[\![1,N_{v}]\!],\;\eta^{n+1}_{i_{0},j}>0$, but this is contradictory with the preliminary result stated in the beginning of this proof. Hence, for all $i\in[\![1,N_{x}]\!]$, $\varphi^{n+1}_{i}\geq 0$. Let $j_{0}$ such that $\eta^{n+1}_{i,j_{0}}\geq 0$. Considering (34) for $(i,j_{0})$, we obtain

[TABLE]

which gives the upper bound for $\varphi_{i}^{n+1}$, and the first point of the proposition

[TABLE]

Eventually, the second point of the proposition can be proved. Indeed, if $\eta^{n+1}_{i,j}\geq 0$, the inequality $0\leq\varphi^{n+1}_{i}+\eta^{n+1}_{i,j}$ follows immediately, and if $\eta^{n+1}_{i,j}<0$, we have

[TABLE]

which proves the positivity of $\varphi^{n+1}_{i}+\eta^{n+1}_{i,j}$. This implies also

[TABLE]

Considering the function $f:x\mapsto x+\Delta t(\mathrm{e}^{x/\varepsilon}-1)$, we remark that since $x/\varepsilon\leq\mathrm{e}^{x/\varepsilon}-1$,

[TABLE]

because $m-\varphi^{n+1}_{i}\geq 0$. Eventually, we get the upper bound of $\varphi^{n+1}_{i}+\eta^{n+1}_{i,j}$

[TABLE]

∎

5 AP character

In this section, we prove that when $\varepsilon$ goes to [math], the quantity $\eta^{n}_{i,j}$, which depends on $\varepsilon$ is such that $\eta^{n}_{i,j}/\varepsilon$ is bounded uniformly with respect to $\varepsilon$, and that the limit scheme for (23) when $\varepsilon$ goes to [math] is consistent with (9). Let us fix $n\in[\![0,N]\!]$, and $(i,j)\in[\![1,N_{x}]\!]\times[\![1,N_{v}]\!]$. Since

[TABLE]

where the discrete integration with respect to velocity has been defined in (21), and $\left\langle M\right\rangle_{N_{v}}=1$, there exists a constant $c$ such that

[TABLE]

Moreover, as

[TABLE]

with,

[TABLE]

we obtain that

[TABLE]

Hence, both $\mathrm{e}^{-\eta^{n}_{i,j}/\varepsilon}$ and $\mathrm{e}^{\eta^{n}_{i,j}/\varepsilon}$ are bounded independently of $\varepsilon$, thanks to the maximum principle of the previous section. As a consequence, $\eta^{n}_{i,j}/\varepsilon$ is bounded independently of $\varepsilon$. In other words, $\eta^{n}_{i,j}$ vanishes when $\varepsilon$ goes to [math], which is the expected asymptotic behavior of $\eta^{\varepsilon}$ at the continuous level (16). To find the limit scheme for (23), it is more convenient to rewrite the third line of (23) as

[TABLE]

These two expressions are equivalent, thanks to the second line of the scheme (23). Hence, when $\varepsilon\to 0$ with fixed $\Delta t$, (23) becomes a scheme for $\varphi^{n}_{i}$ only that writes

[TABLE]

if $r=0$, and

[TABLE]

if $r>0$. In both cases, the maximum principle of the second point of Prop. 2 ensures that $\varphi^{n}_{i}\geq 0$, and $H^{n+1}_{i}$ is defined implicitly for all $n\in[\![1,N-1]\!]$, and $i\in[\![1,N_{x}]\!]$, by

[TABLE]

where the space derivative $\left[v\partial_{x}\varphi\right]^{n}_{i,j}$ is computed with an upwind scheme as in (24).

To prove that the scheme (23) enjoys the AP property, it remains to show that (26) and (27) catch the viscosity solution of (8) and (9). To do so, we use a general result on numerical schemes for Hamilton-Jacobi equations, stated in [14]. It states that, for a scheme written in the general form

[TABLE]

that can also be written with a differenced form

[TABLE]

the following theorem, adapted from [14], holds

Theorem 2**.**

Let $H:\mathbb{R}\longrightarrow\mathbb{R}$ be continuous, $\varphi_{\mathrm{in}}$ be bounded and Lipschitz continuous on $\mathbb{R}$ with Lipschitz constant $L$, and $\varphi$ be the viscosity solution of

[TABLE]

Under the $CFL$ condition

[TABLE]

suppose that the scheme (38) has a differenced form (39), that the numerical Hamiltonian $h$ is consistent with (40) with $H$, and locally Lipschitz continuous. Suppose also that the scheme $F$ defined in (38) is a nondecreasing function of each argument. Then, there exists a constant $c$ depending only on $\sup{|\varphi_{\text{in}}|}$, $L$, $h$ and $T$ such that

[TABLE]

Here, the discretized Hamiltonian $H_{N_{v}}$ and the numerical Hamiltonian $h$ are implicitly defined. Indeed, if $p\in\mathbb{R}$, $H_{N_{v}}(p)$ is solution of

[TABLE]

with $1+H_{N_{v}}(p)-vp>0$, and if $(p,q)\in\mathbb{R}^{2}$, $h(p,q)$ is solution of

[TABLE]

with $1+r+h(p,q)-v_{+}p-v_{-}q>0$. Therefore, the implicit function theorem ensures that $H_{N_{v}}$ (resp. $h$) is a continuous and locally Lipschitz function of $p$ (resp. $(p,q)$). The consistency of $h$ with $H_{N_{v}}$ is a consequence of the fact that $h(p,p)=H_{N_{v}}(p)$. Note that we only consider the discrete integrations $\left\langle\cdot\right\rangle_{N_{v}}$ for this result. The fact that the solution of (8)-(9) with $H_{N_{v}}$ converges to the solution of (8)-(9) with $H$ is a consequence of stability results for Hamilton-Jacobi equations (see [2]). To prove that the schemes (26)-(27) catch the viscosity solution of (8)-(9) with the integrations in velocity done with $\left\langle\cdot\right\rangle_{N_{v}}$, it remains to show that it is nondecreasing in each argument. Considering (26)-(27) in the general form (38), we compute the partial derivatives of $F$ with the implicit function theorem. They read

[TABLE]

which gives bounds on the partial derivatives of $F$,

[TABLE]

The CFL condition (41) ensures the second one to be nonnegative. The schemes (26)-(27) then fulfill the hypothesis of Theorem 2. Therefore, it converges to the viscosity solution of (8)-(9) with the integrations in velocity done with $\left\langle\cdot\right\rangle_{N_{v}}$, when the discretization parameters $\Delta t$ and $\Delta x$ tend to [math].

6 Numerical tests

In this section, we present some numerical tests for the scheme (23). In all the simulations, we will use the discretizations (19)-(20)-(22), with $x_{\text{max}}=1,\Delta v=1.25\cdot 10^{-2}$, and $v_{\max}=1$. The time and space steps $\Delta t$ and $\Delta x$ will be given in each case. Section 6.4 excepted, the equilibrium $M$ we consider is constant on $[-v_{\max},v_{\max}]$, such that $\left\langle M\right\rangle_{N_{v}}=1$.

To compare the results given by (23) for large values of $\varepsilon$, we will compute numerical solutions of (4) using an explicit scheme for (4). Denoting $f^{\varepsilon}(t_{n},x_{i},v_{j})\sim f^{n}_{i,j}$ and $\rho^{\varepsilon}(t_{n},x_{i})\sim\rho^{n}_{i}$, such a scheme writes

[TABLE]

where the gradient $[v\cdot\nabla_{x}f]^{n}_{i,j}$ is computed using an upwind scheme

[TABLE]

The values $\varphi^{n}_{i}$ and $\eta^{n}_{i,j}$ are then deduced from the values of $f^{n}_{i,j},\rho^{n}_{i}$

[TABLE]

For large values of $\varepsilon$, the scheme (42) is accurate, since it is an order $1$ upwind scheme for an equation which contains no stiff terms. However, the computation of $\varphi$ and $\eta$ for small values of $\varepsilon$ with this method requires the explicit resolution of a stiff equation, which is costly from a computational point of view. Indeed, as a layer of width $\varepsilon$ appears in the asymptotic regime, a condition linking $\Delta x$ and $\varepsilon$ must be satisfied to ensure the solution of (42) to be accurate. Roughly, it writes $\Delta x\leq c\varepsilon$. In addition, the CFL condition

[TABLE]

makes the numerical resolution of (42) costly when $\varepsilon$ is small. The solution of (23) will be compared to the numerical solution of the limit equation (8)-(9) computed with (26)-(27). In this limit scheme, the numerical resolution of the nonlinear equation needed to find $H^{n+1}$ is performed with the Newton method proposed in [31]. In this method, the Newton method is used to solve (28). It is however slightly modified to ensure $\varphi^{n+1}_{i}$, defined in (26)-(27), and the denominator of (28), to remain nonnegative along the iterations. When considering the micro-macro scheme (23) for small $\varepsilon>0$, this constraint is automatically satisfied, thanks to the maximum principle.

6.1 Consistency and AP character of the scheme in the linear case $r=0$

In this section, we test the consistency and the AP property of the scheme (23) when $r=0$, and we consider spatially periodic boundary conditions. Firstly, we consider the initial value

[TABLE]

To check the consistency of the scheme (23), we test it for large values of $\varepsilon$, for which the solution of (11) is far from the solution of the limit equation (8)-(9). For $x\in[-1,1]$, we compute the solutions given by (23) for $\varepsilon=1,10^{-1}$, and $10^{-2}$, and the parameters $\Delta t=2.5\cdot 10^{-3}$, and $\Delta x=10^{-2}$. For $\varepsilon=1$ and $\varepsilon=10^{-1}$, the equation (4) can be solved numerically with (42) and the same numerical parameters. The solutions given by (23) are displayed in Fig. 1 and Fig. 2, they match with the solutions given by (42) which are represented on the same figures.

When $\varepsilon$ becomes smaller, the computation of the solution of (4) with the scheme (42) becomes costly, because of the stiffness of the problem, and of the condition (44). The solution of (23) for $\varepsilon=10^{-2}$ is then compared to the solution of the limit scheme (8) in Fig. 3. Once again, both are computed with $\Delta t=2.5\cdot 10^{-3}$ and $\Delta x=10^{-2}$. We observe that for such a small value of $\varepsilon$, the solution of the kinetic problem (11), computed with (23) is close to the solution of the asymptotic problem (8). These tests, done independently for $\varepsilon=1,10^{-1}$ and for $\varepsilon=10^{-2}$, show that the scheme (23) is accurate for the values of $\varepsilon$ for wich the problem is not stiff, and that it catches the viscosity solution of the limit problem (8) when $\varepsilon$ is smaller.

In both cases $r=0$ and $r>0$, the solution of (11) converges, when $\varepsilon$ goes to [math], to the viscosity solution of an Hamilton-Jacobi equation (8)-(9). Depending on the properties of the initial data, the solutions of these equations may not be smooth. Indeed, if there are two distinct minimum points in the initial data, the solutions of (8)-(9) lack the $\mathcal{C}^{1}$ regularity after some time. Theorem 2 ensures that the limit schemes (26)-(27) catch the viscosity solution of (8)-(9). As a consequence, they are robust when such a lack of regularity appears during the computations. However, to ensure that a scheme for (11) enjoys the AP property, it is necessary to make sure that it is robust in the case of an initial data leading to a solution which does not enjoy the $\mathcal{C}^{1}$ regularity, for all the values of $\varepsilon$. In order to test the robustness of (23) in this case, we now consider an initial data with two minimum points. Once again, for large values of $\varepsilon$, the solution of (23) is compared to the solution of (42), to check the agreement between them. The values $\varepsilon=1$ and $\varepsilon=10^{-1}$ are considered in Fig. 4 and Fig. 5, with the parameters $\Delta t=10^{-3}$ and $\Delta x=4\cdot 10^{-3}$.

When $\varepsilon$ is smaller, the solution of the scheme (23) is compared to the solution of (8)-(9). Still because of the stiffness of the problem (4) for small $\varepsilon$, when $\varepsilon$ is smaller it is costly to compare the solution of the micro-macro scheme (23) to the solution of the explicit scheme (42). In addition, it is not robust when the solution is not regular enough, as in the test case we are considering. Instead, we compare the results given by the micro-macro scheme (23) to the solutions given by the limit scheme (26), which catches the viscosity solution of the limit problem (8). The computations are done with the parameters $\Delta t=10^{-3}$, and $\Delta x=4\cdot 10^{-3}$, and the results are displayed in Fig. 6 for $\varepsilon=10^{-2}$ and Fig. 7 for $\varepsilon=10^{-3}$.

6.2 Consistency and AP character in the non-linear case $r>0$

In this section, we study the consistency and the AP property of the micro-macro scheme (23) in the non-linear case $r>0$. As in the linear case, the accuracy of the micro-macro scheme (23) is tested for different values of $\varepsilon$. The results are displayed in Fig. 8, where the phase $\varphi$ given by the scheme is compared to the solution of an explicit scheme (42) for (1) when $\varepsilon=1$ or $\varepsilon=10^{-1}$, and to the solution of the limit scheme (27) for $\varepsilon=10^{-2}$ or $\varepsilon=10^{-3}$. All the results are presented at time $T=0.5$, and have been computed with $\Delta x=10^{-2}$, and $\Delta t=2.5\cdot 10^{-3}$, and spatial periodic boundary conditions. As in the linear case, the comparison with the explicit scheme test shows that the scheme is accurate for large values of $\varepsilon$, and the comparison with the limit scheme highlights its AP character.

When $r>0$, it is expected that the nullset of $\varphi$ spreads at a fixed speed $c^{*}$, which depends on the model. Indeed, the following result has been proved in [6],

Proposition 4**.**

Assume that

[TABLE]

and define

[TABLE]

Then the nullset of the solution $\varphi$ of (9) propagates at speed $c^{*}$:

[TABLE]

The tests we propose in what follows are designed to highlight the propagation of fronts that is expected, thanks to the positivity of $r$. The speed of the propagation of the fronts is also tested at the numerical level. The initial condition we consider is a density $\rho$ such that $\rho=1$ at the left of the spatial domain and [math] elsewhere, with Neumann spatial boundary conditions. To ensure the test is done in the asymptotic regime, we consider $\varepsilon=10^{-4}$ with the parameters $\Delta t=3.125\cdot 10^{-4}$, and $\Delta x=1.25\cdot 10^{-3}$. Such a refined grid is not necessary to observe the propagation of the front, but it provides a better accuracy on the computation of the numerical propagation speed. The density given by the micro-macro scheme (23) at different times is displayed in Fig. 9. The speed of the propagation of the front is computed in the left part of Fig. 10. To determine numerically its theoretical value, the minimum of

[TABLE]

for $p>0$ is computed. It is presented in the right side of Fig. 10. The precision of the result depends on $\Delta x$, as highlighted in Fig. 11. In this figure, the relative error $|c_{\Delta x}^{*}-c^{*}|/c^{*}$, where $c^{*}_{\Delta x}$ denotes the numerical propagation speed, is presented as a function of $\Delta x$. In a logarithmic scale, a line is obtained, with slope $2$.

6.3 Order and uniform accuracy

In this section, we study the order of the scheme when $\varepsilon$ is fixed, and we investigate its uniform accuracy. To study the order of the scheme, we choose the initial data (45), and we define a reference solution as the solution given by (23)-(25) for $\Delta x=2\cdot 10^{-3}$, $\Delta t=5\cdot 10^{-4}$, and $\Delta v=1.25\cdot 10^{-2}$. Keeping $\Delta t$ and $\Delta v$ fixed, we compute the solutions of (23) for different values of $N_{x}$. The error

[TABLE]

is then computed for $\varepsilon=1$. It is displayed in Fig. 12 in logarithmic scale. As the slope of the line is slightly greater than $1$, the scheme is of order $1$ in space, as expected for an upwind scheme.

The scheme enjoys the UA property, with uniform order $1$ if

[TABLE]

with $C$ independent of $\varepsilon$. This property is is highlighed in Fig. 13 where the error $E(\varepsilon,\Delta x)$ is displayed in function of $\varepsilon$ for different values of $\Delta x$. As the lines are stratified, the scheme (23)-(25) is uniformly accurate in $\varepsilon$.

6.4 Case of a singular equilibrium

In this section, we consider the linear case $r=0$, and an equilibrium $M$ vanishing at some point of the velocity domain. It has been proved in [11] that, in this case, the asymptotic model is still a Hamilton-Jacobi equation, but that the Hamiltonian $H$ is not given by the implicit formula (10) anymore. Indeed, if we suppose that the velocity space $V$, which is the support of $M$, writes $V=[-v_{\max},v_{\max}]$ (such an hypothesis simplifies the notations but is not necessary), and if we denote

[TABLE]

and

[TABLE]

the following result holds

Theorem 3**.**

Suppose that $\varphi^{\varepsilon}(0,x,v)=\varphi_{\mathrm{in}}(x)$, then $\varphi^{\varepsilon}$ converges locally uniformly on $\mathbb{R}_{+}\times\mathbb{R}\times V$ towards $\varphi$, where $\varphi$ does not depend on $v$. Moreover, $\varphi$ is the viscosity solution of the following Hamilton-Jacobi equation:

[TABLE]

where the Hamiltonian $H$ is defined as follows: if $p\in\mathrm{Sing}(M)$, then $H(p)=\mu(p)-1$, else $H(p)$ is uniquely determined by the following implicit formula

[TABLE]

This theorem, which has been proved in [11] can be explained formally. Indeed, with the previous notations

[TABLE]

the equation (1) reads, for $r=0$

[TABLE]

which can be reformulated as follows

[TABLE]

where $H^{\varepsilon}=-\partial_{t}\varphi^{\varepsilon}$ and $p^{\varepsilon}=\partial_{x}\varphi^{\varepsilon}$. Since $\left\langle M\mathrm{e}^{-\eta^{\varepsilon}/\varepsilon}\right\rangle=1$, when $\varepsilon\to 0$, $Q^{\varepsilon}=\mathrm{e}^{-\eta^{\varepsilon}/\varepsilon}$ converges formally to the solution $Q$ of the spectral problem

[TABLE]

where $H=-\partial_{t}\varphi$, $p=\partial_{x}\varphi$. It implies that $1+H(p)-vp\geq 0$, and then that $H(p)\geq\mu(p)-1$. We distinguish between two cases. The first one is the same we treated in this paper, namely if $p\in\mathrm{Sing}(M)^{c}$, then

[TABLE]

by monotonicity of

[TABLE]

there exists $H(p)>\mu(p)-1$ such that

[TABLE]

But if $p\in\mathrm{Sing}(M)$, such an $H(p)$ cannot be determined and we have

[TABLE]

Therefore $H(p)$ is continuous but not $\mathcal{C}^{1}$ at $p$ such that

[TABLE]

In addition, the solution of the spectral problem (51) is singular, and a Dirac mass arises at $\pm v_{\max}$. It is even possible to compute the weight of this Dirac mass, see [11, 7]. Note that when the equilibrium satisfies (3), $\text{Sing}(M)=\emptyset$, and that in the case we are considering in this part

[TABLE]

$\text{Sing}(M)$ is not empty. Indeed, $M$ vanishes at $\pm(1-\Delta v/2)$ and Taylor expansions of $M$ at $v=\pm(1-\Delta v/2)$, show that any $p$ large enough belongs to $\text{Sing}(M)$.

In this part, we investigate the behavior of the micro-macro scheme with the equilibrium (52), since it was a priori not designed for such a singular equilibrium function. The consistency of the scheme with the kinetic equation for large values of $\varepsilon\sim 1$ is still clear, since it is written as a finite differences scheme for the kinetic equation. To test the accuracy of the scheme in the asymptotic regime, we first ensure that the reference scheme we use for the limit equation, which has been proposed in [31], is consistent with the limit equation defined with the constrained Hamiltonian (50). The right part of Fig. 14 presents the numerical Hamiltonian $H(p)$ computed with the limit scheme. This limit scheme is defined following the idea of [31] with an Euler scheme for (49), and the constrained Hamiltonian (50) computed with a constrained Newton’s method. As expected, it lacks the $\mathcal{C}^{1}$ character and the condition $H(p)\geq\mu(p)-1$ is satisfied. The left part of Fig. 14 presents the equilibrium function $M$ (52), with $m$ such that $\left\langle M\right\rangle_{N_{v}}=1$.

The AP character of the micro-macro scheme with this equilibrium is tested in Fig 15, where we fixed $\varepsilon=10^{-4}$, $\Delta x=10^{-2},\Delta t=2.5\cdot 10^{-3}$, and $\Delta v=5\cdot 10^{-2}$. It presents $\varphi$, the Hopf-Cole transform of the spatial density $\rho$, given by the micro macro scheme, and by the limit scheme. As the two curves match, the scheme seems to enjoy the AP property in the case of a singular equilibrium. It can actually be seen on the expression of the scheme. Indeed, since it can be written

[TABLE]

the third line implies that the denominator of the second line remains positive in the computations.

Let us now remark that, as it is stated in the continuous case in Thm. 3, $\eta^{n+1}_{i,j}$ vanishes when $\varepsilon$ goes to [math]. Indeed, if there are at least two index $j\in[\![1,N_{v}]\!]$ such that $M(v_{j})\neq 0$, the proof of the discrete maximum principle stated in Section 4 still holds. As a consequence, $\varphi^{n+1}_{i}$ and $\eta^{n+1}_{i,j}$ are uniformly bounded in $\varepsilon$. The third line of (53) provides the following inequality

[TABLE]

that can be recast as follows

[TABLE]

As a consequence, $\eta^{n+1}_{i,j}$ is nonnegative when $\varepsilon$ goes to [math] with the discretization parameters fixed. However, if it is positive, the third point of Prop. 3 cannot be satisfied in the small $\varepsilon$ limit. Eventually, $\eta^{n+1}_{i,j}$ vanishes when $\varepsilon$ goes to [math]. Thus, the following inequality holds when $\varepsilon\to 0$

[TABLE]

the constraint $H(\partial_{x}\varphi)\geq\mu(\partial_{x}\varphi)-1$ is hence fulfilled in the limit scheme.

Moreover, if $v=\pm v_{\max}$, the corrector $\mathrm{e}^{-\eta^{\varepsilon}(t,x,v)/\varepsilon}$ may not be bounded anymore, and a Dirac mass can even appear for such $v$, see [11]. This phenomena is highlighted by Fig. 16, in which the corrector $\mathrm{e}^{-\eta^{\varepsilon}(T,x,v)/\varepsilon}$ is represented at time $T=1$ as a function of $x$ and $v$ for different values of $\varepsilon$. As it is expected, we observe that a Dirac mass arises at some points located at the border of the velocity domain.

As a conclusion, the scheme presented here handles singular measures in the velocity variable which can arise when the equilibrium distribution $M$ vanishes at the boundary of $V$. However, we lack suitable numerical analysis. For instance, the accuracy of the resolution of the non-linear system (53) may be reduced when Dirac masses arise in the corrector.