Analysis of an asymptotic preserving scheme for linear kinetic equations in the diffusion limit

TL;DR

This paper provides a mathematical analysis demonstrating that an asymptotic preserving scheme for linear kinetic equations remains stable and accurate across different regimes, with conditions depending on the mean free path.

Contribution

The paper offers a rigorous proof of the uniform stability and accuracy of the scheme under a CFL condition that adapts to the mean free path, using simple energy estimates.

Findings

Scheme is uniformly stable and accurate across regimes.

CFL condition adapts to mean free path size.

Analysis uses simple energy estimates.

Abstract

We present a mathematical analysis of the asymptotic preserving scheme proposed in [M. Lemou and L. Mieussens, SIAM J. Sci. Comput., 31, pp. 334-368, 2008] for linear transport equations in kinetic and diffusive regimes. We prove that the scheme is uniformly stable and accurate with respect to the mean free path of the particles. This property is satisfied under an explicitly given CFL condition. This condition tends to a parabolic CFL condition for small mean free paths, and is close to a convection CFL condition for large mean free paths. Ou r analysis is based on very simple energy estimates.

Equations

$$\epsilon {\partial }_t\varphi +\Omega \cdot {\nabla }_r\varphi =\frac{\sigma }{\epsilon }L\varphi -\epsilon {\sigma }_A\varphi +\epsilon S,$$

$$\epsilon {\partial }_{t}f+v{\partial }_{x}f=\frac{\sigma }{\epsilon }Lf-\epsilon {\sigma }_{A}f+\epsilon S,$$

$$Lf(v)={\int }_{-1}^{1}s(v,v^{\prime })(f(v^{\prime })-f(v))d{v}^{\prime },$$

$$\left[\varphi L\varphi \right]\le -2s_m\left[{\varphi }^2\right]$$

$${\partial }_t\rho -{\partial }_x\kappa {\partial }_x\rho =-{\sigma }_A\rho +S,$$

$${\partial }_t\rho +{\partial }_x\left[vg\right]=-{\sigma }_A\rho +S,$$

$${\partial }_{t}g+\frac{1}{\epsilon }(I-\left[.\right])(v{\partial }_{x}g)=-\frac{\sigma }{{\epsilon }^2}Lg-\frac{1}{{\epsilon }^2}v{\partial }_x\rho ,$$

$$\frac{{\rho }_i^{n+1}-{\rho }_i^n}{\Delta t}+\left[v\frac{g_{i+\frac{1}{2}}^{n+1}-g_{i-\frac{1}{2}}^{n+1}}{\Delta x}\right]=-{\sigma }_{A,i}{\rho }_i^{n+1}+S_i,$$

$$\frac{g_{i+\frac{1}{2}}^{n+1}-g_{i+\frac{1}{2}}^n}{\Delta t}+\frac{1}{\epsilon \Delta x}(I-\left[.\right])\left(v^{+}(g_{i+\frac{1}{2}}^n-g_{i-\frac{1}{2}}^n)+v^{-}(g_{i+\frac{3}{2}}^n-g_{i+\frac{1}{2}}^n)\right)$$

$$=-\frac{{\sigma }_{i+\frac{1}{2}}}{{\epsilon }^2}L{g}_{i+\frac{1}{2}}^{n+1}-\frac{1}{{\epsilon }^2}v\frac{{\rho }_{i+1}^n-{\rho }_i^n}{\Delta x},$$

$$\left[g_{i+\frac{1}{2}}^n\right]=0.$$

$$\left[v{g}_{i+\frac{1}{2}}^{n+1}\right]=-\frac{1}{{\sigma }_{i+\frac{1}{2}}}\left[v{L}^{-1}v\right]\frac{{\rho }_{i+1}^n-{\rho }_i^n}{\Delta x}+O(\epsilon ),$$

$$\frac{{\rho }_i^{n+1}-{\rho }_i^n}{\Delta t}-\frac{1}{\Delta x}\left({\kappa }_{i+\frac{1}{2}}\frac{{\rho }_{i+1}^n-{\rho }_i^n}{\Delta x}-{\kappa }_{i-\frac{1}{2}}\frac{{\rho }_i^n-{\rho }_{i-1}^n}{\Delta x}\right)=-{\sigma }_{A,i}{\rho }_i^n+S_i,$$

$$\Delta t\le \frac{1}{3}\left(\tilde{\sigma }\Delta x^2+2\epsilon \Delta x\right),$$

$$\sum _{i\in \mathbb{Z}}{\left({\rho }_i^n\right)}^2\Delta x+{\epsilon }^2\sum _{i\in \mathbb{Z}}\left[{\left(g_{i+\frac{1}{2}}^n\right)}^2\right]\Delta x\le \sum _{i\in \mathbb{Z}}{\left({\rho }_i^0\right)}^2\Delta x+{\epsilon }^2\sum _{i\in \mathbb{Z}}\left[{\left(g_{i+\frac{1}{2}}^0\right)}^2\right]\Delta x$$

$$\Delta t\le \min \left(\frac{2}{1+{\sigma }_{A,M}},\frac{3}{3+{\sigma }_{A,M}}\Delta t_S\right),$$

$$\Vert \mu {\Vert }^2=\sum _{i\in \mathbb{Z}}{\mu }_i^2\Delta x.$$

$$\Vert |\varphi \Vert |=\sum _{i\in \mathbb{Z}}\left[{\varphi }_{i+\frac{1}{2}}^2\right]\Delta x.$$

$$⟨\varphi ,\psi ⟩=\sum _{i\in \mathbb{Z}}\left[{\varphi }_{i+\frac{1}{2}}{\psi }_{i+\frac{1}{2}}\right]\Delta x.$$

$$D^{-}{\varphi }_{i+\frac{1}{2}}=\frac{{\varphi }_{i+\frac{1}{2}}-{\varphi }_{i-\frac{1}{2}}}{\Delta x}\text{ and }{D}^{+}{\varphi }_{i+\frac{1}{2}}=\frac{{\varphi }_{i+\frac{3}{2}}-{\varphi }_{i+\frac{1}{2}}}{\Delta x}$$

$$D^c{\varphi }_{i+\frac{1}{2}}=\frac{{\varphi }_{i+\frac{3}{2}}-{\varphi }_{i-\frac{1}{2}}}{2\Delta x}\text{ and }{D}^0{\varphi }_i=\frac{{\varphi }_{i+\frac{1}{2}}-{\varphi }_{i-\frac{1}{2}}}{\Delta x}(=D^{-}{\varphi }_{i+\frac{1}{2}}).$$

$${\delta }^0{\mu }_{i+\frac{1}{2}}=\frac{{\mu }_{i+1}-{\mu }_i}{\Delta x}.$$

$$\left(v^{+}{D}^{-}+v^{-}{D}^{+}\right){\varphi }_{i+\frac{1}{2}}=v{D}^c{\varphi }_{i+\frac{1}{2}}-\frac{\Delta x}{2}|v|D^{-}{D}^{+}{\varphi }_{i+\frac{1}{2}}.$$

$$\sum _{i\in \mathbb{Z}}{\left(D^{+}{\varphi }_{i+\frac{1}{2}}\right)}^2\Delta x\le \frac{4}{\Delta x^2}\sum _i{\varphi }_{i+\frac{1}{2}}^2\Delta x.$$

$$\left|⟨\left(v^{+}{D}^{+}+v^{-}{D}^{-}\right)\psi ,\varphi ⟩\right|\le \alpha \Vert |\varphi \Vert {|}^2+\frac{1}{4\alpha }\Vert ||v|D^{+}\psi \Vert {|}^2.$$

$$\left|⟨\left(v^{+}{D}^{+}+v^{-}{D}^{-}\right)\psi ,\varphi ⟩\right|\le \alpha \Vert |\varphi \Vert {|}^2+\frac{1}{4\alpha }\Vert |\left(v^{+}{D}^{+}+v^{-}{D}^{-}\right)\psi \Vert {|}^2.$$

$$\frac{1}{4\alpha }\Vert |\left(v^{+}{D}^{+}+v^{-}{D}^{-}\right)\psi \Vert {|}^2=\frac{1}{4\alpha }\sum _{i\in \mathbb{Z}}\frac{1}{2}\left({\int }_{-1}^0{v}^2(D^{-}{\psi }_{i+\frac{1}{2}}{)}^{2}dv+{\int }_0^1{v}^2(D^{+}{\psi }_{i+\frac{1}{2}}{)}^{2}dv\right)\Delta x.$$

$$\begin{array}{rl}\frac{1}{4\alpha }\Vert |\left(v^{+}{D}^{+}+v^{-}{D}^{-}\right)\psi \Vert {|}^2& =\frac{1}{4\alpha }\sum _{i\in \mathbb{Z}}\frac{1}{2}{\int }_{-1}^1{v}^2(D^{+}{\psi }_{i+\frac{1}{2}}{)}^{2}dv\Delta x\\ & =\frac{1}{4\alpha }\Vert ||v|D^{+}\psi \Vert {|}^2.\end{array}$$

$$\displaystyle\sum_{i\in{\mathbb{Z}}}\mu_{i}D^{0}\phi_{i}\Delta x=-\sum_{i\in{\mathbb{Z}}}\bigl{(}\delta^{0}\mu_{i+\frac{1}{2}}\bigr{)}\phi_{i+\frac{1}{2}}\Delta x,$$

$$\displaystyle\sum_{i\in{\mathbb{Z}}}\psi_{i+\frac{1}{2}}D^{-}\phi_{i+\frac{1}{2}}\Delta x=-\sum_{i\in{\mathbb{Z}}}\bigl{(}D^{+}\psi_{i+\frac{1}{2}}\bigr{)}\phi_{i+\frac{1}{2}}\Delta x,$$

$$\sum _{i\in \mathbb{Z}}{\varphi }_{i+\frac{1}{2}}{D}^c{\varphi }_{i+\frac{1}{2}}\Delta x=0.$$

$${\left[vg\right]}^2\le \frac{1}{2}\left[|v|g^2\right]$$

$$\begin{array}{rl}{\left[vg\right]}^2& =\frac{1}{4}{\left({\int }_{-1}^{1}vgdv\right)}^2=\frac{1}{4}{\left({\int }_{-1}^1\mathrm{sign}(v)\sqrt{|v|}\sqrt{|v|}gdv\right)}^2\\ & \le \frac{1}{4}{\int }_{-1}^1{\left(\mathrm{sign}(v)\sqrt{|v|}\right)}^{2}dv{\int }_{-1}^1{\left(\sqrt{|v|}g\right)}^{2}dv\text{ by Cauchy-Schwarz inequality}\\ & =\frac{1}{2}\left[|v|g^2\right].\end{array}$$

Full text

*draft, for submission to SINUM, *

**Analysis of an asymptotic preserving scheme for linear kinetic equations in the diffusion limit

**

Jian-Guo Liu111Department of Physics and Department of Mathematics, Duke University, Durham, NC 27708, USA (Jian-Guo.Liu@duke.edu), Luc Mieussens222Institut de Mathématiques de Bordeaux (UMR 5251), Université de Bordeaux, 351, cours de la Libération, 33405 Talence cedex, France (Luc.Mieussens@math.u-bordeaux1.fr)

**Abstract. ** We present a mathematical analysis of the asymptotic preserving scheme proposed in [M. Lemou and L. Mieussens, SIAM J. Sci. Comput., 31, pp. 334–368, 2008] for linear transport equations in kinetic and diffusive regimes. We prove that the scheme is uniformly stable and accurate with respect to the mean free path of the particles. This property is satisfied under an explicitly given CFL condition. This condition tends to a parabolic CFL condition for small mean free paths, and is close to a convection CFL condition for large mean free paths. Our analysis is based on very simple energy estimates.

Key words. transport equations, diffusion limit, asymptotic preserving schemes, stiff terms, stability analysis

AMS subject classifications. 65M06, 35B25, 82C80, 41A60

1 Introduction

Particle systems are often described at the microscopic level by kinetic models (neutron transport, radiative transfer, electrons in semi-conductors, or rarefied gas dynamics). Simulating such systems by using kinetic models can be computationally very expensive, but modern super computers now enable realistic simulations. When the mean free path of the particles is very small as compared to the (macroscopic) size of the computational domain, the kinetic model can be very well approximated by a much simpler macroscopic model (diffusion equation, Rosseland approximation, Euler and Navier-Stokes equations), that can be numerically solved much faster.

However, there are many cases where the ratio mean free path/macroscopic size (the so-called “Knudsen number” in rarefied gas dynamics, denoted by $\varepsilon$ in this paper) is not constant: depending on the geometry of the boundaries, or on the boundary conditions, this ratio may vary in time, and in space. In such multiscale situations, usual kinetic solvers are often useless: for stability and accuracy reasons, they must resolve the microscopic scales, which is computationally to expensive in “fluid” zones (where $\varepsilon$ is small). By contrast, macroscopic solvers are faster but may be inaccurate in “kinetic” zones (where $\varepsilon$ is large).

This is why many people have been working for more than 20 years on a kind of multiscale kinetic schemes: the asymptotic-preserving (AP) schemes. Such schemes are uniformly stable with respect to $\varepsilon$ (thus their computational complexity does not depend on $\varepsilon$), and are consistent with the macroscopic model when $\varepsilon$ goes to [math] (the limit of the scheme is a scheme for the macroscopic model).

Up to our knowledge, AP schemes have first been studied (for steady problems) in neutron transport by Larsen, Morel and Miller [22], Larsen and Morel [21], and then by Jin and Levermore [10, 11]. For unstationary problems, the difficulty is the time stiffness due to the collision operator. To avoid the use of expensive fully implicit schemes, two classes of semi-implicit time discretizations have been proposed by Klar [16] and Jin, Pareschi and Toscani [15] (see preliminary works in [14, 9] and extensions in [13, 12, 25, 17, 18]. Similar ideas have also been used by Buet et al. in [4].

In [23], Lemou and Mieussens have proposed a new AP scheme based on the micro-macro decomposition of the distribution function into microscopic and macroscopic components (similar schemes have also been proposed by Klar and Schmeiser [19], and more recently by Carrillo et al. [5, 6]). A coupled system of equations is obtained for these two components without any linearity assumption. The decomposition only uses basic properties of the collision operator that are common to most of kinetic equations (namely conservation and equilibrium properties). Then this system is solved with a suitable time semi-implicit discretization and space finite differences on staggered grids. While almost all the schemes mentioned before are based on very similar ideas, the approach proposed in [23] has been shown to be very general, since it applies to kinetic equations for both diffusion and hydrodynamic regimes (for the diffusion regime, see [23] for linear transport equations and [3] for the non-linear Kac equation, for the hydrodynamic regime, see [2] for the Boltzmann equation). We also mention the work of Degond, Liu and Mieussens [7] who proposed a similar approach (micro-macro decomposition) to design macroscopic diffusion models with kinetic upscalings: this approach also leads to AP schemes, at least for a semi-discrete time discretization (no AP space discretization was studied in this paper).

While many different schemes have been proposed in the past few years, it appears that the rigorous proof of their AP property looks rather difficult and is seldom investigated. Up to our knowledge, there are only two papers on this subject. Klar and Unterreiter [20] have proved that a scheme similar to that of [16] and [15] for the linear transport equation is uniformly stable. However, their proof is based on a von Neumann analysis, and hence is restricted to a one dimensional equation with constant coefficients in a periodic domain. Gosse and Toscani have proposed in [8] an AP scheme based on a rather different idea (discretization of the collision term as a non-conservative product and use of well-balanced Godunov schemes): they have been able to prove a uniform stability property, still in the linear case, but also a strong positivity property. However, their scheme is based on techniques (like approximation of the steady solution) that are difficult to generalize to other equations.

In this paper, we propose a very simple stability proof for the AP scheme of [23] and we exhibit an explicit CFL condition. This condition is uniform with respect to $\varepsilon$, and gives a standard parabolic CFL condition $\Delta t=O(\Delta x^{2})$ when $\varepsilon$ goes to 0. While a stability property was already proved in [23] for this scheme, this was only for a simple two-velocity model (telegraph equation), by using von Neumann analysis. Here, our proof is based on energy estimates that are more general than von Neumann analysis, and hence is valid for one-dimensional linear equations with non-constant coefficients, continuous velocity variable, in the whole space. By the same technique, we are also able to prove uniform error estimates.

Our paper is organized as follows. In section 2, we introduce a general linear equation, we present its discretization by the AP scheme, and we give the main features of this scheme. Then, in section 3, we give our main stability result and its proof. The error estimates are given in section 4.

2 An AP scheme for the linear transport equation

Linear transport equation is a model for the evolution of particles in some medium (neutron transport, linear radiative transfer, …). Generally, this model reads, in scaled variables

[TABLE]

where $\phi(t,r,\Omega)$ is the number density of particles in the position-direction phase space that depends on time $t$, position $r=(x,y,z)\in{\mathbb{R}}^{3}$, and angular direction of propagation of particles $\Omega\in S^{2}$. Moreover, $\sigma$ is the total cross section, $\sigma_{A}$ is the absorption cross section, and $S$ is an internal source of particles, which is independent of $\Omega$. The linear operator $L$ models the scattering of the particles by the medium and acts only on the angular dependence of $\phi$. This simple model does not allow for particles of possibly different energy (or frequency); it is called “one-group” or “monoenergetic” equation. The parameter $\varepsilon$ is a scale factor that measures the ratio between a typical microscopic length (the mean free path of a particle, for instance) to a typical macroscopic length (the size of the computational domain, for instance). See [22] for details.

In this paper, we consider this one-group equation in the slab geometry: we assume that $\phi$ depends only on the slab axis variable $x\in{\mathbb{R}}$. Then it can be shown that the average of $\phi$ with respect to the $(y,z)$ cosine directions of $\Omega$, denoted by $f(t,x,v)$, satisfies the one-dimensional equation

[TABLE]

where $v\in[-1,1]$ is the $x$ cosine direction of $\Omega$. At $t=0$, we have the initial data $f(0,x,v)=f^{0}(x,v)$. We assume that the cross sections satisfy the inequalities $0<\sigma_{m}\leq\sigma(x)\leq\sigma_{M}$ and $0\leq\sigma_{A}(x)\leq\sigma_{AM}$ for every $x$. We do not consider boundary conditions in this paper, and hence $x\in{\mathbb{R}}$.

The linear operator $L$ is given by

[TABLE]

where the kernel $s$ is such that $0<s_{m}\leq s(v,v^{\prime})\leq s_{M}$ for every $v,v^{\prime}$ in $[-1,1]$. We also assume that $s$ satisfies $\int_{-1}^{1}s(v,v^{\prime})\,dv^{\prime}=1$, and that it is symmetric: $s(v,v^{\prime})=s(v,v^{\prime})$. For the sequel, it is useful to define the operator $\left[.\right]$ such that $\left[\phi\right]=\frac{1}{2}\int_{-1}^{1}\phi(v)\,dv$ is the average of every velocity dependent function $\phi$. With these assumptions, it is standard [1] to state the following properties of $L$:

Proposition 2.1**.**

• •

$\left[L\phi\right]=0$* for every $\phi$ in $L^{2}([-1,1])$*

• •

the null space of $L$ is ${\cal N}(L)=\{\phi=\left[\phi\right]\}$ (constant functions)

• •

the rank of $L$ is ${\cal R}(L)={\cal N}^{\perp}(L)=\{\phi\text{ s.t }\left[\phi\right]=0\}$

• •

$L$* is non-positive self-adjoint in $L^{2}([-1,1])$ and we have*

[TABLE]

for every $\phi\in{\cal N}^{\perp}(L)$

• •

$L$* admits a pseudo-inverse from ${\cal N}^{\perp}(L)$ onto ${\cal N}^{\perp}(L)$, denoted by $L^{-1}$*

• •

the orthogonal projection from $L^{2}([-1,1])$ onto ${\cal N}^{\perp}(L)$ is $\left[.\right]$

When $\varepsilon$ becomes small (the “diffusion” regime), it is well known that the solution $f$ of (1) tends to its own average density $\rho=\left[f\right]$, which is a solution of the asymptotic diffusion limit

[TABLE]

where the diffusion coefficient is $\kappa(x)=-\frac{\left[vL^{-1}v\right]}{\sigma(x)}$. An asymptotic preserving scheme for the linear kinetic equation (1) is a numerical scheme that discretizes (1) in such a way that it leads to a correct discretization of the diffusion limit (4) when $\varepsilon$ is small.

Now, we summarize the results obtained in [23]. By using the micro-macro decomposition $f=\rho+\varepsilon g$, where $\rho=\left[f\right]$ and $g$ is such that $\left[g\right]=0$, we derived the micro-macro model for (1) that reads

[TABLE]

with initial data $\rho^{0}=\left[f^{0}\right]$ and $\varepsilon g^{0}=f^{0}-\rho^{0}$. This system if formally equivalent to (1).

Then, we proposed the following numerical scheme for this system. We choose a time step $\Delta t$ and times $t_{n}=n\Delta t$, and two staggered grids of step $\Delta x$ and nodes $x_{i}=i\Delta x$ and $x_{i+\frac{1}{2}}=(i+\frac{1}{2})\Delta x.$ We use the approximated values $\rho^{n}_{i}\approx\rho(t_{n},x_{i})$ and $g^{n}_{i+\frac{1}{2}}(v)\approx g(t_{n},x_{i+\frac{1}{2}},v)$, and the scheme reads

[TABLE]

where $v^{\pm}=\frac{v\pm|v|}{2}$. Note that as in the continuous case, this scheme preserves the zero average of $g$:

Proposition 2.2**.**

If the initial data $g^{0}$ satisfies $\left[g^{0}_{i+\frac{1}{2}}\right]=0$ for every $i$, then for every $n$ and $i$, we have

[TABLE]

Proof.

Apply the average operator $\left[.\right]$ to (6b): since we have $\left[I-\left[.\right]\right]=0$ (obvious), $\left[L\right]=0$ (proposition 2.1), and $\left[v\right]=0$ (obvious), this relation yields $\left[g^{n+1}_{i+\frac{1}{2}}\right]-\left[g^{n}_{i+\frac{1}{2}}\right]=0$, which gives the result. ∎

In scheme (6), the upwind discretization of $(I-\left[.\right])(v\partial_{x}g)$ is to insure stability in the kinetic regime, while the centered approximations of $\partial_{x}\left[vg\right]$ and $v\partial_{x}\rho$ are to capture the diffusion limit. Indeed, it is clear that when $\varepsilon$ goes to 0, we have from (6b) $g^{n+1}_{i+\frac{1}{2}}=-\frac{1}{\sigma_{i+\frac{1}{2}}}L^{-1}\left(v\frac{\rho^{n}_{i+1}-\rho^{n}_{i}}{\Delta x}\right)+O(\varepsilon)$. Consequently, the flux of $g^{n+1}_{i+\frac{1}{2}}$ is

[TABLE]

where we have used the facts that $\rho^{n}_{i}$ does not depend on $v$ and that $L$—and hence $L^{-1}$ too—only applies to functions of $v$. Then, using this relation in (6a), we get

[TABLE]

with $\kappa_{i+\frac{1}{2}}=-\frac{\left[vL^{-1}v\right]}{\sigma_{i+\frac{1}{2}}}$, which is the usual 3-points stencil explicit scheme for the diffusion equation (4).

Of course, this property is true only if $\Delta t$ can be chosen independently of $\varepsilon$, or in other words, if the scheme is uniformly stable with respect to $\varepsilon$. In [23], we have proved that this scheme is indeed uniformly stable under some CFL condition, in the simpler case of the telegraph equation (in which we have only two discrete velocities $v=\pm 1$, and $s\equiv 1$, $\sigma\equiv 1$, $\sigma_{A}=S\equiv 0$). In the following section, we extend this result to the general equation (1) for scheme (6).

3 Uniform stability

We give our main result of stability for scheme (6). Without lost of generality, we assume $S=0$ (source free case).

Theorem 3.1**.**

If $\Delta t$ satisfies the following CFL condition

[TABLE]

with $\tilde{\sigma}=2s_{m}\sigma_{m}$, then the sequences $\rho^{n}$ and $g^{n}$ defined by scheme (6) satisfy the energy estimate

[TABLE]

for every $n$, and hence the scheme (6) is stable.

Note that CFL condition (9) can be viewed as an average of a diffusive CFL condition $\Delta t\leq\tilde{\sigma}\Delta x^{2}$ (needed for the diffusion scheme (8)) and of a convection CFL $\Delta t\leq\varepsilon\Delta x$. It shows that the scheme is stable uniformly in $\varepsilon$, that is to say a diffusive CFL condition $\Delta t\leq C\Delta x^{2}$ is sufficient for stability for small $\varepsilon$, while a convection CFL is sufficient for $\varepsilon=O(1)$.

Moreover, while our CFL condition (9) is valid for every $\varepsilon$, we point out that it is not optimal for each $\varepsilon$. This might be the price to be paid for getting a uniform condition. For instance, in the simpler case of constant cross sections ($\sigma\equiv\sigma_{m}$) and kernel ($s=s_{m}\equiv\frac{1}{2}$), the diffusion coefficient of the diffusion equation (8) is $\kappa=\frac{1}{3\sigma}$, and the optimal CFL condition for scheme (8) is $\Delta t\leq\frac{\Delta x^{2}}{2\kappa}$. However, CFL condition (9) for our AP scheme reads (for $\varepsilon=0$) $\Delta t\leq\frac{2}{9}\frac{\Delta x^{2}}{2\kappa}$, which is $0.2$ as small as the optimal CFL.

Remark 3.1*.*

If we have a source term $S\neq 0$, then the same CFL condition naturally gives a linear growth of the energy. Since this is standard and does not lead to any additional difficulty, we do not consider this case here.

Remark 3.2*.*

In the case of a time explicit discretization of the absorption term (that is to say when $-\sigma_{A,i}\rho^{n+1}_{i}$ is replaced by $-\sigma_{A,i}\rho^{n}_{i}$ in (6), we have the same result with a CFL condition which is a bit more complicated. Namely, the scheme is stable if

[TABLE]

where $\Delta t_{S}$ is the maximum time step allowed by CFL condition (9) for the scheme with implicit (or zero) absorption term.

This theorem is proved in the following sections: section 3.1 contains compact notations and useful lemma, and section 3.2 contains the derivation of our energy estimate.

3.1 Notations and useful lemma

We give some useful notations for norms and inner products that are used in our analysis. For every grid function $\mu=(\mu_{i})_{i\in{\mathbb{Z}}}$ we define:

[TABLE]

For every velocity dependent grid function $v\in[-1,1]\mapsto\phi(v)=(\phi_{i+\frac{1}{2}}(v))_{i\in{\mathbb{Z}}}$, we define:

[TABLE]

If $\phi$ and $\psi$ are two velocity dependent grid functions, we define their inner product:

[TABLE]

Now we give some notations for the finite difference operators that are used in scheme (6). For every grid function $\phi=(\phi_{i+\frac{1}{2}})_{i\in{\mathbb{Z}}}$, we define the following one-sided operators:

[TABLE]

We also define the following centered operators:

[TABLE]

Finally, for every grid function $\mu=(\mu_{i})_{i\in{\mathbb{Z}}}$, we define the following centered operator:

[TABLE]

Lemma 3.1** (Centered form of the upwind operator).**

For every grid function $\phi=(\phi_{i+\frac{1}{2}})_{i\in{\mathbb{Z}}}$, we have:

[TABLE]

Proof.

This result is easily obtained by using the relations $v^{\pm}=\frac{v\pm|v|}{2}$ and the identity $D^{-}+D^{+}=2D^{c}$. ∎

Lemma 3.2**.**

For every grid function $\phi=(\phi_{i+\frac{1}{2}})_{i\in{\mathbb{Z}}}$, we have:

[TABLE]

Proof.

Expand the left-hand side, use the inequality $(a+b)^{2}\leq 2(a^{2}+b^{2})$, and then use the change of index $i+\frac{3}{2}\rightarrow i+\frac{1}{2}$. ∎

Lemma 3.3** (Estimate for the adjoint upwind operator).**

For every positive real number $\alpha$ and for every velocity dependent grid functions $\phi$ and $\psi$, we have:

[TABLE]

Proof.

From Young’s inequality, we obtain for any positive real number $\alpha$:

[TABLE]

Now, the second term of the right-hand side of this inequality can be written:

[TABLE]

Then, with a simple changing of index, $D^{-}$ can be replaced by $D^{+}$ in the first integral, which gives:

[TABLE]

Finally, using this inequality in (17) gives the result. ∎

Lemma 3.4** (Discrete integration by parts).**

For every grid functions $\phi=(\phi_{i+\frac{1}{2}})_{i\in{\mathbb{Z}}}$, $\psi=(\psi_{i+\frac{1}{2}})_{i\in{\mathbb{Z}}}$, and $\mu=(\mu_{i})_{i\in{\mathbb{Z}}}$, we have:

[TABLE]

Proof.

These results are simply obtained by using obvious changing of indexes and the definitions of the finite difference operators given in (14–16). ∎

Lemma 3.5**.**

If $g\in L^{2}([-1,1])$ then

[TABLE]

Proof.

We just note that

[TABLE]

∎

3.2 Energy estimates

Since the time discretization of absorption term $-\sigma_{A}\rho$ is implicit, it plays no role in the energy estimate. Consequently, to simplify the proof, we consider the case where there is no absorption (see remark 3.3 at the end of this section). We proceed in five short steps.

Step 1.

Here we derive a first energy relation. With the finite difference operators defined in (14)–(16), the scheme can be written in the following compact form:

[TABLE]

We define the energy of system (5) as $\int_{\mathbb{R}}\rho^{2}\,dx+\varepsilon^{2}\int\left[g^{2}\right]\,dx$. It is clear that the scheme can be proved to be stable if the energy at time $n+1$ can be controlled by the energy at time $n$. Consequently, we first multiply (18a) by $\rho^{n+1}_{i}$, then we take the sum over $i\in{\mathbb{Z}}$, and finally, we use the standard equality $a(a-b)=\frac{1}{2}(a^{2}-b^{2}+|a-b|^{2})$ to get:

[TABLE]

Now we use relation (7): since $\left[g^{n+1}_{i+\frac{1}{2}}\right]=0$ for every $i$, a simple expansion of the inner product of the left-hand side of (19b) shows that it can be reduced to:

[TABLE]

Moreover, we can use relation (3) and the assumptions of the cross section and the kernel to estimate the inner product of the right-hand side of (19b) as follows:

[TABLE]

Consequently, we add up (19a) and $\varepsilon^{2}$ times (19b), then we use (20), (21), and the discrete integration by parts of lemma 3.4 to get our preliminary energy estimate:

[TABLE]

Step 2.

In this step, we show how the $\rho^{n+1}-\rho^{n}$ term can be eliminated in (22). First, it is useful to write $\rho^{n}_{i}$ in the right-hand side of (22) as $(\rho^{n}_{i}-\rho^{n+1}_{i})+\rho^{n+1}_{i}$: indeed, the terms $\sum_{i\in{\mathbb{Z}}}\rho^{n+1}_{i}D^{0}\left[vg^{n+1}_{i}\right]$ and $\sum_{i\in{\mathbb{Z}}}\left[vD^{0}g^{n+1}_{i}\right]\rho^{n+1}_{i}$ in the left and right-hand sides cancel out and we obtain:

[TABLE]

Now, we use the following Young inequality:

[TABLE]

Then the $\rho^{n+1}-\rho^{n}$ terms cancel out in (23) if $\alpha=\frac{1}{2\Delta t}$ and we get

[TABLE]

Step 3.

Here, we work on the inner product of (25) to show that the $g^{n+1}-g^{n}$ terms can also be eliminated. First, we insert $g^{n+1}$ in this inner product to get:

[TABLE]

and we rewrite terms $A$ and $B$ as follows. For $A$, we use the centered form of the upwind operator (lemma 3.1) and the discrete integration by parts of lemma 3.4 to get:

[TABLE]

For $B$, we also use the discrete integration by parts of lemma 3.4 to get:

[TABLE]

Then we apply the inequality of lemma 3.3 to $B$ to get

[TABLE]

Therefore, using (25), (26), (27) and (29), we see that the $g^{n+1}-g^{n}$ terms cancel out in (25) if $\alpha=\frac{\varepsilon}{2\Delta t}$, and we get

[TABLE]

Step 4.

Now, we show how all the $D^{+}g^{n+1}$ and the $D^{0}g^{n+1}$ terms can be controlled by $\||g^{n+1}\||$. First, note that the term $\frac{\Delta t}{2}\|||v|D^{+}g^{n+1}\||^{2}$ of the left-hand side of (30) can be estimated as follows:

[TABLE]

since $|v|\leq 1$. Moreover, using lemma 3.5 and a change of indices shows that the last term of the right-hand side of (30) satisfies

[TABLE]

Finally, we use these two estimates in (30) to obtain:

[TABLE]

Now, taking the positive part of the factor $(\frac{3\Delta t}{4}-\varepsilon\frac{\Delta x}{2})$ of the right-hand side of (33), we have the estimate

[TABLE]

where we have used $|v|\leq 1$ and the estimate of lemma 3.2.

Step 5.

Finally, estimates (33) and (34) show that

[TABLE]

This means that we have the final energy estimate

[TABLE]

if $\Delta t$ is such that

[TABLE]

Since $\tilde{\sigma}\geq 0$, an equivalent condition is $(\frac{3\Delta t}{4}-\varepsilon\frac{\Delta x}{2})\frac{4}{\Delta x^{2}}\leq\sigma$, which gives the sufficient condition

[TABLE]

which proves the theorem.

Remark 3.3*.*

As explained at the beginning of this section, when the absorption term is non zero and is discretized implicitly, its contribution $-\sum_{i\in{\mathbb{Z}}}\sigma_{A,i}(\rho^{n+1}_{i})^{2}\Delta x$ to the energy estimate is non-positive and plays no role in the previous analysis. However, if we use instead an explicit discretization, then our analysis has to be modified. The difference now is that there is the additional term $-\sum_{i\in{\mathbb{Z}}}\sigma_{A,i}\rho^{n+1}_{i}\rho^{n}_{i}\Delta x$ in the right-hand side of (22). The idea of step 2 can be applied to this term (replace $\rho^{n}$ by $(\rho^{n}-\rho^{n+1})+\rho^{n+1}$ and use a Young inequality) so that the $\rho^{n+1}-\rho^{n}$ terms cancel out, but now with $\alpha=\frac{1}{2\Delta t(1+\sigma_{A,M})}$. The other steps are the same, except that some coefficients are different. In order to shorten the paper, these details are left to the reader.

4 Error estimates

In this section, we simplify the presentation by taking constant total cross section $\sigma$ and kernel $s\equiv\frac{1}{2}$, no absorption ($\sigma_{A}\equiv 0$), and no source term ($S\equiv 0$). These assumptions are not restrictive at all, and our analysis could be directly applied in the general case.

Let $T>0$ be some finite time. For this study, we assume that the exact solution $(\rho,g)$ of (5) has the following regularity:

[TABLE]

Note that this assumption implies that $g$ is uniformly bounded with respect to $\varepsilon$, which is quite strong. In particular, this excludes the case of initial or transition layers. Indeed, if, for instance, $f$ is not isotropic at $t=0$, then $\rho(t=0)\neq f(t=0)$, and hence $g(t=0)=\frac{1}{\varepsilon}(f-\rho)|_{t=0}=O(\varepsilon^{-1})$ cannot be uniformly bounded with respect to $\varepsilon$.

With this assumption, we can obtain the following result.

Theorem 4.1**.**

If $\Delta t$ satisfies the following condition

[TABLE]

then for every time $T>0$, there exists a constant $C$ independent of $\Delta t$, $\Delta x$, and $\varepsilon$, such that the numerical solution obtained by scheme (6) satisfies the following error estimate

[TABLE]

This theorem is proved in the following sections: in section 4.1, we first derive the truncation error of the scheme, then in section 4.2, we apply the same analysis as for the stability result to prove the theorem.

4.1 Truncation error

Let $a^{n}_{i}$ and $\frac{1}{\varepsilon^{2}}b^{n}_{i}$ be the truncation errors of scheme (6), that is to say the reminders obtained by inserting the exact solution of (5) in relations (6):

[TABLE]

We can prove the following estimate of these truncation errors:

Lemma 4.1**.**

There exists a constant $\tilde{C}$ independent of $\Delta t$, $\Delta x$ and $\varepsilon$, such that

[TABLE]

for every $n$

This lemma is proved by using standard simple techniques (Taylor-Lagrange formula of different orders). But to simplify the paper, the proof—which is a bit long—is given in appendix A.

Now, let $\tilde{\rho}^{n}_{i}$ and $\tilde{g}^{n}_{i+\frac{1}{2}}$ be the convergence errors, that is to say the sequences defined by

[TABLE]

Then these sequences satisfy the ”perturbed” scheme, written in the following compact form:

[TABLE]

with the homogeneous initial data $\tilde{\rho}^{0}_{i}=0$ and $\tilde{g}^{0}_{i+\frac{1}{2}}=0$ for every $i$.

4.2 Analysis of the convergence error

In this section, we apply the same analysis as for the stability result to prove that $\tilde{\rho}^{n}$ and $\tilde{g}^{n}$ can be controlled by the truncation errors.

Step 1.

We multiply (38a) by $\tilde{\rho}^{n+1}_{i}$ and take the sum over $i$ to get

[TABLE]

Finally, we add up (39a) and $\varepsilon^{2}$ times (39b), we use (20) and lemma 3.4 to get

[TABLE]

Step 2.

Here, we can copy the steps 2 to 4 of the stability analysis (see section 3) for relation (40). Therefore, skipping the details, we just give the resulting energy estimate:

[TABLE]

Step 3.

We estimate the scalar products of the right-hand side of (41) by using two different Young inequalities:

[TABLE]

Moreover, by taking the positive part of the factor $(\frac{3\Delta t}{4}-\varepsilon\frac{\Delta x}{2})$