An asymptotic preserving scheme for Lévy-Fokker-Planck equation with fractional diffusion limit

Wuzhe Xu and Li Wang School of Mathematics, University of Minnesota, Twin cities, MN 55455. W.Xu ([email protected]), L.Wang ([email protected]).

Abstract

In this paper, we develop a numerical method for the Lévy-Fokker-Planck equation with the fractional diffusive scaling. There are two main challenges. One comes from a two-fold nonlocality, that is, the need to apply the fractional Laplacian operator to a power law decay distribution. The other arises from long-time/small mean-free-path scaling, which introduces stiffness to the equation. To resolve the first difficulty, we use a change of variable to convert the unbounded domain into a bounded one and then apply the Chebyshev polynomial based pseudo-spectral method. To treat the multiple scales, we propose an asymptotic preserving scheme based on a novel micro-macro decomposition that uses the structure of the test function in proving the fractional diffusion limit analytically. Finally, the efficiency and accuracy of our scheme are illustrated by a suite of numerical examples.

1 Introduction

We consider the Lévy-Fokker-Planck (LFP) equation

\begin{cases}{}\partial_{t}f+v\cdot\nabla_{x}f=\nabla_{v}\cdot(vf)-(-\Delta_{v})^{s}f:=\mathcal{L}^{s}(f)\,,\qquad s\in(0,1)\,,\\ f(0,x,v)=f_{in}(x,v)\,,\end{cases}

(1.1)

where $f(t,x,v):(0,\infty)\times\mathbb{R}^{d}\times\mathbb{R}^{d}\mapsto\mathbb{R}^{+}$ is the distribution function of a cloud of particles in plasma, which undergoes a free transport describing by the convection on the left hand side, and an interaction with the background, described by the LFP operator on the right. Here $(-\Delta_{v})^{s}$ is the fractional Lapacian operator that models the Lévy processes at the microscopic level. Among various equivalent ways of defining the fractional Laplace operator [14], we only mention the one that will be used throughout the paper:

(-\Delta_{v})^{s}f(v):=C_{s,d}\text{P.V.}\int_{\mathbb{R}^{d}}\frac{f(v)-f(w)}{|v-w|^{d+2s}}\mathrm{d}w,

(1.2)

where $P.V.$ denotes Cauchy principal value and $C_{s,d}=\frac{4^{s}\Gamma(d/2+s)}{\pi^{d/2}|\Gamma(-s)|}$ .

When $s=1$ , the right hand side of (1.1) reduces to the to classical Fokker-Planck operator that models the Brownian motion of the microscopic particles. In contrast, $s\in(0,1)$ allows particles to make long jumps at the microscopic scale, and hence leads to the nonlocal effect at the mesoscopic scale. Consequently, as opposed to the Gaussian Maxwellian to the Fokker-Planck operator, the Lévy-Fokker-Planck operator admits an equilibrium that has a fat tail. More precisely, there is a unique normalized distribution $\mathcal{M}(v)$ , such that [9]:

\mathcal{L}^{s}(\mathcal{M})=0,\quad\int_{\mathbb{R}^{d}}\mathcal{M}(v)\mathrm{d}v=1,\quad\mathcal{M}(v)\sim\frac{C}{|v|^{d+2s}}~{}\text{as}~{}|v|\rightarrow\infty\,.

(1.3)

Moreover, the convergence toward the equilibrium is shown to be exponential in the sense of relative entropy. In particular, let $\Phi$ be a convex smooth function, and define the relative entropy to the equilibrium as

H^{\Phi}(f|\mathcal{M}):=\int\Phi(f)\mathcal{M}\mathrm{d}v-\Phi\left(\int f\mathcal{M}\mathrm{d}v\right)\,,

then from the Theorem 2 in [9], one has

H^{\Phi}(f|\mathcal{M})\left(\frac{f(t)}{\mathcal{M}}\right)\leq e^{-\frac{t}{C}}H^{\Phi}(f_{in}|\mathcal{M})\left(\frac{f(t)}{\mathcal{M}}\right),\quad t\geq 0\,,

for some constant $C$ .

In the long time and small mean free path scaling, that is, rescaling the time and space as

x\mapsto\varepsilon x,\quad t\mapsto\varepsilon^{2s}t,

(1.4)

equation (1.1) can be rewritten into the following form:

\begin{cases}{}\varepsilon^{2s}\partial_{t}f^{\varepsilon}+\varepsilon v\cdot\nabla_{x}f^{\varepsilon}=\mathcal{L}^{s}(f^{\varepsilon}),\\ f^{\varepsilon}(0,x,v)=f_{in}(x,v)\,.\end{cases}

(1.5)

Sending $\varepsilon$ to zero will give rise to a fractional diffusion equation for the density of particles. More precisely, we cite the following theorem from [6].

Theorem 1.

Assume that $f_{0}\in L^{2}(\mathbb{R}^{N},\mathcal{M}(v)^{-1}\mathrm{d}v\mathrm{d}x)$ , where $\mathcal{M}(v)$ is the unique normalized equilibrium distribution that satisfies (1.3). Then, up to a subsequence, the solution $f^{\varepsilon}$ of (1.5) converges weakly in $L^{\infty}(0,T;L^{2}(\mathbb{R}^{2d},\mathcal{M}(v)^{-1}\mathrm{d}v\mathrm{d}x))$ to $\rho(t,x)\mathcal{M}(v)$ as $\varepsilon\rightarrow 0$ , where $\rho(t,x)$ solves

\begin{cases}{}\partial_{t}\rho+(-\Delta_{x})^{s}\rho=0\,,\\ \rho(0,x)=\rho_{in}(x):=\int_{\mathbb{R}^{d}}f_{in}(x,v)\mathrm{d}v\,.\end{cases}

(1.6)

In the classical case (i.e., $s$ =1) when $\mathcal{M}$ is a fast decaying function such as Gaussian, one rescales $t$ as $t\mapsto\varepsilon^{2}t$ and the resulting macroscopic equation is the diffusion equation [20]:

\partial_{t}\rho+\nabla_{x}\cdot(D\nabla_{x}\rho)=0,

where $D$ is the diffusion matrix

D=\int v\otimes v\mathcal{M}\mathrm{d}v\,.

Clearly the fat tail equilibrium (1.3) renders the above integral unbounded and therefore invalids the classical diffusion limit. Conversely, the anomalous scaling (1.4) is necessary. Similar scaling has also been investigated in the framework of linear Boltzmann equation, see [19, 1] for a reference.

Numerically computing (1.5) has two major challenges. One comes from necessity to apply the fractional Laplacian operator to a slow decay function, in which case one needs to consider infinite computational domain, since any truncation would lose the important information carried by the tail and leads to erroneous result. To this end, two kinds of numerical methods have been developed to approximate the fractional Laplacian operator. One hinges upon the integral form of $(-\Delta_{v})^{s}$ in (1.2) and splits the infinite domain into a computable body part and a compensate tail part that can be integrated exactly, see especially [11, 12]. This method heavily relies on analytic expression of the tail, which is not known for our case except when the solution has reached equilibrium. But since our goal is to simulate the dynamics, this method cannot be adopted without a major modification. The other uses the spectral method with nonlocal basis. For instance, the Hermite spectral method is developed in [18] which takes advantage of the fact that Hermite polynomials are invariant under Fourier transform and therefore can be efficiently computed when taking the fractional Laplacian. In [21], mapped Chebyshev polynomials are used as basis, then the fractional Laplacian is calculated via the celebrated Dunford-Taylor formula, and therefore is very efficient in higher dimensions. Another approach, proposed in [5], also starts from the Chebyshev polynomial basis, but it then uses a change of variable and even extension of the function, and therefore boils down the problem to computing the Fractional Laplacian of the Fourier basis on a finite domain, which can be approximated numerically with high accuracy. All these spectral methods were developed to address the nonlocality of the fractional Lapalcian operator, but when they apply to a slow decay function, additional nonlocality is introduced and a large number of basis are expected in order to meet a certain accuracy. See also [16, 3] for a review on numerical issues related to fractional diffusion. For our problem, we find that the approach in [5] gives the best result within the computational budget when some tuning parameters are chosen appropriately.

Another challenge arises from the diffusive scaling, which introduces stiffness to the system. Our goal is to develop a numerical solver with uniform performance across different regimes, i.e., $\varepsilon$ varies in magnitude by several orders. In particular, the scheme for (1.5) is expected to reduce to the solver for (1.6) automatically with unresolved mesh. This is the so-called asymptotic preserving (AP) scheme. There has been an extensive study on the AP scheme for kinetic equations with various scalings, see [13, 10] for a review. When it comes to anomalous diffusive scaling, we cite specially [23, 7, 8, 24], which all deal with the linear Boltzmann type equation. These works and the current paper share the same equilibrium (1.3) and diffusion limit (1.6), but the different nature between Fokker Planck type and Boltzmann type operator lead to very difference convergence mechanism and therefore hinders the application of the methods developed there. In fact, when the collision is of the Boltzmann type, a strong convergence toward the anomalous diffusion is obtained [19, 1], which plays a significant role in designing the AP method. On the contrary, with the Lévy-Fokker-Planck operator in our case, only a weak convergence is available, which gives very limited knowledge on how the scheme can be constructed. Nevertheless, the special choice of the test function that aids the proof of Theorem 1 in [6] inspires our macro-micro decomposition, which sets the base of our AP scheme.

The rest of the paper is organized as follows. In the next section, we detail the computation of the collision operator $\mathcal{L}^{s}$ and combine it with backward Euler scheme to solve the spatially homogeneous case. Section 3 is devoted to the design of AP scheme, along with a rigorous proof of the AP property and a detailed guide in implementation. In section 4, extensive numerical examples are given to test the performance of our AP scheme. Finally the paper is concluded in section 5.

2 Computation of the collision operator $\mathcal{L}^{s}$

Aside from multiple scales that appear in equation (1.5), the collision operator $\mathcal{L}^{s}$ itself poses severe computational difficulties, which is attributed to the non-locality of the operator $(-\Delta_{v})^{s}$ . As already mentioned in [11, 12], if $f$ is compactly supported in $v$ , then the computation of $(-\Delta_{v})^{s}f$ can be fulfilled by the Fourier transform. However, in our case, the interplay between the two operators in $\mathcal{L}^{s}$ leads to an equilibrium that has only a power law decay, see (1.3). As a result, a more sophisticated treatment is needed for fractional Laplacian, as will be detailed in the following section. Here for notation simplicity, we assume $f$ only depends on $v$ throughout this section.

2.1 Change of variable

As mentioned above, one of the major difficulties of fractional Laplacian operator $(-\Delta_{v})^{s}$ is its non-locality, especially when it applies to a slow decaying function like $\mathcal{M}$ in (1.3). In order to treat the fat tail distribution on an unbounded domain, two approaches can be taken: one is to truncate the computation domain and introduce suitable boundary conditions; the other is to use a change of variable that maps the infinite domain into a finite one and then use a spectral method. In this paper, we take the latter approach, which is also termed as the rational spectral method. Below we briefly introduce the rational Chebyshev spectral method, more details can be found in [5].

Consider an algebraic mapping that maps the unbounded domain $(-\infty,\infty)$ into $[-1,1]$ , i.e. ,

\xi=\frac{v}{\sqrt{L_{v}^{2}+v^{2}}}\in(-1,1)\Longleftrightarrow v=\frac{L_{v}\xi}{\sqrt{1-\xi^{2}}}\in(-\infty,\infty)\,,

where $L_{v}$ is a scaling parameter that is chosen for the sake of accuracy and mass conservation. In principle, it can be chosen adaptive as mentioned in [17], see also Remark 2 for more discussion. Take the Chebyshev polynomials of the first kind as the basis on $[-1,1]$ ,

T_{k}(\xi)=\cos(k\arccos(\xi))\,,\qquad\xi\in[-1,1]\,,

(2.1)

then

TB_{k}(v)=T_{k}(\frac{v}{\sqrt{L_{v}^{2}+v^{2}}})\,,\qquad v\in\mathbb{R}\,,

is the so-called Chebyshev rational polynomials on infinite domain. It has been pointed out in [4] that Chebyshev rational polynomials are appropriate for approximating the algebraically decay function, and has also been used in [21] to compute the fractional diffusion operator.

We now concentrate on the finite domain $[-1,1]$ and employ a further change of variable. Let $q=\arccos(\xi)\in[0,\pi]$ , then

v=\frac{L_{v}\xi}{\sqrt{1-\xi^{2}}}=L_{v}\cot(q),

(2.2)

and $T_{k}(\xi)$ in (2.1) reduces to

T_{k}(\xi)=\cos(kq).

Therefore, expanding $f(q)$ in terms of Chebyshev functions is equivalent to a cosine expansion.

In the new variable $q\in[0,\pi]$ , the fractional Laplacian can be rewritten as [5]

(-\Delta_{q})^{s}f(q)=\left\{\begin{array}[]{cc}-\frac{1}{L_{v}\pi}\int_{0}^{\pi}\frac{f^{\prime}(p)}{\cot(q)-\cot(p)}\mathrm{d}p,&s=\frac{1}{2}\,,\\ \frac{C_{s,d}}{2L_{v}^{2s}s(1-2s)}\int_{0}^{\pi}\frac{\sin^{2}(p)f^{\prime\prime}(p)+2\sin(p)\cos(p)f^{\prime}(p)}{|\cot(q)-\cot(p)|^{2s-1}}\mathrm{d}p,&s\neq\frac{1}{2}\,.\end{array}\right.

And one can reformulate (1.1) into

\begin{cases}{}\partial_{t}f+L_{v}\cot(q)\partial_{x}f=f-\cos(q)\sin(q)\partial_{q}f(q)-(-\Delta_{q})^{s}f:=\mathcal{L}^{s}_{q}(f)\,,\\ f(0,x,q)=f_{in}(x,q)\,,\end{cases}

(2.3)

where $f$ now depends on $t$ , $x$ and $q\in[0,\pi]$ . In the rest of this section, we will use (2.3) as our target equation, and discretize $q$ in the following way:

q_{j}=\frac{\pi(2j+1)}{2N_{v}},\qquad 0\leq j\leq N_{v}-1\,,

(2.4)

with $\Delta q=\frac{\pi}{N_{v}}$ .

2.2 Computation of $(-\Delta_{q})^{s}f$

To further ease the computation of $(-\Delta_{q})^{s}f$ , we conduct the even extension of $f$ at $\pi$ , i.e.,

\tilde{f}(q)=\left\{\begin{array}[]{cc}f(q),&q\in[0,\pi]\,,\\ f(2\pi-q),&q\in[\pi,2\pi]\,.\end{array}\right.

This way, according to the relation

\int_{0}^{2\pi}\tilde{f}(q)\cos(kq)\mathrm{d}q=\int_{0}^{2\pi}\tilde{f}(q)e^{-ikq}\mathrm{d}q\,,

one can compute the coefficients of cosine expansion of $\tilde{f}$ via the Fast Fourier Transform. More specifically, denote

\tilde{\mathsf{f}}=[\tilde{f}(q_{0}),\tilde{f}(q_{1}),\cdots,\tilde{f}(q_{N_{v}-1}),f(q_{2N_{v}-1}),\tilde{f}(q_{2N_{v}-2}),\cdots,\tilde{f}(q_{N_{v}})]^{T},

(2.5)

then its discrete Fourier transform takes the form

f(q_{j})=\sum_{k=-N_{v}}^{N_{v}-1}\hat{\tilde{f}}_{k}e^{ikq_{j}}\,,\qquad j=0,\cdots,2N_{v}-1\,,

where $\hat{\tilde{f}}_{k}=\frac{1}{2N_{v}}\sum_{j=-N_{v}}^{N_{v}-1}\tilde{f}(q_{j})e^{-ikq_{j}}$ , and hence

(-\Delta_{q})^{s}f(q_{j})=\sum_{k=-N_{v}}^{N_{v}-1}\hat{\tilde{f}}_{k}(-\Delta_{q})^{s}e^{ikq_{j}}\,,\qquad j=0,\cdots,2N_{v}-1\,.

Then the question boils down to calculating $(-\Delta_{q})^{s}e^{ikq_{j}}$ , and we directly cite the result from [5].

Theorem 2.

Let $s\in(0,0.5)\cup(0.5,1)$ , then

(-\Delta_{q})^{s}\left(e^{imq}\right)=\left\{\begin{array}[]{cc}\frac{c_{s,1}|\sin(q)|^{2s-1}}{8L^{2s}\tan\left(\pi s\right)}\sum_{l=-\infty}^{\infty}e^{i2lq}\left((1-2s)m^{2}-4ml\right)\\ \times\frac{\Gamma\left(\frac{-1+2s}{2}+|l|\right)\Gamma\left(\frac{-1-2s}{2}+\left|\frac{m}{2}-l\right|\right)}{\Gamma\left(\frac{3-2s}{2}+|l|\right)\Gamma\left(\frac{3+2s}{2}+\left|\frac{m}{2}-l\right|\right)},&m\text{ even}\,,\\ i\frac{c_{s,1}|\sin(q)|^{2s-1}}{8L^{2s}}\sum_{l=-\infty}^{\infty}e^{i2lq}\left((1-2s)m^{2}-4ml\right)\\ \times\operatorname{sgn}\left(\frac{m}{2}-l\right)\frac{\Gamma\left(\frac{-1+2s}{2}+|l|\right)\Gamma\left(\frac{-1-2s}{2}+\left|\frac{m}{2}-l\right|\right)}{\Gamma\left(\frac{3-2s}{2}+|l|\right)\Gamma\left(\frac{3+2s}{2}+\left|\frac{m}{2}-l\right|\right)},&m\text{ odd }\,.\end{array}\right.

(2.6)

Moreover, when $s=0.5$ ,

(-\Delta_{q})^{0.5}\left(e^{imq}\right)=\left\{\begin{array}[]{ll}\frac{|m|\sin^{2}(q)}{L}e^{imq},&m\text{ even}\,,\\ \frac{im}{L\pi}\left(\frac{-2}{m^{2}-4}-\sum_{l=-\infty}^{\infty}\frac{4\operatorname{sgn}(l)e^{i2lq}}{(m-2l)\left((m-2l)^{2}-4\right)}\right),&m\text{ odd }\,.\end{array}\right.

To implement it numerically, one truncates $l$ by setting $l=l_{1}N_{v}+l_{2}$ , where $l_{2}\in\{-N_{v}/2,\ldots,N_{v}/2-1\}$ , and $l_{1}\in\left\{-l_{lim},\ldots,l_{lim}\right\}$ . Then (2.6) is approximated as

(-\Delta_{q})^{s}\left(e^{imq_{j}}\right)\approx\left\{\begin{array}[]{l}\frac{c_{s,1}\left|\sin\left(q_{j}\right)\right|^{2s-1}}{8L^{2s}\tan\left(\frac{\pi 2s}{2}\right)}\sum_{l_{2}=-N_{v}/2}^{N_{v}/2-1}\left[\sum_{l_{1}=-l_{lim}}^{l_{lim}}a_{m,l_{1},l_{2}}\right]e^{i2l_{2}q_{j}},\quad m\text{ even }\,,\\ i\frac{c_{s,1}\left|\sin\left(q_{j}\right)\right|^{2s-1}}{8L^{2s}}\sum_{l_{2}=-N_{v}/2}^{N_{v}/2-1}\left[\sum_{l_{1}=-l_{lim}}^{l_{lim}}a_{m,l_{1},l_{2}}\right]e^{i2l_{2}q_{j}},\quad m\text{ odd }\,,\end{array}\right.

(2.7)

where

a_{m,l_{1},l_{2}}=\left\{\begin{array}[]{ll}(-1)^{l_{1}}\left((1-2s)m^{2}-4m\left(l_{1}N_{v}+l_{2}\right)\right)\\ \times\frac{\Gamma\left(\frac{-1+2s}{2}+\left|l_{1}N_{v}+l_{2}\right|\right)\Gamma\left(\frac{-1-2s}{2}+\left|\frac{m}{2}-l_{1}N_{v}-l_{2}\right|\right)}{\Gamma\left(\frac{3-2s}{2}+\left|l_{1}N_{v}+l_{2}\right|\right)\Gamma\left(\frac{3+2s}{2}+\left|\frac{m}{2}-l_{1}N_{v}-l_{2}\right|\right)},&m\text{ even }\,,\\ (-1)^{l_{1}}\left((1-2s)m^{2}-4m\left(l_{1}N_{v}+l_{2}\right)\right)\operatorname{sgn}\left(\frac{m}{2}-l_{1}N_{v}-l_{2}\right)\\ \times\frac{\Gamma\left(\frac{-1+2s}{2}+\left|l_{1}N_{v}+l_{2}\right|\right)\Gamma\left(\frac{-1-2s}{2}+\left|\frac{m}{2}-l_{1}N_{v}-l_{2}\right|\right)}{\Gamma\left(\frac{3-2s}{2}+\left|l_{1}N_{v}+l_{2}\right|\right)\Gamma\left(\frac{3+2s}{2}+\left|\frac{m}{2}-l_{1}N_{v}-l_{2}\right|\right)},&m\text{ odd }\,.\end{array}\right.

Note that $l_{lim}$ here is an adjustable parameter, and it is obvious that larger $l_{lim}$ gives better approximation. For all our numerical examples in this paper, we use $l_{lim}=300$ .

Now let the Matrix $\mathsf{M}\in\mathbb{C}^{2N_{v}\times 2N_{v}}$ be

\mathsf{M}_{m,n}=\left\{\begin{array}[]{cc}(-\Delta_{q})^{s}e^{i(m-N_{v})q_{n-1}}&n\leq N_{v}\,,\\ (-\Delta_{q})^{s}e^{i(m-N_{v})q_{2N_{v}-n}}&N_{v}+1\leq n\leq 2N_{v}\,,\end{array}\right.

then we have

(-\Delta_{q})^{s}\tilde{\mathsf{f}}=(\mathsf{M}\times\mathsf{F})~{}\tilde{\mathsf{f}}\,,

where $\tilde{\mathsf{f}}$ is defined in (2.5), and $\mathsf{F}$ denotes the $2N_{v}$ -periodic discrete Fourier transform. Confining the above calculation of $\tilde{\mathsf{f}}$ to $\mathsf{f}$ , i.e.,

\mathsf{f}=(f(q_{0}),f(q_{1}),\cdots,f(q_{N_{v}-1}))^{T}.

(2.8)

we can write down the matrix representation of $(-\Delta_{q})^{s}\mathsf{f}$ as follows:

(-\Delta_{q})^{s}\mathsf{f}=\mathsf{L}_{s}\mathsf{f}\,,

(2.9)

where $\mathsf{L}_{s}=(\mathsf{M}\times\mathsf{F})(1:N_{v},1:N_{v})+(\mathsf{M}\times\mathsf{F})(1:N_{v},N_{v}+1:2N_{v})$ .

As already mentioned in the introduction, the way we treat the fractional Laplacian is not unique. We just choose the one that performs the best in our case in terms of accuracy and efficiency.

2.3 Spatially homogeneous case

In this section, we detail the computation of the spatially homogeneous case of (2.3):

\partial_{t}f=\mathcal{L}^{s}_{q}(f):=f-\cos(q)\sin(q)\partial_{q}f-(-\Delta_{q})^{s}f\,.

Here the fractional Laplacian term is treated via the aforementioned pseudospectral method, and $\partial_{q}f$ is discretized using the Fourier spectral method. Still using the vector form of the discrete $f$ defined in (2.8), we have the discretization of $\mathcal{L}^{s}_{q}f$ as

\mathsf{P}^{s}\mathsf{f}:=(\mathsf{C}\cdot\mathsf{F}^{-1}\mathsf{K}\mathsf{F}+\mathsf{I}-\mathsf{L}_{s})\mathsf{f}\,,

(2.10)

where

\mathsf{K}=i\begin{pmatrix}-N_{v}&&\\ &\ddots&\\ &&N_{v}-1\end{pmatrix},\quad\mathsf{C}=\begin{pmatrix}-\cos(q_{0})\sin(q_{0})&&\\ &\ddots&\\ &&-\cos(q_{N_{v}-1})\sin(q_{N_{v}-1})\end{pmatrix}\,.

Note specifically that even though the computation of $\mathsf{L}_{s}$ can be expensive, it only needs to be computed once.

For time discretization, we choose the backward Euler scheme, namely,

\frac{1}{\Delta t}(\mathsf{f}^{n+1}-\mathsf{f}^{n})=\mathsf{P}^{s}\mathsf{f}^{n+1}\,,

(2.11)

as it will be used in the spatially in-homogeneous case in the next section. In practice, the ODE system is solved by Matlab builtin solver ODE45 or ODE15s.

Remark 1 (Positivity).

The scheme (2.11) is not positivity preserving, but it will be so after a slight modification. In fact, the positivity is lost when $q_{j}$ close to the boundary, $0$ or $\pi$ . This is because the function value $f$ near the boundary are already small, then any small numerical error would render it negative. To resolve this issue, our idea is to use the tail information to reassign the value of $\mathsf{f}$ at the boundary. More precisely, since the equilibrium is proportional to $(1+L_{v}\cot(q))^{-(1+2s)}$ at its tail, then if there is an index $l$ close to the right boundary (i.e., $q=\pi$ ) such that $\mathsf{f}_{l}<0$ for the first time, namely $\mathsf{f}_{j}>0$ for $j<l$ , then we set $\mathsf{f}_{j}=\frac{(1+L_{v}\cot(q(j-1)))^{1+2s}}{(1+L_{v}\cot(qj))^{1+2s}}$ for $j\geq l$ .

Remark 2 (Choice of $L_{v}$ and $N_{v}$ ).

In the scheme $\eqref{eqn: homo_scheme_ori}$ , $N_{v}$ and $L_{v}$ should be chosen according to $s$ for the sake of mass conservation. Since the tail of equilibrium distribution relates to $s$ via $\mathcal{M}\sim\frac{1}{|v|^{1+2s}}$ as $|v|\rightarrow\infty$ , larger $N_{v}$ (with fix $L_{v}$ ) is required for smaller $s$ to capture the tail information and conserve the total mass. One should also choose the parameter $L_{v}$ properly. On one hand, for fixed $N_{v}$ , $L_{v}$ should not be too large, otherwise accuracy is lost in approximating the body part of the distribution function. On the other hand, if $L_{v}$ is too small, the method lose the capability to capture the tail information. Therefore, in principle, the smaller the $s$ is, the larger number of mode $N_{v}$ and $L_{v}$ are needed. For instance, when $s\geq 0.5$ , $N_{v}=64$ and $L_{v}=3$ is enough; whereas for $s=0.4$ and if $L_{v}=3$ , $N_{v}=128$ is required to ensure the mass is conserved up to error $\mathcal{O}(10^{-3})$ , see Table 1.

3 Asymptotic preserving scheme

We introduce an asymptotic preserving scheme for the spatially inhomogeneous system and its implementation in this section. The main difficulty of capturing the anomalous diffusion limit is to acquire operator $(-\Delta_{x})^{s}$ when only $(-\Delta_{v})^{s}$ appears in original system. In the proof of Theorem 1, the authors use a test function in the form of $\phi(x+\varepsilon v)$ to get $(-\Delta_{x})^{s}$ from $(-\Delta_{v})^{s}$ weakly via integration by parts, which indicates a kind of symmetry between $x$ and $v$ . Inspired by this symmetry, we propose a scheme based on a novel micro-macro decomposition by requiring the macro part to respect such symmetry.

Let us first introduce some notation. Denote $\Omega_{x}=[-L_{x},L_{x}]$ as our spatial domain and partition it into $N_{x}$ uniform grids, i.e. $x_{i}=-L_{x}+(i+\frac{1}{2})\Delta x$ , where $\Delta x=\frac{2L_{x}}{N_{x}}$ . Periodic boundary condition in the spatial domain will be used. For velocity domain, we will work with the variable $q$ defined in (2.2) and use the same discretization as in (2.4). Then $f(t,x,q)=f(t,x,v(q))$ , and we denote numerical approximation of $f(t_{n},x_{i},q_{j})$ as $f^{n}_{i,j}$ , where $t_{n}=n\Delta t$ , $0\leq j\leq N_{v}-1$ , $0\leq i\leq N_{x}-1$ .

To start, we introduce the following lemmas.

Lemma 1.

$(-\Delta_{v})^{s}(fg)=g(-\Delta_{v})^{s}f+f(-\Delta_{v})^{s}g+I(f,g)$ , where

I(f,g)=C_{1,s}\int_{\mathbb{R}^{d}}\frac{(f(v)-f(w))(g(w)-g(v))}{|v-w|^{d+2s}}\mathrm{d}w\,.

Proof.

$\displaystyle(-\Delta_{v})^{s}(fg)$	$\displaystyle=$	$\displaystyle C_{1,s}\int_{\mathbb{R}}\frac{f(v)g(v)-f(w)g(w)}{\|v-w\|^{d+2s}}\mathrm{d}w\,$
	$\displaystyle=$	$\displaystyle C_{1,s}\int_{\mathbb{R}}\frac{f(v)g(v)-f(v)g(w)}{\|v-w\|^{d+2s}}\mathrm{d}w+C_{1,s}\int_{\mathbb{R}}\frac{f(v)g(v)-f(w)g(v)}{\|v-w\|^{d+2s}}\mathrm{d}w\,$
		$\displaystyle+C_{1,s}\int_{\mathbb{R}}\frac{(f(v)-f(w))(g(w)-g(v))}{\|v-w\|^{d+2s}}\mathrm{d}w\,$
	$\displaystyle=$	$\displaystyle g(-\Delta_{v})^{s}f+f(-\Delta_{v})^{s}g+I(f,g).$

∎

Lemma 2.

Suppose $h(x,v)\in L^{1}_{v}(\mathbb{R})W^{2,1}_{x}(\mathbb{R})\cap L_{v}^{1}(\mathbb{R})C^{2}_{x}(\mathbb{R})$ , then $\left\langle(-\Delta_{x})^{s}h(x,v)\right\rangle=(-\Delta_{x})^{s}\left\langle h\right\rangle$ , where $\left\langle f\right\rangle:=\int_{\mathbb{R}}f(x,v)\mathrm{d}v$ .

Proof.

First, let us rewrite fractional Laplacian in following finite difference form,

(-\Delta_{x})^{s}h(x,v)=C_{1,s}\int_{0}^{\infty}(2h(x,v)-h(x-y,v)-h(x+y,v))\nu(y)\mathrm{d}y\,,

where $\nu(y)=|y|^{-(1+2s)}$ , which is a direct consequence of (1.2). Denote

	$\displaystyle I(x,v)$	$\displaystyle=\int_{0}^{\infty}\|(2h(x,v)-h(x-y,v)-h(x+y,v))\nu(y)\|dy$
		$\displaystyle=\int_{0}^{\delta}\|(2h(x,v)-h(x-y,v)-h(x+y,v))\nu(y)\|dy$
		$\displaystyle\qquad\qquad+\int_{\delta}^{\infty}\|(2h(x,v)-h(x-y,v)-h(x+y,v))\nu(y)\|dy$
		$\displaystyle:=I_{1}(x,v)+I_{2}(x,v)\,,$

for $\delta<1$ . By Taylor expansion we have

2h(x,v)-h(x+y,v)-h(x-y,v)=-\partial_{x}^{2}h(\xi,v)y^{2}\,,

where $\xi\in(x-\delta,x+\delta)$ . Since $s\in(0,1)$ , we have $\int_{0}^{\delta}y^{n}\nu(y)\mathrm{d}y<\infty$ . The assumption that $h(x,v)\in L^{1}_{v}(\mathbb{R})W^{2,1}_{x}(\mathbb{R})$ then leads to $\int_{\mathbb{R}}|\partial_{x}^{2}h(\xi,v)|\int_{0}^{\delta}y^{2}\nu(y)\mathrm{d}y\mathrm{d}v<\infty\,.$ On the other hand, it’s obvious that

\int_{\mathbb{R}}\int_{\delta}^{\infty}|(2h(x,v)-h(x-y,v)-h(x+y,v))\nu(y)|\mathrm{d}y\mathrm{d}v<\infty.

Then we conclude the result by the Fubini’s theorem.

∎

3.1 A novel micro-macro decomposition

A typical approach in designing an asymptotic preserving method is via a micro-macro decomposition. That is, write

f(t,x,v)=\rho(t,x)\mathcal{M}(v)+g(t,x,v)\,,

(3.1)

where $\rho(t,x)$ is the macroscopic density, and $g$ is viewed as the microscopic fluctuation. See for instance [15, 23, 24, 7]. However, directly applying this decomposition to our case fails as it is not easy to obtain $(-\Delta_{x})^{s}$ when only $(-\Delta_{v})^{s}$ appears in $\eqref{eqn:111}$ . Inspired by the proof of Theorem 1 in [6], we see that the operator $(-\Delta_{x})^{s}$ and operator $\varepsilon^{2s}(-\Delta_{v})^{s}$ are related by considering a test function in variable $(x+\varepsilon v)$ . Therefore, we propose the following novel micro-macro decomposition of distribution $f$ :

f(t,x,v)=\eta(t,x,v)\mathcal{M}(v)+g(t,x,v)\,,

(3.2)

where $\eta(t,x,v)$ takes the form

\eta(t,x,v)=h(t,x+\varepsilon v)\,,

(3.3)

for some function $h(t,x)$ , and $\mathcal{M}$ satisfies

\partial_{v}(v\mathcal{M})-(-\Delta_{v})^{s}\mathcal{M}=0,\qquad\int_{\mathbb{R}}\mathcal{M}(v)\mathrm{d}v=1\,.

(3.4)

Here unlike the classical micro-macro decomposition (3.1), we allow $\eta$ to depend on $v$ , but intrinsically $\eta$ lives on a lower dimensional manifold than $f$ does, as required from (3.3). As a result, $\eta$ satisfies

\varepsilon\partial_{x}\eta=\partial_{v}\eta,\qquad(-\Delta_{v})^{s}\eta=\varepsilon^{2s}(-\Delta_{x})^{s}\eta\,.

(3.5)

Plug $\eqref{eqn: decomp}$ into $\eqref{eqn:111}$ , we get:

\varepsilon^{2s}\partial_{t}(\eta\mathcal{M}+g)+\varepsilon v\partial_{x}(\eta\mathcal{M}+g)=\partial_{v}(v(\eta\mathcal{M}+g))-(-\Delta_{v})^{s}(\eta\mathcal{M}+g)\,.

(3.6)

Using (3.4),(3.5) and Lemma 1, the above equation simplifies to

\varepsilon^{2s}\partial_{t}(\eta\mathcal{M}+g)+\varepsilon v\partial_{x}g=\mathcal{L}^{s}(g)-\varepsilon^{2s}(-\Delta_{x})^{s}\eta\mathcal{M}-I(\eta,\mathcal{M}).

(3.7)

To solve (3.7), we split it into the following system

\begin{cases}{}\varepsilon^{2s}\partial_{t}g+\varepsilon v\partial_{x}g=\mathcal{L}^{s}(g)-I(\eta,\mathcal{M})\,,\\ \partial_{t}\eta=-(-\Delta_{x})^{s}\eta\,,\end{cases}

(3.8)

equipped with the initial condition

\eta_{in}(x,v)=\rho_{in}(x+\varepsilon v)\,\qquad g_{in}(x,v)=f_{in}(x,v)-\eta_{in}(x,v)\mathcal{M}(v)\,.

(3.9)

Upon solving (3.8), one can recover $f$ using (3.2).

Note that the reduction of (3.7) from (3.6) is possible only when $\eta$ has the form (3.3). Therefore, it is important to make sure that such property is preserved along the dynamics of (3.8). To this end, we have the following lemma.

Lemma 3.

Let $\eta(t,x,v)$ solves

\partial_{t}\eta=-(-\Delta_{x})^{s}\eta\,,\qquad\eta(0,x,v)=\rho_{in}(x+\varepsilon v)

(3.10)

then there exists a function h(t,x) such that $\eta(t,x,v)=h(t,x+\varepsilon v)$ .

Proof.

Let $h$ satisfies

\partial_{t}h=-(-\Delta_{x})^{s}h\,,\qquad h(0,x,v)=h_{in}(x)\,.

(3.11)

Then we claim that $\eta(t,x,v)=h(t,x+\varepsilon v)$ is the solution to (3.10). Indeed, denote $y=x+\varepsilon v$ for any fixed $v$ , we have

	$\displaystyle\partial_{t}h(t,x+\varepsilon v)=\partial_{t}h(t,y)=-(-\Delta_{y})^{s}h$	$\displaystyle=-\int_{\mathbb{R}}\frac{h(t,y)-h(t,y^{\prime})}{\|y-y^{\prime}\|^{2s+1}}\mathrm{d}y^{\prime}$
		$\displaystyle=-\int_{\mathbb{R}}\frac{h(t,x+\varepsilon v)-h(t,x^{\prime}+\varepsilon v)}{\|(x+\varepsilon v)-(x^{\prime}+\varepsilon v)\|^{2s+1}}\mathrm{d}(x^{\prime}+\varepsilon v)$
		$\displaystyle=-(-\Delta_{x})^{s}h(t,x+\varepsilon v)\,.$

∎

Remark 3.

From the above lemma, one sees that in order to solve the last equation in (3.8), one can solve the low dimensional problem (3.11), and then obtain $\eta$ by shifting $h$ , i.e., $\eta(t,x,v)=h(t,x+\varepsilon v)$ .

Next, we show that the system (3.8) is energy stable.

Proposition 1.

If $(\eta,g)$ solves (3.8) with initial data (3.9), then $f=\eta\mathcal{M}+g$ solves (1.5). Both system has the energy dissipation property. That is, define the total energy

E_{f}^{2}=\int_{\mathbb{R}}\int_{\mathbb{R}}\frac{f^{2}}{\mathcal{M}}\mathrm{d}v\mathrm{d}x=\int_{\mathbb{R}}\int_{\mathbb{R}}\frac{(\eta\mathcal{M}+g)^{2}}{\mathcal{M}}\mathrm{d}v\mathrm{d}x\,,

(3.12)

then $\frac{dE_{f}}{dt}\leq 0$ .

Proof.

It is easy to show that if $(\eta,g)$ solves (3.8) with initial data (3.9), then by directly adding the two equation, $f=\eta\mathcal{M}+g$ solves (1.5). The energy dissipation of the original system (1.5) follows from Proposition 2.1 in [6]. One just needs to check the same property for the split system (3.8). To this end, multiply the first equation in (3.8) by $\frac{g}{\mathcal{M}}$ , the second equation by $\frac{\eta}{\mathcal{M}}$ , integrate in $x$ and $v$ , and add them together, we get (here we omit the integration domain, which is $\mathbb{R}$ for both $x$ and $v$ ):

	$\displaystyle\varepsilon^{2a}\partial_{t}\int\int\left[\frac{1}{2}\frac{g^{2}}{\mathcal{M}}+\frac{1}{2}\eta^{2}\mathcal{M}+g\eta\right]\mathrm{d}x\mathrm{d}v+\varepsilon\int\int v\partial_{x}g\frac{g}{\mathcal{M}}\mathrm{d}x\mathrm{d}v$
	$\displaystyle=\int\int\mathcal{L}^{s}(g)\left(\frac{g}{\mathcal{M}}+\eta\right)-I(\eta,\mathcal{M})\left(\frac{g}{\mathcal{M}}+\eta\right)-(\eta\mathcal{M}+g)(-\Delta_{v})^{s}\eta-\varepsilon v\partial_{x}g\eta\mathrm{d}x\mathrm{d}v$
	$\displaystyle=\int\int(\mathcal{L}^{s}(g)+\mathcal{L}^{s}(\eta\mathcal{M}))\left(\frac{g}{\mathcal{M}}+\eta\right)-\partial_{v}(v\eta\mathcal{M})\left(\frac{g}{\mathcal{M}}+\eta\right)+\eta(-\Delta_{v})^{s}\mathcal{M}\left(\frac{g}{\mathcal{M}}+\eta\right)-\varepsilon v\partial_{x}g\eta\mathrm{d}x\mathrm{d}v$
	$\displaystyle=\int\int(\mathcal{L}^{s}(g)+\mathcal{L}^{s}(\eta\mathcal{M}))\left(\frac{g}{\mathcal{M}}+\eta\right)-v\mathcal{M}\partial_{v}\eta\left(\frac{g}{\mathcal{M}}+\eta\right)-\varepsilon v\partial_{x}g\eta\mathrm{d}x\mathrm{d}v\,,$

where the second equality uses Lemma 1 and the third equality uses the fact that $\mathcal{M}$ satisfies $\mathcal{L}^{s}(\mathcal{M})=0$ . Since $\partial_{v}\eta=\varepsilon\partial_{x}\eta$ and $\int\int v\partial_{x}\eta g+\partial_{x}g\eta\mathrm{d}x\mathrm{d}v=\int\int v\mathcal{M}\partial_{x}\eta\eta\mathrm{d}x\mathrm{d}v=\int\int v\partial_{x}g\frac{g}{\mathcal{M}}\mathrm{d}x\mathrm{d}v=0$ , we immediately get

\varepsilon^{2a}\partial_{t}\int\int\left[\frac{1}{2}\frac{g^{2}}{\mathcal{M}}+\frac{1}{2}\eta^{2}\mathcal{M}+g\eta\right]\mathrm{d}x\mathrm{d}v=\int\int(\mathcal{L}^{s}(g)+\mathcal{L}^{s}(\eta\mathcal{M}))\left(\frac{g}{\mathcal{M}}+\eta\right)\leq 0\,.

∎

Moreover, we can bound the energy for $\eta$ and $\rho$ separately.

Proposition 2.

If $(\eta,g)$ solves (3.8) with initial data (3.9), then

E_{\eta}^{2}=\int_{\mathbb{R}}\int_{\mathbb{R}}\eta^{2}\mathcal{M}\mathrm{d}x\mathrm{d}v,\qquad E_{g}^{2}=\int_{\mathbb{R}}\int_{\mathbb{R}}\frac{g^{2}}{\mathcal{M}}\mathrm{d}x\mathrm{d}v

(3.13)

are both uniformly bounded in time.

Proof.

Multiply the second equation in (3.8) with $\eta\mathcal{M}$ and integrate in $x$ and $v$ , we get

\displaystyle\frac{\varepsilon^{2s}}{2}\partial_{t}\int\int\eta^{2}\mathcal{M}\mathrm{d}x\mathrm{d}v=-\int\int\eta(-\Delta_{x})^{s}\eta\mathrm{d}x\mathcal{M}\mathrm{d}v

From the definition (1.2), it is straightforward to see that $\int\eta(-\Delta_{x})^{s}\eta\mathrm{d}x=\frac{1}{2}\int\int\frac{(\eta(x)-\eta(y))^{2}}{|v-w|^{1+2s}}rdx\mathrm{d}y\geq 0$ . Therefore, $\frac{d}{dt}E_{\eta}\leq 0$ . Then from

E_{g}^{2}=\int\int\frac{(f-\eta\mathcal{M})^{2}}{\mathcal{M}}\mathrm{d}x\mathrm{d}v\leq 2\int\int\frac{f^{2}+\eta^{2}\mathcal{M}^{2}}{\mathcal{M}}\mathrm{d}x\mathrm{d}v=2E_{f}^{2}+2E_{\eta}^{2}\,,

we get the uniform boundedness of $E_{g}$ . ∎

Numerically, we propose the following semi-discrete scheme to (3.8):


		$\displaystyle\frac{\varepsilon^{2s}}{\Delta t}(g^{}-g^{n})=\mathcal{L}^{s}(g^{})-\gamma g^{*}-I(\eta^{n},\mathcal{M}),$			(3.14a)
		$\displaystyle\frac{\varepsilon^{2s}}{\Delta t}(g^{n+1}-g^{*})+\varepsilon v\partial_{x}g^{n+1}=\gamma g^{n+1}\,,$			(3.14b)
		$\displaystyle\frac{1}{\Delta t}(\eta^{n+1}-\eta^{n})=-(-\Delta_{x})^{s}\eta^{n}\,,$			(3.14c)

where $\gamma$ is a positive constant. As opposed to directly applying an implicit-explicit discretization to the first equation in (3.8), i.e.,

\frac{\varepsilon^{2s}}{\Delta t}(g^{n+1}-g^{n})+\varepsilon v\partial_{x}g^{n+1}=\mathcal{L}^{s}(g^{n+1})-I(\eta^{n},\mathcal{M})\,,

we conduct an operator splitting here, due to three reasons. One is that the Lévy-Fokker-Planck operator $\mathcal{L}^{s}$ has nonzero null space, and therefore the inversion in this equation will become stiff for small $\varepsilon$ (see also Table 2). The augmented term $-\gamma g^{*}$ will then shift the spectrum of the to-be-inverted operator and therefore remove the ill-conditioning. The second is that we need the magnitude of $g$ to remain small for small $\varepsilon$ in order to preserve the asymptotic property, and this requirement is fulfilled in (3.14b), which warrants $g^{n+1}$ to be of order $\varepsilon^{2s}$ as $\varepsilon\rightarrow 0$ (More details is shown in the proof of Proposition 3). The third is the computational efficiency. Thanks to the splitting, one no longer needs to invert operators in $x$ and $v$ simultaneously. Instead, only an inversion in $v$ is needed in solving (3.14a), whereas in (3.14b) only an inversion in $x$ is needed, and this inversion can be efficiently accomplished via either sweeping or the fast Fourier transform.

The rest of this subsection is devoted to proving the asymptotic property of (3.14). First let us introduce the following lemma that describes the smoothing effect of fractional diffusion equation.

Lemma 4.

Consider the initial value problem

\partial_{t}u+(-\Delta_{x})^{s}u=0\,\qquad u(0,x)=u_{in}(x)\,.

(3.15)

If $u_{in}(x)\in W^{2,1}(\mathbb{R})\cap C^{2}(\mathbb{R})$ , then $u(t,\cdot)\in W^{2,1}(\mathbb{R})\cap C^{\infty}(\mathbb{R})$ for $t>0$ .

Proof.

Using the Fourier transform for $x$ , one writes down the solution to (3.15) as

\hat{u}(t,\xi)=\hat{u}_{in}(\xi)e^{-|\xi|^{2s}t}:=\hat{u}_{in}(\xi)\hat{K}_{s}(t,\xi)\,.

Changing back to $x$ , one has

u(t,x)=\int_{\mathbb{R}}K_{s}(t,x-y)u_{0}(y)\mathrm{d}y\,,\qquad K_{s}(t,x)=\frac{1}{t^{1/2s}}F\left(\frac{|x|}{t^{1/2s}}\right)\,,

where $F$ is positive and decreasing, and it behaves like $F(r)\sim r^{-(1+2s)}$ at infinity (see also [22] for further discussion). It is easy to see $K_{s}(t,\cdot)\in L^{1}(\mathbb{R})\cap C^{\infty}(\mathbb{R})$ for $t>0$ . Therefore $u(t,\cdot)\in C^{\infty}(\mathbb{R})$ . Further, by Young’s convolution inequality, we have

\|u(t,\cdot)\|_{1}\leq\|K_{s}\|_{1}\|u_{0}\|_{1}\,,\quad\|\partial^{2}_{x}u(t,\cdot)\|_{1}\leq\|K_{s}\|_{1}\|\partial^{2}_{x}u_{0}\|_{1}\,.

∎

The proposition on the asymptotic property of the splitting scheme is in order.

Proposition 3.

Consider system $\eqref{eqn:111}$ with initial data $\left\langle f_{in}\right\rangle\in W^{2,1}(\mathbb{R})$ . Let $\eta^{n}$ and $g^{n}$ be the solution to (3.14), where $\eta^{0}=\rho_{in}(x+\varepsilon v)$ and $g^{0}=f_{in}-\eta_{in}\mathcal{M}$ . Then the numerical solution

\rho^{n}=\left\langle f^{n}\right\rangle=\left\langle\eta^{n}\mathcal{M}+g^{n}\right\rangle

satisfies

\frac{\rho^{n+1}-\rho^{n}}{\Delta t}=(-\Delta_{x})^{s}\rho^{n}

(3.16)

as $\varepsilon\rightarrow 0$ .

Proof.

From the reconstruction formula, we have

	$\displaystyle\frac{\rho^{n+1}-\rho^{n}}{\Delta t}$	$\displaystyle=\left\langle\frac{\eta^{n+1}-\eta^{n}}{\Delta t}\mathcal{M}\right\rangle+\left\langle\frac{g^{n+1}-g^{n}}{\Delta t}\right\rangle$
		$\displaystyle=-\left\langle(-\Delta_{x})^{s}\eta^{n}\mathcal{M}\right\rangle+\left\langle\frac{g^{n+1}-g^{n}}{\Delta t}\right\rangle$
		$\displaystyle=(-\Delta_{x})^{s}\left\langle f^{n}-g^{n}\right\rangle+\left\langle\frac{g^{n+1}-g^{n}}{\Delta t}\right\rangle,$
		$\displaystyle=(-\Delta_{x})^{s}\rho^{n}-(-\Delta_{x})^{s}\left\langle g^{n}\right\rangle+\left\langle\frac{g^{n+1}-g^{n}}{\Delta t}\right\rangle\,,$

where the third equality uses lemma 2 and lemma 4, namely, $\left\langle(-\Delta_{x})^{s}\eta^{n}\mathcal{M}\right\rangle=(-\Delta_{x})^{s}\left\langle\eta^{n}\mathcal{M}\right\rangle$ . Then to show (3.16), it amounts to show that the magnitude of $g^{n}$ vanishes as $\varepsilon$ approaches zero. First, from (3.14a), one sees that

\left(\frac{\varepsilon^{2s}}{\Delta t}+\gamma+\mathcal{L}^{s}\right)g^{*}=\frac{\varepsilon^{2s}}{\Delta t}g^{n}-I(\eta^{n},\mathcal{M})\,.

From the contractive estimate in Corollary 3.1 of [2] and Hille-Yosida Theorem, we have, for positive $\gamma$ ,

\|g^{*}(x,\cdot)\|_{L^{\infty}_{v}}\lesssim\|I(\eta^{n},\mathcal{M})\|_{L_{v}^{\infty}}+\frac{\varepsilon^{2s}}{\Delta t}\|g^{n}\|_{L_{v}^{\infty}}\,.

(3.17)

Further, from (3.14b), we have, for $v<0$ ,

\displaystyle g^{n+1}(x,v)=\frac{\varepsilon^{2s-1}}{v\Delta t}\int_{-\infty}^{x}e^{\frac{\frac{\varepsilon^{2s}}{\Delta t}-\gamma}{\varepsilon v}(y-x)}g^{*}(y,v)\mathrm{d}y

Therefore,

\displaystyle|g^{n+1}(x,v)|\leq\frac{\varepsilon^{2s-1}}{|v|\Delta t}\|g^{*}(\cdot,v)\|_{L^{\infty}_{x}}\int_{-\infty}^{x}e^{\frac{\frac{\varepsilon^{2s}}{\Delta t}-\gamma}{\varepsilon v}(y-x)}\mathrm{d}y\leq\frac{2\varepsilon^{2s}}{\Delta t}\|g^{*}(\cdot,v)\|_{L^{\infty}_{x}}\,.

(3.18)

The case with $v>0$ can be estimated similarly. Then a combination of (3.17) and (3.18) immediately leads to $\|g^{n+1}\|_{L^{\infty}_{x,v}}\lesssim\mathcal{O}(\varepsilon^{2s})$ . ∎

3.2 Fully discrete scheme and implementation

Now we briefly discuss the implementation of scheme $\eqref{eqn: semi_scheme}$ . Using the same notation as mentioned at the beginning of this section, we denote the numerical approximations $\eta^{n}_{i,j}$ , $g^{n}_{i,j}$ as $\eta^{n}_{i,j}\approx\eta(t^{n},x_{i},q_{j})$ and $g^{n}_{i,j}\approx g(t^{n},x_{i},q_{j})$ , with $0\leq i\leq N_{x}-1$ , $0\leq j\leq N_{v}-1$ . For a fixed $x$ -index $i$ , let

	$\displaystyle{\bf\eta^{n}_{i}}=(\eta^{n}_{i,1},\eta^{n}_{i,2},\cdots,\eta^{n}_{i,N_{v}-1})^{T},\qquad{\bf g^{n}_{i}}=(g^{n}_{i,1},g^{n}_{i,2},\cdots,g^{n}_{i,N_{v}-1})^{T}\,,\qquad$
	$\displaystyle{\bf\mathcal{M}}=(\mathcal{M}_{1},\cdots,\mathcal{M}_{N_{v}-1})^{T}\,,\qquad{\bf\eta^{n}_{i}\mathcal{M}}=(\eta^{n}_{i,1}\mathcal{M}_{1},\cdots,\mathcal{\eta}^{n}_{i,N_{v}-1}M_{N_{v}-1})^{T}\,.$

For a fixed $q$ -index $j$ , let

{\bf\eta^{n}_{j}}=(\eta^{n}_{1,j},\eta^{n}_{2,j},\cdots,\eta^{n}_{N_{x}-1,j})^{T},\qquad{\bf g^{n}_{j}}=(g^{n}_{1,j},g^{n}_{2,j},\cdots,g^{n}_{N_{x}-1,j})^{T}\,.

First we compute $I(\eta^{n},\mathcal{M})$ . According to its definition

I(\eta^{n},\mathcal{M})=(-\Delta_{v})^{s}(\eta^{n}\mathcal{M})-\mathcal{M}(-\Delta_{v})^{s}\eta^{n}-\eta^{n}(-\Delta_{v})^{s}\mathcal{M},

this is simply accomplished by applying $\mathsf{L}^{s}$ defined in (2.9) to ${\bf\eta_{i}^{n}\mathcal{M}}$ , ${\bf\eta_{i}^{n}}$ , and ${\bf\mathcal{M}}$ , respectively, for a fixed spatial index $i$ . Then to treat the spatial discretization in (3.14b) and (3.14c), we use the Fourier spectral method with periodic boundary condition. Note that the transport term in (3.14b) is treated implicitly, so the stability is guaranteed. Moreover, unlike in $v$ direction, where one needs to consider a fat tail distribution due to the kernel of $\mathcal{L}^{s}$ , we can consider a sufficiently decaying profile in $x$ and therefore the Fourier spectral method can be used here.

To summarize, for given $\eta^{n}_{i,j}$ and $g^{n}_{i,j}$ , we have

Step 1 Compute $\bf g_{i}^{*}$ by

{\bf g_{i}^{*}}=\left[\left(\frac{\varepsilon^{2s}}{\Delta t}+\gamma\right)\mathsf{I}-\mathsf{P}^{s}\right]^{-1}\left(\frac{\varepsilon^{2s}}{\Delta t}{\bf g_{i}^{n}}-I({\bf\eta_{i}}^{n}\mathcal{M})\right)\,,\qquad{\bf i}=0,1,\cdots,N_{x}-1\,;

Step 2 Compute ${\bf g_{j}^{n}}$ by inverting

\varepsilon^{2s}\frac{{\bf\hat{g}_{j}^{n+1}}-{\bf\hat{g}_{j}^{*}}}{\Delta t}+\varepsilon v(q_{j})\mathsf{D}_{x}{\bf\hat{g}_{j}^{n+1}}=\frac{1}{2}{\bf\hat{g}_{j}^{n+1}}\,,\qquad{\bf j}=0,1,\cdots,N_{v}-1\,,

where ${\bf\hat{g}_{j}}$ is the Fourier transform of ${\bf g_{j}}$ , and $\mathsf{D}_{x}$ is the diagonal matrix with elements $k=0,1,\cdots,N_{x}-1$ .

Step 3 Compute ${\bf\eta_{0}^{n+1}}$ via

{\bf\hat{\eta}_{0}^{n+1}}=(\mathsf{I}+\Delta t\mathsf{D}_{s})^{-1}{\bf\hat{\eta}_{0}^{n}},\qquad{\bf j}=0,1,\cdots,N_{v}-1\,,

where ${\bf\hat{\eta}_{0}}$ is the Fourier transform of ${\bf\eta_{0}}$ , and $\mathsf{D}_{s}$ is the diagonal matrix with elements $-|k|^{2s}$ with $k$ ranging from 0 to $N_{x}-1$ .

4 Numerical examples

In this section, we present several numerical examples to check the accuracy and efficiency of schemes $\eqref{eqn: homo_scheme_ori}$ and $\eqref{eqn: semi_scheme}$ . Periodic boundary condition is always used in $x$ direction. Unless otherwise specified, we choose $L_{v}=3$ and $l_{lim}=300$ in computing (2.7), and $\gamma=1$ in (3.14).

4.1 Computation of $(-\Delta)^{s}f$

We first check the performance of the pseudospectral method in computing $(-\Delta)^{s}f$ via (2.9). Two specific examples will be considered here: an exponential decay function and a polynomial decay function, both of which have exact form when taken the fractional Laplacian.

$\displaystyle 1)\qquad$	$\displaystyle f(v)=(1+v^{2})^{-\frac{1-2s}{2}},\quad s\in(0,0.5),$	(4.1)
	$\displaystyle(-\Delta_{v})^{s}f(v)=-2^{2s}\Gamma\left(\frac{1+2s}{2}\right)\Gamma\left(\frac{1-2s}{2}\right)^{-1}(1+v^{2})^{-\frac{1+2s}{2}}\,.$
$\displaystyle 2)\qquad$	$\displaystyle f(v)=e^{-v^{2}}\,,$	(4.2)
	$\displaystyle(-\Delta_{v})^{s}f(v)=-\frac{1}{\sqrt{\pi}}2^{2s}\Gamma\left(\frac{1+2s}{2}\right){}_{1}F_{1}\left(\frac{1+2s}{2},0.5;-v^{2}\right)\,.$

Here $\Gamma$ is the gamma function and ${}_{1}F_{1}$ is the hypergeometric function.

First, a comparison of the numerical solutions with exact solutions is gathered in Fig. 1 and 2, for (4.1) and (4.2), respectively, where good agreements are observed in both cases.

Refer to caption — Figure 1: Comparison of numerical approximation and exact expression of $(-\Delta)^{s}f$ for $f$ in (4.1). $N_{v}=128$ are used for both cases.

Second, we show relationship between accuracy and the number of modes $N_{v}$ . Given fixed $L_{v}=3$ and $l_{lim}=300$ , Fig. 3 displays the error versus $N_{v}$ for different $s$ , with exponential decay function on the left, and power law decay function on the right. The error is measured in $l^{\infty}$ norm, i.e.

e_{\infty}=\max_{j}\|u_{j}-(-\Delta_{v})^{s}f(v_{j})\|_{\infty}

(4.3)

where $u_{j}$ is the numerical approximation of $(-\Delta_{v})^{s}f$ at $v_{j}$ . One sees that a small number of mode is adequate for the exponential decay case, whereas for the power law decay case, increasing the number of modes leads to better approximation. Moreover, a comparison between these two figures also gives a visual indication on why computing the fractional Laplacian of a slow decaying function is significantly harder than a fast decaying function, as a lot more modes are needed to reach an even lower accuracy criteria.

4.2 Spatially homogeneous case

We restrict our attention to the spatially homogeneous case in this section and consider the following specific example:

\partial_{t}f=\partial_{v}(vf)-(-\Delta_{v})^{s}f\,,\qquad f(0,v)=e^{-v^{2}}\,.

(4.4)

In order to check the performance of the scheme, two properties of the solution will be considered: 1) convergence towards equilibrium in the long time limit and 2) mass conservation.

4.2.1 Long time behavior

As pointed out in [9], $f(t,v)$ will converge toward the equilibrium $\mathcal{M}$ exponentially fast. To observe this dynamics numerically, we compute the relative entropy $\int_{\mathbb{R}}f\ln\frac{f}{\mathcal{M}}\mathrm{d}v$ as

H^{n}=\sum_{j=0}^{N_{v}-1}f^{n}(q_{j})\ln\frac{f^{n}(q_{j})}{\mathcal{M}(q_{j})}w_{j},\qquad w_{j}=\frac{L_{v}}{(\sin(q_{j}))^{2}}\,.

(4.5)

The numerical equilibrium, denoted as $f_{\infty}$ , is obtained when the variation in the solution is negligible, i.e.,

\sum_{j}|f_{j}^{n+1}-f_{j}^{n}|w_{j}\Delta q\leq\delta,\qquad w_{j}=\frac{L_{v}}{(\sin(q_{j}))^{2}}\,,

(4.6)

and then set $f_{\infty}=f^{n+1}$ . Here $\delta$ is a small parameter and chosen to be $10^{-6}$ in our examples.

We first test the case with $s=0.5$ , where an explicit form of equilibrium is available: $\mathcal{M}(v)=\pi^{-0.5}(1+v^{2})^{-1}$ . The results are collected in Fig 4. In the left and middle figures, one sees that the numerical equilibrium coincides with analytic one, both in the bulk area and in the tail. On the right, the evolution of relative entropy (4.5) in time is plotted, and exponential rate of convergence is confirmed.

Next, we test two other cases $s=0.6,~{}0.8$ , where no explicit form of equilibrium is available, therefore we only check the tail behavior, as predicted in (1.3). As shown in Fig. 4, correct power law decay rate at tail is captured numerically and exponential convergence toward equilibrium is observed.

4.2.2 Mass conservation

In this section, we check the total mass versus time. As mentioned in Remark 2, the total mass is not exactly conserved by our scheme, but with proper choice of $L_{v}$ and $N_{v}$ , the error can be controlled. Consider the same initial value problem as in (4.4), where the total mass is $M_{0}=\pi^{\frac{1}{2}}$ , we define the following error in mass:

M^{n}=|\sum_{j=1}^{N_{v}}f_{N}^{n}(q_{j})w_{j}\Delta q-M_{0}|,\qquad w_{j}=\frac{L_{v}}{(\sin(q_{j}))^{2}}\,,

(4.7)

and plot it with $s=0.4,~{}0.6,~{}0.8$ in Fig. 6. It is shown that, with the same choice of $L_{v}=3$ , $N_{v}=128$ , the mass is conserved better for larger $s$ , which further indicates that slower decaying function is harder to compute numerically. We also check the mass error at the numerical equilibrium (defined in (4.6)) for different $s$ with increasing $N_{v}$ . The results are shown in Table 1. As expected, with fixed $L_{v}$ , larger $N_{v}$ leads to better conservation.

Table 1: Mass error at the numerical equilibrium

f_{\infty}^{N_{v}}

$N_{v}$	64	128	256	512
$M^{\infty}$ with $s=0.4$	7.8e-2	5.9e-3	2.1e-3	9.5e-4
$M^{\infty}$ with $s=0.5$	5.7e-2	2.2e-3	5.3e-4	1.3e-4
$M^{\infty}$ with $s=0.6$	3.3e-3	5.9e-4	9.2e-5	5.3e-6
$M^{\infty}$ with $s=0.8$	2.9e-4	3.2e-5	9.4e-7	1.1e-6

4.3 Spatially inhomogeneous case

Throughout this section, we compute $\eqref{eqn:111}$ with different choices of $\varepsilon$ . The following two initial conditions are considered:

f(0,x,v)=\pi^{-0.5}(1+\sin(\frac{\pi}{L_{x}}x))e^{-v^{2}},\quad x\in[-\pi,\pi]\,,

(4.8)

and

f(0,x,v)=\pi^{-0.5}e^{-15x^{2}}e^{-v^{2}}\,,\quad x\in[-5,5]\,.

(4.9)

The equilibrium $\mathcal{M}$ , except for $s=1/2$ , is obtained numerically by running the spatially homogeneous solver until converge, as described in Section 4.2.1. Note that the periodic initial data (4.8) does not exactly fall into the assumption of Lemma 2, but the conclusion there still holds. Indeed, for function $f(x,v)$ of the form (4.8), one can write it into Fourier series, then it amounts to check $\int_{\mathbb{R}}\int_{\Omega}f(v,x)e^{-i\xi x}\mathrm{d}x\mathrm{d}v=\int_{\Omega}\int_{\mathbb{R}}f(v,x)e^{-i\xi x}\mathrm{d}v\mathrm{d}x$ , which is true as long as $f\in L^{1}(\mathbb{R}\times\Omega)$ .

4.3.1 Advantage of splitting in (3.14)

Before presenting specific numerical examples, we first verify the advantage of splitting and adding the term $-\gamma g^{*}$ in (3.14a). In particular, from (3.14a), one updates $g^{*}$ as $g^{*}=\mathsf{A}_{\gamma}^{-1}(g^{n}-\Delta t\varepsilon^{-2s}I(\tilde{\rho}^{n},\mathcal{M}))$ , where $\mathsf{A}_{\gamma}=[\mathsf{I}-\Delta t\varepsilon^{-2s}(\mathsf{L}^{s}-\gamma\mathsf{I})]$ . Then the success of this step hinges on the condition number of $\mathsf{A}_{\gamma}$ . In Table 2, we compute the condition number of $\mathsf{A}_{\gamma}$ with $\gamma=0,\frac{1}{2},1,2$ . The other parameters used are $L_{x}=\pi$ , $N_{x}=50$ , $L_{v}=3$ , $N_{v}=64$ and $\Delta t=0.1$ . One sees that larger $\gamma$ leads to a better conditioned $\mathsf{A}_{\gamma}$ , and this improvement is more pronounced for smaller $\varepsilon$ . Therefore, in the fractional diffusive regime, $\gamma$ plays a non-negligible role.

Table 2: Condition number of

\mathsf{A}_{\gamma}

	$Cond(\mathsf{A}_{0})$	$Cond(\mathsf{A}_{\frac{1}{2}})$	$Cond(\mathsf{A}_{1})$	$Cond(\mathsf{A}_{2})$
$s=0.4$ , $\varepsilon=1$	4.73	4.54	4.36	4.06
$s=0.6$ , $\varepsilon=1$	4.93	4.73	4.55	4.24
$s=0.8$ , $\varepsilon=1$	13.28	12.66	12.11	11.14
$s=0.4$ , $\varepsilon=1e-3$	7.74e3	158	52	20.90
$s=0.6$ , $\varepsilon=1e-3$	1.24e5	187.42	63.44	25.45
$s=0.8$ , $\varepsilon=1e-3$	5.97e6	564.32	190.75	75.81
$s=0.4$ , $\varepsilon=1e-5$	3.57e5	181.66	55.84	21.43
$s=0.6$ , $\varepsilon=1e-5$	3.15e7	188.98	63.67	25.49
$s=0.8$ , $\varepsilon=1e-5$	9.46e9	564.62	190.79	75.82

4.3.2 Uniform accuracy in time

In this section, we show the first order accuracy in time by computing the following error

e_{\Delta t}=\sum_{i}^{N_{x}}\sum_{j}^{N_{v}}|f^{\Delta t}_{i,j}(T)-f^{\Delta t/2}_{i,j}(T)|w_{j}\Delta q\Delta x.

In Fig. 7, first order accuracy in time is shown among different choices of $eps$ : $\varepsilon=1$ , $\varepsilon=1e-3$ and $\varepsilon=1e-5$ .

4.3.3 Energy stability

We compute the total energy, and the individual energy corresponding to the macro and micro part in this section. More precisely, the numerical approximation to (3.12) and (3.13) are

E_{f}=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{v}}\frac{f_{i,j}^{2}}{\mathcal{M}_{j}}w_{j}\Delta q\Delta x\,,\quad E_{g}=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{v}}\frac{g_{i,j}^{2}}{\mathcal{M}_{j}}w_{j}\Delta q\Delta x\,,\quad E_{\eta}=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{v}}\eta_{i,j}^{2}\mathcal{M}_{j}w_{j}\Delta q\Delta x\,.

In Fig. 8, we compute (1.5) with initial condition $\eqref{eqn: IC2}$ and various $\varepsilon$ , ranging from $1$ to $1e-5$ . The solution is computed up to $T=0.1$ . As shown, $E_{f}$ decays with time, and $E_{g}$ and $E_{\eta}$ are uniformly bounded in time, which coincides with the Propositions 1 and 2 and confirms the stability of our scheme. When $\varepsilon=1e-5$ , the sudden change at the first step is due to fact that (4.9) is not at the equilibrium.

4.3.4 Kinetic regime $\varepsilon=1$

We then check the performance of our scheme in the kinetic regime with $\varepsilon=1$ . As a comparison, the following implicit-explicit (IMEX) scheme is used on finer mesh to produce the reference solution:

\displaystyle\begin{cases}{}\frac{f^{*}-f^{n}}{\Delta t}+v\partial_{x}f^{n}=0,\\ \frac{f^{n+1}-f^{*}}{\Delta t}=\partial_{v}(vf^{n+1})-(-\Delta_{v})^{s}f^{n+1}.\end{cases}

Here Fig. 9 and Fig. 10 correspond to initial conditions (4.8) and (4.9), respectively. Different choices of $s$ are considered. One sees that the numerical solution from our AP scheme agrees well with the reference solution.

4.3.5 AP property and diffusive regime

In this section, we check the AP property of the scheme and test its performance in the diffusive regime. To check the AP property, we compute the following asymptotic error

Error=\|\rho^{\varepsilon}\mathcal{M}-{f}^{\varepsilon}\|_{1}=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{v}}|\rho^{\varepsilon}_{i}\mathcal{M}_{j}-{f^{\varepsilon}_{i,j}}|w_{j}\Delta q\Delta x\,,

(4.10)

where $\mathcal{M}$ is the equilibrium and $\rho^{\varepsilon}=\int_{\mathbb{R}}f^{\varepsilon}\mathrm{d}v$ .

When $\varepsilon=1e-5$ , we compare the solution to our AP scheme with the solution to the diffusive equation (1.6), which is computed using the Fourier spectral method.

The results are collected in Fig. 11 and Fig. 12, for initial data (4.8) and (4.9), respectively. The left column of each figure represents the asymptotic error (4.10) in time with different choice of $\varepsilon$ . It is clearly that the error decreases with vanishing $\varepsilon$ , and it is at a magnitude of approximately $\mathcal{O}(\varepsilon)$ . On the right, a good match between the solution to the AP scheme with the solution to the fractional limit is observed. Two different $s$ are tested for each cases.

5 Conclusion

We designed an asymptotic preserving scheme for Lévy-Fokker-Planck equation with fractional diffusion limit. This limit emerges due to the fat tail equilibrium of the Lévy-Fokker-Planck operator, which breaks down the classical diffusion limit as it renders the diffusion matrix unbounded. Similar anomalous diffusion was considered for the linear Boltzmann case [19, 1], for which asymptotic preserving schemes have been designed [23, 7, 8, 24]. Comparing to the linear Boltzmann case, there are two major difficulties here in constructing numerical schemes. One is that the fat tail equilibrium does not appear explicitly in the collision operator, but exits implicitly as the kernel of the collision operator. Therefore, the idea of truncating the infinite domain into a finite computational one with a tail compensation can not be directly apply, as the tail behavior is not known unless the solution reaches the equilibrium. The other comes from the derivation of the fractional diffusion limit. In the linear Boltzmann case, a reshuffled Hilbert expansion is performed to show the strong convergence of the kinetic equation to the anomalous diffusion limit and it is the building block of the design of AP scheme. In contrast, only a weak convergence is known for our case. To resolve the first difficulty, we adopt a pseudo spectral method based on rational Chebyshev polynomial, which transforms an infinite domain into a finite one, and therefore no domain decomposition is needed. For the second difficulty, we propose a novel macro-micro decomposition, with a unique macro part that is inspired by the special choice of of the test function in proving the weak convergence. The stability of the split system is obtained. We also propose an operator splitting discretization to the split system, which removes the ill-posedness due to the stiffness and reduces the computational cost from a direct implicit treatment. A rigorous asymptotic preserving property of our scheme is established. Several numerical results are carried out to demonstrate the properties of our scheme, including asymptotic-preservation, uniform accuracy and energy dissipation.

Acknowledgment: This work is partially supported by NSF grant DMS-1846854. L.W. would like to thank Dr. Min Tang and Dr. Jingwei Hu for the discussion on computing the fractional Laplacian operator.

References

[1] N. B. Abdallah, A. Mellet, and M. Puel, Fractional diffusion limit for collisional kinetic equations: a hilbert expansion approach, Kinetic & Related Models, 4 (2011), p. 873.
[2] P. Biler and G. Karch, Generalized Fokker-Planck equations and convergence to their equilibria, Banach Center Publications, 60 (2003), pp. 307–318.
[3] A. Bonito, J. P. Borthagaray, R. H. Nochetto, E. Otárola, and A. J. Salgado, Numerical methods for fractional diffusion, Computing and Visualization in Science, 19 (2018), pp. 19–46.
[4] J. P. Boyd, Spectral methods using rational basis functions on an infinite interval, Journal of Computational Physics, 69 (1987), pp. 112–142.
[5] J. Cayama, C. M. Cuesta, and F. de la Hoz, A pseudospectral method for the one-dimensional fractional laplacian on R, 2019.
[6] L. Cesbron, A. Mellet, and K. Trivisa, Anomalous transport of particles in plasma physics, Applied Mathematics Letters, 25 (2012), pp. 2344–2348.
[7] N. Crouseilles, H. Hivert, and M. Lemou, Numerical schemes for kinetic equations in the anomalous diffusion limit. part i: The case of heavy-tailed equilibrium, SIAM Journal on Scientific Computing, 38 (2016), pp. A737–A764.
[8] , Numerical schemes for kinetic equations in the anomalous diffusion limit. part ii: Degenerate collision frequency, SIAM Journal on Scientific Computing, 38 (2016), pp. A2464–A2491.
[9] I. Gentil and C. Imbert, The Lévy–Fokker–Planck equation: $\phi$ -entropies and convergence to equilibrium, Asymptotic Analysis, 59 (2008), pp. 125–138.
[10] J. Hu, S. Jin, and Q. Li, Asymptotic-preserving schemes for multiscale hyperbolic and kinetic equations, in Handbook of Numerical Analysis, vol. 18, Elsevier, 2017, pp. 103–129.
[11] Y. Huang and A. Oberman, Numerical methods for the fractional laplacian: A finite difference-quadrature approach, SIAM Journal on Numerical Analysis, 52 (2014), pp. 3056–3084.
[12] , Finite difference methods for fractional Laplacians, arXiv preprint arXiv:1611.00164, (2016).
[13] S. Jin, Asymptotic preserving (ap) schemes for multiscale kinetic and hyperbolic equations: a review, Lecture notes for summer school on methods and models of kinetic theory (M&MKT), Porto Ercole (Grosseto, Italy), (2010), pp. 177–216.
[14] M. Kwaśnicki, Ten equivalent definitions of the fractional laplace operator, Fractional Calculus and Applied Analysis, 20 (2017), pp. 7–51.
[15] M. Lemou and L. Mieussens, A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit, SIAM Journal on Scientific Computing, 31 (2008), pp. 334–368.
[16] A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M. M. Meerschaert, M. Ainsworth, et al., What is the fractional laplacian? a comparative review with new results, Journal of Computational Physics, 404 (2020), p. 109009.
[17] H. Ma, W. Sun, and T. Tang, Hermite spectral methods with a time-dependent scaling for parabolic equations in unbounded domains, SIAM journal on numerical analysis, 43 (2005), pp. 58–75.
[18] Z. Mao and J. Shen, Hermite spectral methods for fractional PDEs in unbounded domains, SIAM Journal on Scientific Computing, 39 (2017), pp. A1928–A1950.
[19] A. Mellet, S. Mischler, and C. Mouhot, Fractional diffusion limit for collisional kinetic equations, Archive for Rational Mechanics and Analysis, 199 (2011), pp. 493–525.
[20] F. Poupaud and J. Soler, Parabolic limit and stability of the vlasov–fokker–planck system, Mathematical Models and Methods in Applied Sciences, 10 (2000), pp. 1027–1045.
[21] C. Sheng, J. Shen, T. Tang, L.-L. Wang, and H. Yuan, Fast Fourier-like mapped Chebyshev spectral-Galerkin methods for pdes with integral fractional Laplacian in unbounded domains, SIAM Journal on Numerical Analysis, 58 (2020), pp. 2435–2464.
[22] J. L. Vázquez, Recent progress in the theory of nonlinear diffusion with fractional Laplacian operators, arXiv preprint arXiv:1401.3640, (2014).
[23] L. Wang and B. Yan, An asymptotic-preserving scheme for linear kinetic equation with fractional diffusion limit, Journal of Computational Physics, 312 (2016), pp. 157–174.
[24] , An asymptotic-preserving scheme for the kinetic equation with anisotropic scattering: Heavy tail equilibrium and degenerate collision frequency, SIAM Journal on Scientific Computing, 41 (2019), pp. A422–A451.

	$\displaystyle I(x,v)$	$\displaystyle=\int_{0}^{\infty}\|(2h(x,v)-h(x-y,v)-h(x+y,v))\nu(y)\|dy$
		$\displaystyle=\int_{0}^{\delta}\|(2h(x,v)-h(x-y,v)-h(x+y,v))\nu(y)\|dy$
		$\displaystyle\qquad\qquad+\int_{\delta}^{\infty}\|(2h(x,v)-h(x-y,v)-h(x+y,v))\nu(y)\|dy$
		$\displaystyle:=I_{1}(x,v)+I_{2}(x,v)\,,$

An asymptotic preserving scheme for Lévy-Fokker-Planck equation with fractional diffusion limit

Abstract

1 Introduction

Theorem 1.

2 Computation of the collision operator ℒs\mathcal{L}^{s}

2.1 Change of variable

2.2 Computation of (−Δq)s​f(-\Delta_{q})^{s}f

Theorem 2.

2.3 Spatially homogeneous case

Remark 1 (Positivity).

Remark 2 (Choice of LvL_{v} and NvN_{v}).

3 Asymptotic preserving scheme

Lemma 1.

Proof.

Lemma 2.

Proof.

3.1 A novel micro-macro decomposition

Lemma 3.

Proof.

Remark 3.

Proposition 1.

Proof.

Proposition 2.

Proof.

Lemma 4.

Proof.

Proposition 3.

Proof.

3.2 Fully discrete scheme and implementation

4 Numerical examples

4.1 Computation of (−Δ)s​f(-\Delta)^{s}f

4.2 Spatially homogeneous case

4.2.1 Long time behavior

4.2.2 Mass conservation

4.3 Spatially inhomogeneous case

4.3.1 Advantage of splitting in (3.14)

4.3.2 Uniform accuracy in time

4.3.3 Energy stability

4.3.4 Kinetic regime ε=1\varepsilon=1

4.3.5 AP property and diffusive regime

5 Conclusion

References

2 Computation of the collision operator $\mathcal{L}^{s}$

2.2 Computation of $(-\Delta_{q})^{s}f$

Remark 2 (Choice of $L_{v}$ and $N_{v}$ ).

4.1 Computation of $(-\Delta)^{s}f$

4.3.4 Kinetic regime $\varepsilon=1$