A finite-volume scheme for fractional diffusion on bounded domains

Abstract

We propose a new fractional Laplacian for bounded domains, expressed as a conservation law and thus particularly suited to finite-volume schemes. Our approach permits the direct prescription of no-flux boundary conditions. We first show the well-posedness theory for the fractional heat equation. We also develop a numerical scheme, which correctly captures the action of the fractional Laplacian and its anomalous diffusion effect. We benchmark numerical solutions for the Lévy–Fokker–Planck equation against known analytical solutions. We conclude by numerically exploring properties of these equations with respect to their stationary states and long-time asymptotics.

1. Introduction

The aim of this work is the design of a finite-volume numerical scheme to approximate the solution of the non-local diffusion problem given by the fractional heat equation and the related Lévy–Fokker–Planck equation. The fractional heat equation is defined in $\mathbb{R}^d$ as

(1.1) \begin{equation} \frac{\partial \rho }{\partial t} = - ({-\Delta })^{\frac{\alpha }{2}}\rho \end{equation}

for $0 \lt \alpha \leq 2$ . The so-called fractional Laplacian, $({-\Delta })^{\frac{\alpha }{2}}\rho$ , can be formally defined by its Fourier symbol $|{\xi }|^\alpha \hat \rho$ , although it admits up to ten equivalent definitions (see [28]). There is a suitable self-similar change of variables that leads to $\frac{\partial \rho }{\partial t} = \nabla \cdot (x\rho ) - ({-\Delta })^{\frac{\alpha }{2}}\rho$ , a particular case of Lévy–Fokker–Planck equation given by

(1.2) \begin{equation} \frac{\partial \rho }{\partial t} = \nabla \cdot (\beta x\rho ) - ({-\Delta })^{\frac{\alpha }{2}}\rho \end{equation}

for $0 \lt \alpha \leq 2$ and $\beta \geq 0$ . Notice that this equation generalises the usual Fokker–Planck equation $\frac{\partial \rho }{\partial t} = \nabla \cdot (\beta x\rho ) + \Delta \rho$ by replacing the Laplacian with a fractional operator, see [8, 24].

Fractional diffusion (in particular, the fractional Laplacian) has been shown to be the mean field limit of Lévy walks under certain scalings [40]. This kind of stochastic process consists in the random movement of particles in space, subject to a probability that allows long jumps with a polynomial tail. Such random walks are long range stochastic processes, and they are generally considered more realistic in the modelling of certain biological phenomena [11, 21, 22, 29, 32].

The inverse Fourier transform of the symbol $|{\xi }|^\alpha \hat \rho$ yields, after some work, the Riesz or singular integral definition of the fractional Laplacian:

\begin{align*} ({-\Delta })^{\frac{\alpha }{2}} \rho (x) \,:\!=\, \mathcal{C}({d,\alpha }) \;\textrm{p.v.} \int _{\mathbb{R}^d} \frac{\rho (x)-\rho (y)}{|{x-y}|^{d+\alpha }} \mathrm{d} y, \end{align*}

where the integral is understood in the Cauchy principal value sense in order to overcome the singularity. The constant $\mathcal{C}({d,\alpha })$ , a term which arises in the computation of the inverse transform of $|{\xi }|^\alpha$ , is given by

\begin{align*} \mathcal{C}({d,\alpha }) = \frac{ 2^\alpha \Gamma \left ({\frac{d+\alpha }{2}}\right ) }{ \pi ^{\frac{d}{2}} |{\Gamma \left ({-\frac{\alpha }{2}}\right )}| }. \end{align*}

Using the Riesz potential, equation (1.2) can be formally written in divergence form as

(1.3) \begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot F = 0, \quad \text{where } F = -\beta x\rho + \nabla \left [{({-\Delta })^{\frac{\alpha }{2}-1} \rho }\right ]. \end{equation}

The advantage of this form is that the fractional operator now appears with a negative exponent. In this case, the inverse Fourier transform of the symbol $|{\xi }|^{-\alpha } \hat \rho$ yields (see [38, Chapter 5])

(1.4) \begin{align} ({-\Delta })^{-\frac{\alpha }{2}} \rho (x) = \mathcal{C}({d,-\alpha }) \int _{\mathbb{R}^d} \frac{\rho (y)}{|{x-y}|^{d-\alpha }} \mathrm{d} y \end{align}

whenever $0\lt \alpha \lt d$ . This form for the inverse operator bypasses the singularity altogether. Equation 1.3 can therefore be rewritten as

(1.5) \begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot F = 0, \quad \text{where } F = -\beta x\rho + \mathcal{C}({d,\alpha -2}) \nabla \int _{\mathbb{R}^d} \frac{\rho (y)}{|{x-y}|^{d+\alpha -2}} \mathrm{d} y \end{equation}

whenever $\alpha \gt 2-d$ . Therefore, in dimension one, this formulation is only valid for $1 \lt \alpha \leq 2$ ; in higher dimensions, for $0 \lt \alpha \leq 2$ . This new form of the equation has two advantages: the first, that the fractional operator is no longer singular; and the second, that an equation in divergence form lends itself to be discretised in the finite-volume fashion. Finite-volume schemes have been used with success to produce structure preserving schemes for equations in divergence form of gradient flow type and related systems, see [5, 6, 15, 16]. This is a departure from the numerical methods for fractional diffusions that have been developed in the past, where the literature has been focused on finite-element and finite-difference methods [3, 4, 10, 18, 20, 25, 33, 35]. We also highlight several spectral methods [17, 34, 37, 41], which deal exclusively with problems on unbounded domains.

For the sake of computation, we would like to pose equation (1.2) on an open bounded domain $\Omega \subset{\mathbb{R}^d}$ . There are several nonequivalent definitions of fractional-type Laplacians on bounded domains that can be obtained as suitable restrictions of the definitions in $\mathbb{R}^d$ (see [1, Section 1.2] and the references therein). In this work, we construct a new fractional Laplacian by restricting equation (1.5) to the domain $\Omega$ , prescribing zero-flux conditions for the divergence, and extending the density as $\rho \equiv 0$ on ${\mathbb{R}^d}\setminus \Omega$ in order to ensure that the non-local operator is well-defined. Thus, our interpretation of the Lévy–Fokker–Planck equation on a bounded domain is

(1.6) \begin{equation} \left \{ \begin{aligned} & \frac{\partial \rho }{\partial t} + \nabla \cdot F = 0, \\ & F = -\beta x\rho + \mathcal{C}({d,\alpha }) \nabla \int _\Omega \frac{\rho (y)}{|{x-y}|^{d+\alpha -2}} \mathrm{d} y, \\ & F \cdot n_{\partial \Omega } = 0. \end{aligned} \right. \end{equation}

In the absence of a drift term (when $\beta =0$ ), we also obtain an interpretation of the fractional heat equation on a bounded domain. Notice, however, that the steady states of this problem will satisfy

\begin{equation*} \int _\Omega \frac{\rho _\infty (y)}{|{x-y}|^{d+\alpha -2}} \mathrm{d} y = C, \end{equation*}

which causes $\rho _\infty$ to be singular on the boundary. However, for $\beta \gt 0$ , our numerical results on suitably scaled quadrangular domains show that the numerical steady state is very similar to the self-similar profile of the $\mathbb{R}^d$ case (see subsections 4.1.3 and 4.2.1).

The rest of this work is organised as follows: in section 2, we study the well-posedness of equation (1.6), distinguishing the cases $\beta =0$ and $\beta \gt 0$ ; in section 3, we introduce a finite-volume numerical scheme for equation (1.6) in one dimension, and then generalise it to higher dimensions via dimensional splitting; we conclude in section 4 by validating our schemes against known analytical results.

2. A new fractional Laplacian in bounded domains

The aim of this section is to establish a well-posedness theory for the new fractional operator on a bounded domain introduced above. We first define the Riesz kernel and the extension operator

\begin{equation*} \mathcal I_{\gamma } [u] (x) = \mathcal{C}(d, -\gamma ) \int _{\mathbb{R}^{d}} \frac{u(y)}{|x-y|^{d-\gamma }}d y, \qquad \qquad \mathcal E [u] (x) = \begin{cases} u(x) & \text{if } x \in \Omega, \\ 0 & \text{otherwise}, \end{cases} \end{equation*}

where we assume that $0 \lt \gamma \lt d$ for the operator to be well defined. The operator $\mathcal I_{\gamma } \,:\, L^p (\mathbb R^d) \to L^q(\mathbb R^d)$ has been widely studied. Note that the flux in equation (1.6) reduces, whenever $\beta = 0$ , to

\begin{equation*} F = \nabla \Big ( \mathcal I_{2\left(1-\frac{\alpha }{2}\right)} \mathcal E u \Big ). \end{equation*}

Thus, besides $\alpha \in (0,2)$ , we require $0 \lt 2\left(1-\frac{\alpha }{2}\right) \lt d$ for well-posedness, i.e. $\alpha \gt 2-d$ . This is only restrictive in dimension $d = 1$ .

Let $\mathcal B = \mathcal I_{2\left(1-\frac{\alpha }{2}\right)} \mathcal E$ . If $\Omega ={\mathbb{R}^d}$ , then $\mathcal B = ({-}\Delta )^{-(1-\frac \alpha 2)}$ , an inverse fractional Laplacian. Our new diffusion operator is therefore $\mathcal A = -\Delta \mathcal B u$ , and we denote its domain by $D(\mathcal{A})$ . Equation 1.6 can be now written in the $\beta = 0$ case as

(2.1) \begin{equation} \begin{cases} \dfrac{\partial u}{\partial t} = \Delta \mathcal B u & \quad \text{in } (0,\infty ) \times \Omega, \\[15pt] \dfrac{\partial (\mathcal B u) }{\partial n} = 0 & \quad \text{on } (0,\infty ) \times \Omega. \end{cases} \end{equation}

Theorem 2.1. There exists a unique semigroup $\mathcal S(t) \,:\, L^2 (\Omega ) \to L^2 (\Omega )$ of solutions of (2.1). In fact, if $\mathcal B u_0 \in H^2(\Omega )$ then $\mathcal B u(t) \in H^2 (\Omega )$ for all times and the equation is satisfied in the operator sense.

Proof. We will first justify the well-posedness of this problem using the Hille–Yosida theorem applied to the operator $\mathcal A$ with the Neumann condition in $L^2 (\Omega )$ (see [12, Theorem 7.4]).

Let us construct $D(\mathcal A)$ . We begin by remarking that $\mathcal B\,:\, L^2 (\Omega ) \to L^2(\Omega )$ is self-adjoint, since

\begin{equation*} \int _\Omega v(x) \mathcal B u (x) \mathrm{d} x = \int _\Omega \int _\Omega \frac{u(y) v(x)}{|x-y|^{d-2\left(1-\frac{\alpha }{2}\right)}} \mathrm{d} x \mathrm{d} y = \int _\Omega u(y) \mathcal B v (y) \mathrm{d} y. \end{equation*}

Since $\mathcal I_\alpha$ is a compact operator on $L^2_c ({\mathbb{R}^d})$ , so is $\mathcal B$ in $L^2 (\Omega )$ . Thus, by the spectral theorem, there exists a basis of $L^2(\Omega )$ of orthonormal eigenfunctions $\varphi _i$ of $\mathcal B$ with eigenvalues $\lambda _i \to 0$ , and, defining $u_i = \int _\Omega u \varphi _i \mathrm{d} x$ , it holds

\begin{equation*} u(x) = \sum _{i=1}^\infty u_i \varphi _i(x), \qquad \text{and} \qquad \mathcal B u(x) = \sum _{i=1}^\infty \lambda _i u_i \varphi _i(x). \end{equation*}

Furthermore, with this construction $\|u\|_{L^2} = \sum _i |u_i|^2$ .

Let us show that $\lambda _i \ge 0$ . Defining $U \,:\!=\, \mathcal B u$ , notice that $({-}\Delta )^{1-\frac{\alpha }{2}} U = \mathcal E (u)$ in $\mathbb{R}^d$ . Therefore, for $u \in C_c^\infty (\Omega )$ , we have

\begin{equation*} \int _\Omega u \mathcal B u = \int _\Omega U ({-}\Delta )^{1-\frac{\alpha }{2}} U = \int _{\mathbb{R}^d} U ({-}\Delta )^{1-\frac{\alpha }{2}} U = \int _{\mathbb{R}^d} \left |({-}\Delta )^{ \frac{1-\frac{\alpha }{2}}{2} } U \right |^2 \ge 0. \end{equation*}

Hence, $\lambda _i \ge 0$ .

We now show that $\lambda _i \gt 0$ for all $i$ . Suppose, to the contrary, that $\lambda _i = 0$ for some $i$ . We therefore have that $U \,:\!=\, \mathcal B \varphi _i = 0$ . But then, a.e. in $\Omega$ we have that $\varphi _i = \mathcal E(\varphi _i) = ({-}\Delta )^{1-\frac{\alpha }{2}} U = 0$ , a contradiction.

Therefore, we can formally define the operator $\mathcal B^{-1}$ through the series $\mathcal B^{-1}u (x) = \sum _i \lambda _i^{-1} u_i \varphi _i(x)$ . We define

\begin{equation*} D(\mathcal A) = \Big\{ u = \mathcal B^{-1} v \,:\, v \in H^2 (\Omega ), \nabla v \cdot n = 0 \text{ on } \partial \Omega \text{, and } \sum _i \lambda _i^{-2} |v_i|^2 \lt \infty \Big\}. \end{equation*}

Notice, by construction, that $D(\mathcal A) = L^2 (\Omega ) \cap \mathcal B^{-1} (H^2(\Omega ))$ . Since $\lambda _i \gt 0$ , this set is not empty. Then $\mathcal A \,:\, D(\mathcal A) \subset L^2 (\Omega ) \to L^2 (\Omega )$ .

Now we check that $\mathcal A$ is monotone. Take $u \in D(\mathcal A)$ . Due the spectral decomposition of $\mathcal B$ it admits a square root and inverse square root $\mathcal B^{\pm \frac 1 2}$ ; take $w = \mathcal B^{\frac 1 2} u$ . Then, due to the Neumann boundary condition

\begin{equation*} \int _\Omega u (\mathcal Au) = \int _\Omega \nabla u \cdot \nabla \mathcal Au = \int _\Omega |\nabla \mathcal B^{\frac 1 2} w|^2 \ge 0. \end{equation*}

Lastly, we check that $\mathcal A$ is maximal monotone. Take $f \in L^2(\Omega )$ ; we want to show there exists $u \in D(\mathcal A)$ such that $u + \mathcal A u = f$ . Consider the weak formulation

\begin{equation*} \int _\Omega u \varphi + \int _\Omega \nabla \mathcal B u \cdot \nabla \varphi = \int _\Omega f \varphi, \qquad \forall \varphi \in H^1 (\Omega ). \end{equation*}

Letting again $w = \mathcal B^{\frac 1 2} u$ and $\psi = \mathcal B^{-\frac 1 2} \varphi$ , we obtain

\begin{equation*} \int _\Omega w \psi + \int _\Omega \nabla \mathcal B^{\frac 1 2} w \cdot \nabla \mathcal B^{\frac 1 2} \psi = \int _\Omega \mathcal B^{-\frac 1 2} w \mathcal B^{\frac 1 2} \psi + \int _\Omega \nabla \mathcal B^{\frac 1 2} w \cdot \nabla \mathcal B^{\frac 1 2} \psi = \int _\Omega f \mathcal B^{\frac 1 2} \psi. \end{equation*}

This is a problem of the form $a(w,\psi ) = L(\psi )$ where $w,\psi \in V = \mathcal B^{-\frac 1 2}(H^1 (\Omega )) \cap L^2 (\Omega )$ , where the bilinear form $a$ is symmetric and continuous in $V$ . Hence, it can be solved using the Lax–Milgram theorem. A posteriori, it is trivial to verify that $u \in D(\mathcal A)$ .

We now satisfy all the hypotheses of the Hille–Yosida theorem. Thus, if $u_0 \in D(\mathcal A)$ , then $\mathcal S(t) u_0$ is a solution in the strong sense, i.e. $u \in C([0,\infty ), D(A)) \cap C^1([0,\infty ), L^2(\Omega ))$ , and the equation is satisfied.

The problem (1.6) is not purely diffusive when $\beta \gt 0$ , hence the existence does not follow directly from the Hille–Yosida (or Lumer–Phillips) theorems. Unlike in $\mathbb{R}^d$ , it cannot be deduced from the diffusive problem by a change of variables that would lead to a domain $\Omega _t$ that evolves in time. Thus, the theory of well-posedness for (1.6) when $\beta \gt 0$ is an open problem. A sensible approach would be to prove the convergence of our numerical scheme below.

3. Numerical schemes

The thrust of this work is the discretisation of the fractional Laplacian term in equation (1.6). First, we introduce the scheme in one spatial dimension, in order to highlight the technique used to approximate the fractional term, and then we generalise it to higher dimensions. Our discretisation of the advection term follows previous finite-volume works for generalised Fokker–Planck equations.

3.1. One dimension

In one dimension, we construct a scheme for equation (1.6) in the range $1\lt \alpha \lt 2$ . We consider, without loss of generality, a domain $\Omega = ({-}R,R)$ and divide it into $N$ cells $C_i = \big[{x_{i-{\frac{1}{2}}}, x_{i+{\frac{1}{2}}}}\big]$ , for $i=1,\cdots,N$ . Each cell is centred at the points $x_i$ , where $x_i=-R+(i-1/2)\Delta x$ . For simplicity, we assume a uniform grid with cell size $\Delta x = 2RN^{-1}$ .

We denote by $\bar{\rho }_i(t)$ the average of the solution $\rho (t,x)$ over the $i$ th cell:

\begin{equation*} \bar {\rho }_i(t) = \frac {1}{\Delta x} \int _{C_i} \rho (t, x) \mathrm {d} x. \end{equation*}

Then, equation (1.6) can be integrated on each cell $C_i$ to yield

\begin{equation*}{\frac{\mathrm{d}{{\bar{\rho }}_i(t)}}{\mathrm{d}{t}}} +{\frac{F\big[{\rho \big({t,x_{i+{\frac{1}{2}}}}\big)}\big] - F\big[{\rho \big({t,x_{i-{\frac{1}{2}}}}\big)}\big]}{\Delta x}} = 0, \quad i=1, \cdots, N, \end{equation*}

which we approximate as

(3.1a) \begin{align}{\frac{\mathrm{d}{\bar{\rho }_i(t)}}{\mathrm{d}{t}}} + \frac{F_{i+{\frac{1}{2}}}(t) - F_{i-{\frac{1}{2}}}(t)}{\Delta x} = 0, \quad i=1, \cdots, N. \end{align}

The flux $F$ is split into an advective part $F^{\textrm{ad}}$ and a diffusive part $F^{\textrm{dif}}$ :

(3.1b) \begin{align} F_{i+{\frac{1}{2}}}(t) = F^{\textrm{ad}}_{i+{\frac{1}{2}}}(t) + F^{\textrm{dif}}_{i+{\frac{1}{2}}}(t). \end{align}

The advection flux corresponds to the discretisation of the Fokker–Planck term $\beta \nabla \cdot ({\rho x})$ ; here we follow the discretisation of [5, 15]:

(3.1c) \begin{align} F^{\textrm{ad}}_{i+{\frac{1}{2}}}(t) &= \bar{\rho }_i(t)\big({v_{i+{\frac{1}{2}}}}\big)^{+} + \bar{\rho }_{i+1}(t)\neg{v_{i+{\frac{1}{2}}}}, \end{align}

where

(3.1d) \begin{align} v_{i+{\frac{1}{2}}} = -\frac{\xi _{i+1}-\xi _i}{\Delta x} \text{ and } \xi _i = \beta \frac{|{x_i}|^2}{2}, \end{align}

for $({s})^{+}=\max \{{0,s}\}$ and $\neg{s}=\min \{0,s\}$ .

To treat the diffusion, the gradient term $ \mathcal{C}({1,\alpha -2}) \nabla \int _{\mathbb{R}^d} \rho (y)|{x-y}|^{1-\alpha } \mathrm{d} y$ is replaced by the difference

(3.1e) \begin{align} F^{\textrm{dif}}_{i+{\frac{1}{2}}}(t) = \frac{I(t,x_{i+1}) - I(t,x_i)}{\Delta x}. \end{align}

The term $I(t,x_i)$ is the approximation of the integral $\mathcal I_{2-\alpha }[\rho ](x_i)$ , given as a discrete sum by

(3.1f) \begin{align} I(t,x_i) & = \sum _{k = 1}^{N} \bar{\rho }_{k}(t) I_{k}(x_i), \quad \textrm{where}\quad I_{k}(x_i) \,:\!=\, \mathcal{C}(1, \alpha - 2) \int _{C_{k}} |{x_i - y}|^{1 - \alpha } \mathrm{d} y. \end{align}

To conclude, we impose no-flux boundary conditions:

(3.1g) \begin{align} F_{1-{\frac{1}{2}}}(t) = F_{N+{\frac{1}{2}}}(t) \equiv 0. \end{align}

Remark 3.1 (Linearity). Scheme (3.1) is linear. We can rewrite equation (3.1a) as a system

\begin{equation*}{\frac{\mathrm{d}{\boldsymbol{\bar{\rho }}(t)}}{\mathrm{d}{t}}} + A \boldsymbol{\bar{\rho }}(t) = 0, \end{equation*}

for a vector $\boldsymbol{\bar{\rho }} = \begin{pmatrix} \bar{\rho }_1 & \cdots & \bar{\rho }_i & \cdots & \bar{\rho }_N \end{pmatrix}^\top$ . As with the flux, the constant matrix $A$ can be split into an advective part and a diffusive part, $A=A^{\textrm{ad}}+A^{\textrm{dif}}$ . The advection matrix is simply

(3.2) \begin{align} A^{\textrm{ad}} & = \frac{\beta }{\Delta x} \begin{pmatrix} \big({v_{1+{\frac{1}{2}}}}\big)^{+} & \quad \neg{v_{1+{\frac{1}{2}}}} \\ \ddots & \quad \ddots & \quad \ddots \\ & \quad -\big({v_{i-{\frac{1}{2}}}}\big)^{+} & \quad -\neg{v_{i-{\frac{1}{2}}}}+\big({v_{i+{\frac{1}{2}}}}\big)^{+} & \quad \neg{v_{i+{\frac{1}{2}}}} \\ & \quad & \quad \ddots & \quad \ddots & \quad \ddots \\ & \quad & \quad & \quad -\big({v_{N-{\frac{1}{2}}}}\big)^{+} & \quad -\neg{v_{N-{\frac{1}{2}}}} \end{pmatrix}. \end{align}

The diffusion matrix can be written as the product of a discrete Laplacian and a dense matrix, $A^{\textrm{dif}}=LD$ , where

(3.3) \begin{align} L & = -\frac{1}{\Delta x^2} \begin{pmatrix} -1 & \quad 1 \\[3pt] 1 & \quad -2 & \quad 1 \\[3pt] & \quad \ddots & \quad \ddots & \quad \ddots \\[3pt] & \quad & \quad 1 & \quad -2 & \quad 1 \\[3pt] & \quad & \quad & \quad 1 & \quad -1 \end{pmatrix}, \quad D = \begin{pmatrix} I_1(x_1) & \quad I_2(x_1) & \quad \ldots & \quad I_N(x_1) \\[3pt] I_1(x_2) & \quad I_2(x_2) & \quad \ldots & \quad I_N(x_2) \\[3pt] \vdots & \quad \vdots & \quad \ddots & \quad \vdots \\[3pt] I_1(x_N) & \quad I_2(x_N) & \quad \ldots & \quad I_N(x_N) \end{pmatrix}. \end{align}

Remark 3.2 (Symmetry). The terms in equation (3.1f) are symmetric, $I_{k}(x_i)=I_i(x_{k})$ . The matrix $D$ is thus symmetric, and it can be constructed by evaluating only $N$ terms.

Remark 3.3 (Higher order advection). The discretisation of the fractional diffusion term is consistent to second order, a fact which will be verified in section 4. However, the treatment of the advection term described above is only first-order accurate. A higher order discretisation (e.g. flux limiters, MUSCL) may be used, at the expense of the linearity of the scheme. An example is presented in the Appendix.

Remark 3.4 (Time discretisation). In practice, we discretise scheme (3.1) on the interval $t\in ({0,T})$ using a uniform step $\Delta t$ . For the sake of stability, we employ the implicit time discretisation

\begin{equation*} \boldsymbol{\bar{\rho }}^{m+1} = ({\textrm{Id}+\Delta t A})^{-1}\boldsymbol{\bar{\rho }}^m, \end{equation*}

where $\textrm{Id}$ is the identity matrix. The update matrix $({\textrm{Id}+\Delta t A})^{-1}$ is computed once, offline, for each mesh size $(\Delta t,\Delta x)$ , and then stored for successive use.

Remark 3.5 (The range $\alpha \leq 1$ ). The scheme presented in subsection 3.1 is only valid in the range $1\lt \alpha \lt 2$ due to the inversion formula (1.4) used to rewrite the Lévy–Fokker–Planck equation as (1.5). The range $0\lt \alpha \leq 1$ can be handled using instead the inversion formulae for the Poisson problem on the ball. The relevant kernels are given in [13, Section 3]:

\begin{align*} I_{k}(x_i) & = \frac{1}{\pi } \int _{x_{k-{\frac{1}{2}}}}^{x_{k+{\frac{1}{2}}}} \log \left ({ \frac{ R^2 - x_i y + \sqrt{\left ({R^2-x_i^2}\right )\left ({R^2-y^2}\right )} }{ R|{x_i - y}| } }\right )\mathrm{d} y,\quad \text{for }\alpha =1; \\ I_{k}(x_i) & = \kappa ({1, \alpha }) \int _{x_{k-{\frac{1}{2}}}}^{x_{k+{\frac{1}{2}}}} |{x_i - y}|^{\alpha - 1} \int _0^{r_0({x_i, y})} \frac{t^{\frac{\alpha }{2}-1}}{\big({t+1}^{\frac{1}{2}}\big)} \mathrm{d} t\,\mathrm{d} y,\quad \text{for } 0 \lt \alpha \lt 1; \end{align*}

where

\begin{align*} r_0({x_i,y}) = \frac{ \left ({R^2-|{x_i}|^2}\right )\left ({R^2-|{y}|^2}\right ) }{ R^2|{x_i - y}|^2 }, \quad \kappa ({1, \alpha }) = \frac{\Gamma (\frac{1}{2})}{2^{\alpha } \pi ^{\frac{1}{2}} \Gamma ^{2}\left(\frac{\alpha }{2}\right)}. \end{align*}

Unfortunately, this approach renders the matrix $D$ no longer symmetric. We will not address this case directly.

3.2. Two dimensions

In two dimensions, the inversion formula (1.4) no longer restricts the fractional exponent. Therefore, we construct a scheme for equation (1.6) in the range $0\lt \alpha \lt 2$ .

We consider without loss of generality a square domain $\Omega = ({-}R,R)^2$ , divided into $N^2$ cells given by $C_{i,\,j} = [{x_{i-{\frac{1}{2}}}, x_{i+{\frac{1}{2}}}}] \times [{y_{j-{\frac{1}{2}}}, y_{j+{\frac{1}{2}}}}]$ , for $i,j=1,\cdots,N$ . Each cell is centred at the points $({x_i,y_{j}})$ , where $x_i=-R+(i-1/2)\Delta x$ , $y_j=-R+(j-1/2)\Delta y$ and $\Delta x = \Delta y = 2RN^{-1}$ . As in subsection 3.1, we approximate the cell averages of equation (1.6) to arrive at

(3.4a) \begin{align} {\frac{\mathrm{d}{\bar{\rho }_{i,\,j}(t)}}{\mathrm{d}{t}}}+ \frac{F_{i+{\frac{1}{2}},\,j}(t) - F_{i-{\frac{1}{2}},\,j}(t)}{\Delta x} + \frac{G_{i,\,j+{\frac{1}{2}}}(t) - G_{i,\,j-{\frac{1}{2}}}(t)}{\Delta y} = 0, \quad i, j = 1, \cdots, N. \end{align}

Once again, the fluxes are split into an advective part $F^{\textrm{ad}}$ and a diffusive part $F^{\textrm{dif}}$ :

(3.4b) \begin{align} F_{i+{\frac{1}{2}},\,j}(t) = F^{\textrm{ad}}_{i+{\frac{1}{2}},\,j}(t) + F^{\textrm{dif}}_{i+{\frac{1}{2}},\,j}(t), \quad G_{i,\,j+{\frac{1}{2}}}(t) = G^{\textrm{ad}}_{i,\,j+{\frac{1}{2}}}(t) + G^{\textrm{dif}}_{i,\,j+{\frac{1}{2}}}(t). \end{align}

The advection terms are now

(3.4c) \begin{align} F^{\textrm{ad}}_{i+{\frac{1}{2}},\,j}(t) & = \bar{\rho }_{i,\,j}(t)\big({v_{i+{\frac{1}{2}},\,j}}\big)^{+} + \bar{\rho }_{i+1,\,j}(t)\neg{v_{i+{\frac{1}{2}},\,j}}, \end{align}

(3.4d) \begin{align} G^{\textrm{ad}}_{i,\,j+{\frac{1}{2}}}(t) & = \bar{\rho }_{i,\,j}(t)\big({w_{i,\,j+{\frac{1}{2}}}}\big)^{+} + \bar{\rho }_{i,\,j+1}(t)\neg{w_{i,\,j+{\frac{1}{2}}}}, \end{align}

where

(3.4e) \begin{align} v_{i+{\frac{1}{2}},\,j} = -\frac{\xi _{i+1,\,j}-\xi _{i,\,j}}{\Delta x},\quad w_{i,\,j+{\frac{1}{2}}} = -\frac{\xi _{i,\,j+1}-\xi _{i,\,j}}{\Delta y},\quad \xi _{i,\,j} = \beta \frac{|{x_i}|^2+|{y_j}|^2}{2}, \end{align}

following [5, 15].

The treatment of the diffusive part described in the previous section generalises to two dimensions:

(3.4f) \begin{align} F_{i+{\frac{1}{2}},\,j}^{\textrm{dif}}(t) = \frac{I(t, x_{i+1}, y_j) - I(t, x_i, y_j)}{\Delta x} \quad \textrm{and} \quad G_{i,\,j+{\frac{1}{2}}}^{\textrm{dif}}(t) = \frac{I(t, x_i, y_jp) - I(t, x_i, y_j)}{\Delta y}. \end{align}

The discrete integrals $I(t, x_i, y_j)$ are given by the sum

(3.4g) \begin{align} I(t, x_i, y_j) = \sum _{k, l = 1}^{N} \bar{\rho }_{k,\,l}(t) I_{k,\,l}(x_i, y_j), \end{align}

where

(3.4h) \begin{align} I_{k,\,l}(x_i, y_j) \,:\!=\, \mathcal{C}(2, \alpha - 2) \int _{C_{k,\,l}} |{({x_i, y_j}) - ({u,v})}|^{-\alpha }{\mathrm{d} u}\,\mathrm{d} v. \end{align}

Once again, we impose no-flux boundary conditions:

(3.4i) \begin{align} F_{1-{\frac{1}{2}},\,j}(t) = F_{N+{\frac{1}{2}},\,j}(t) \equiv 0, \quad G_{i,\,1-{\frac{1}{2}}}(t) = G_{i,\,N+{\frac{1}{2}}}(t) \equiv 0, \quad i, j = 1, \cdots, N. \end{align}

3.2.1. Dimensional splitting

As in the one-dimensional case, scheme (3.4) is linear. Writing

(3.5) \begin{equation} {\frac{\mathrm{d}{\boldsymbol{\bar{\rho }}(t)}}{\mathrm{d}(t)}} + A \boldsymbol{\bar{\rho }}(t) = 0 \end{equation}

for a vector $\boldsymbol{\bar{\rho }} = \begin{pmatrix} \bar{\rho }_{1,1} & \bar{\rho }_{2,1} & \cdots & \bar{\rho }_{N,1} & \bar{\rho }_{1,2} & \cdots & \bar{\rho }_{i,j} & \cdots & \bar{\rho }_{N,N} \end{pmatrix}^\top$ , we may split the matrix of the scheme into an advective part and a diffusive part, $A=A^{\textrm{ad}}+A^{\textrm{dif}}$ , which are two-dimensional generalisations of (3.2) and (3.3).

However, the linear treatment of the scheme becomes impractical here, as the storage required for matrix of the scheme grows exponentially. While $A^{\textrm{ad}}$ is a banded matrix, $A^{\textrm{dif}}$ is dense. For $N=2^7$ , a dense $N^2\times N^2$ matrix would require 2 gigabytes of RAM. For $N=2^8$ , the size would be 16 gigabytes, or 128 gigabytes for $N=2^9$ . In order to handle the computation, an approach which does not require the direct inversion of the matrix $({\textrm{Id}+\Delta t A})$ is required. The Krylov subspace methods, such as GMRES or BiCGSTAB [36], are among the available options. Instead, we resort to a dimensional splitting strategy.

The (matrix) operator $A$ in equation (3.5) is decomposed as $A = A_2 + A_1$ , where $A_1$ corresponds to the transport terms along the $x$ -direction (those in equation (3.4a) which arise from the fluxes $F_{i+{\frac{1}{2}},\,j}$ ), and $A_2$ corresponds to the transport along the $y$ -direction (terms related to $G_{i,\,j+{\frac{1}{2}}}$ ); naturally, $A_1$ and $A_2$ are independent of $\boldsymbol{\bar{\rho }}$ , as is $A$ . Formally, the solution to (3.5) can be written as $\boldsymbol{\bar{\rho }}(t) = \exp (tA) \boldsymbol{\bar{\rho }}_0$ . One would like to approximate this by $\exp (tA_2)\exp (tA_1) \boldsymbol{\bar{\rho }}_0$ (i.e. by solving the problem one dimension at a time) but, in general, the solution operator cannot be factored in that way: $\exp (tA_2+tA_1) \neq \exp (tA_2)\exp (tA_1)$ . However, the Lie–Trotter (or Trotter–Kato) formula

\begin{align*} \exp ({tA_2+tA_1}) = \lim \limits _{n\rightarrow \infty }{\left ({\exp \left ({\frac{t}{n}A_2}\right )\exp \left ({\frac{t}{n}A_1}\right )}\right )}^n \end{align*}

does hold for general square matrices [39] and some linear operators [27] and has been used to study the convergence of dimensional splitting in the case of linear semi-groups [26]. Choosing $\Delta t = t n^{-1}$ , we see that the exact solution operator $\exp ({tA_2+tA_1})$ can be approximated by applying the operators $\exp ({\Delta t A_1})$ and $\exp ({\Delta t A_2})$ in an alternating sequence, i.e. the exact solution $\boldsymbol{\bar{\rho }}(t)$ can be approximated by performing a sequence of intermediate updates of $\boldsymbol{\bar{\rho }}_0$ , each involving a short time, alternating the $x$ -direction and $y$ -direction sub-problems. Upon discretising time, the approximate solution $\boldsymbol{\bar{\rho }}^{m+1}$ at time $(m+1)\Delta t$ is computed from $\boldsymbol{\bar{\rho }}^m$ (that at time $m\Delta t$ ) via $\boldsymbol{\bar{\rho }}^{m+{\frac{1}{2}}}$ , an intermediate step; $\boldsymbol{\bar{\rho }}^{m+{\frac{1}{2}}}$ is computed from $\boldsymbol{\bar{\rho }}^m$ by solving the $x$ -direction problem and $\boldsymbol{\bar{\rho }}^{m+1}$ is found from $\boldsymbol{\bar{\rho }}^{m+{\frac{1}{2}}}$ by solving the $y$ -problem.

At this stage, the advantage of the dimensional splitting approach is not clear: the matrices $A_1$ and $A_2$ are dense, as was $A$ , so the memory requirement has effectively doubled. However, one further approximation is possible: each of the dimensional updates can be approximately decomposed row-wise or column-wise. For instance, to compute $\boldsymbol{\bar{\rho }}^{m+{\frac{1}{2}}}$ from $\boldsymbol{\bar{\rho }}^m$ : for each row $j$ , compute $\bar{\rho }^{m+{\frac{1}{2}}}_{i,\,j}$ by solving the one-dimensional implicit problem within the row, assuming the value of the density will not change outside of it (i.e. $\bar{\rho }^{m+{\frac{1}{2}}}_{i,\,k} \equiv \bar{\rho }^{m}_{i,\,k}$ whenever $k\neq j$ ). This is done independently on each row, and therefore can be trivially parallelised. A schematic diagram of the update is shown in Figure 1. Each update now involves the inversion of an $N\times N$ matrix, rather than $N^2\times N^2$ , though the matrices are no longer independent of $\boldsymbol{\bar{\rho }}^m$ . To obtain $\boldsymbol{\bar{\rho }}^{m+1}$ , the process is repeated along the $y$ -direction, mutatis mutandis.

While this approach to dimensional splitting is partially justified by Lie–Trotter formula above, we will nevertheless justify it numerically in section 4, both in terms of checking the convergence of the scheme and its long-time behaviour.

Figure 1. Dimensional splitting, row update. The split implicit problem considers information on the whole domain, but the density is allowed to change only within a single row. These updates take place independently for each row in parallel and can be parallelised.

Remark 3.6 (Sweeping dimensional splitting). A valid alternative is the sweeping dimensional splitting described in [5]. In that approach, the row and column updates take place sequentially, each considering the updated information from the previous step. This approach can be beneficial in some settings (it was used in [5] to prove structural properties of the scheme), but was discarded here because it cannot be parallelised.

4. Numerical experiments

We now demonstrate the accuracy and performance of our scheme in a variety of test cases, both in one and two dimensions. We will refer to the fractional heat equation (1.1) and the Lévy–Fokker–Plank equation (1.2) in the discussion; however, for the numerics, these are always understood as equation (1.6) with $\beta =0$ and $\beta =1$ , respectively.

In one dimension, we employ scheme (3.1); in two dimensions, we employ scheme (3.4) with the dimensional splitting described in subsection 3.2.1. Experiments use the first-order upwind fluxes (3.1c) and (3.4c), unless otherwise stated. The experiments that compute the order of accuracy of the scheme use instead the second-order minmod flux (A.1) presented in the Appendix and discussed in Remark 3.3.

4.1. One dimension

4.1.1. Fractional diffusion

We first consider the fractional heat equation (1.1). As in the classical heat equation, an explicit self-similar solution on the whole space is known when $\alpha =1$ :

\begin{equation*} \phi{(t,x)} = C(d){\frac{t}{\left ({t^{2}+|x|^{2}}\right )}}^{\frac{d+1}{2}}. \end{equation*}

Notice that the problem is linear. We pick $C(d)$ so that $\| \phi \|_{L^1} = 1$ , i.e., $C(1) = \frac 1 \pi$ and $C(2) = \frac 1{2\pi }$ . We shall use this explicit solution to validate our numerical scheme. Technically, scheme (3.1) is not valid for $\alpha = 1$ ; however, we can set $\alpha =1+\varepsilon$ and perform the comparison regardless. In practice, we choose $\varepsilon = 10^{-11}$ .

Figure 2 shows a comparison of the numerical solution ( $\alpha =1+\varepsilon$ , $R=100$ , $\Delta x=0.1$ , $\Delta t=0.1$ ) on $\Omega =({-}R,R)$ and the restriction of $\phi$ to $\Omega$ . The initial datum is taken as $\phi (\Delta t, x)$ . Both solutions match well on the interior of the domain; however, there is a clear discrepancy on the boundary, where the numerical solution behaves singularly. The discrepancy is explained by the fact that the self-similar profile $\phi$ is leptokurtic (i.e. has higher kurtosis, or thicker tails, than a Gaussian); therefore, the amount of mass that is ignored by considering $\phi$ on a bounded domain is never exponentially small. The singular behaviour at the boundary is a known effect of certain fractional operators [2]. This effect is explored further in the next experiment.

Figure 2. Fractional heat equation (1.1) in one dimension. Numerical solution $\boldsymbol{\bar{\rho }}^m$ on $\Omega =({-}R,R)$ and explicit solution $\phi$ on $\mathbb{R}$ . Scheme (3.1), $\alpha = 1 + \varepsilon$ , $R=100$ , $\Delta x=0.1$ , $\Delta t=0.1$ . Good agreement is shown on the interior of the domain; boundary effects are visible.

4.1.2. Singular behaviour at the boundary

We consider here the steady states of the fractional heat equation (1.1) on a bounded domain in order to explore the singular behaviour at the boundary. In one dimension, the steady state $\rho _\infty$ of (1.6) with $\beta =0$ satisfies

\begin{align*} \partial _{xx}\left ({\int _{-R}^{R} \frac{\rho _\infty (y)}{|x-y|^{\alpha -1}} \mathrm{d} y}\right ) = 0, \end{align*}

which, upon considering the boundary conditions, reduces to

\begin{equation*} \int _{-R}^{R} \frac{\rho _\infty (y)}{|x-y|^{\alpha -1}} \mathrm{d} y = C \end{equation*}

for some constant $C$ . For $\alpha \gt 1$ , the steady profile can be found explicitly:

(4.1) \begin{align} \rho _\infty (x) & = \cos ^{2} \left ({\frac{\pi (\alpha - 1)}{2}}\right ) \frac{ C(R+x)^{\frac{\alpha }{2} - 1} }{ \pi ^{2}(R-x)^{1 - \frac{\alpha }{2}} } \int _{-R}^{R} \frac{ (R-y)^{1-\frac{\alpha }{2}} (R+y)^{\frac{\alpha }{2} - 1} }{ y-x } \mathrm{d} y + \frac{ A\sin ({\pi (\alpha - 1)}) }{ 2 \pi (R+x)^{2 - \alpha } }; \end{align}

see [23] for details.

Figure 3 shows the numerical solution ( $\alpha =1.5$ , $R=50$ , $\Delta x=0.1$ , $\Delta t=0.5$ ) on $\Omega =({-}R,R)$ as it tends to the stationary profile (4.1). The datum is taken as in the previous section. The explicit steady state is captured by the numerical solution as time grows. Note, however, that we have to run the simulation for a long time before the match is apparent; this is in contrast to the experiment in the next section. The slow convergence may be due to the singular behaviour at the boundary.

Figure 3. Fractional heat equation (1.1) in one dimension. Numerical solution $\boldsymbol{\bar{\rho }}^m$ and explicit steady state $\rho _\infty$ on $\Omega =({-}R,R)$ . Scheme (3.1), $\alpha = 1.5$ , $R=50$ , $\Delta x=0.1$ , $\Delta t=0.5$ . The numerical solution tends to $\rho _\infty$ .

4.1.3. Steady states as a function of domain size

We now turn to the Lévy–Fokker–Planck equation (1.2). First, we consider the case $\alpha =1$ , where an explicit solution on the whole line is known:

\begin{equation*} \rho ^{\ast }(t,x) = \frac{1}{\pi } \frac{ e^{t}(e^{t}-1) }{ (1+x^{2})e^{2t} - 2e^{t} + 1 }; \end{equation*}

as $t\rightarrow \infty$ , this solution tends to the steady state

(4.2) \begin{equation} \rho _{\infty }(x) = \frac{1}{\pi } \frac{1}{1+x^{2}}. \end{equation}

Figure 4 shows the numerical solution ( $\alpha =1+\varepsilon$ , $R=50$ , $\Delta x=0.05$ , $\Delta t=0.01$ ) on $\Omega =({-}R,R)$ compared to the explicit steady state (4.3). The datum for the numerical solution is a uniform distribution with unit mass. Once again, the explicit steady state is captured well by the numerical solution as time grows. Unlike in the previous experiment, this solution approaches the corresponding steady state very rapidly.

Figure 4. Lévy–Fokker–Planck equation (1.2) in one dimension. Numerical solution $\boldsymbol{\bar{\rho }}^m$ , exact solution $\rho ^*$ , and explicit steady state $\rho _\infty$ . Scheme (3.1), $\alpha = 1+\varepsilon$ , $R=50$ , $\Delta x=0.05$ , $\Delta t=0.01$ . The numerical solution clearly tends to $\rho _\infty$ .

As was the case with the fractional heat equation, the typical solution of the Lévy–Fokker–Planck equation is leptokurtic, as it has algebraic tails. Thus, the error committed when a whole-space solution is restricted to a bounded domain is not exponentially small, even if the presence of the Fokker–Planck term prevents singularities from developing at the boundary. We therefore expect that the steady state in a bounded domain will differ from (4.3) by a non-trivial amount.

Figure 5 shows the $L^{1}({\Omega })$ distance between the numerical steady state of the Lévy–Fokker–Planck equation ( $\alpha =1+\varepsilon$ , $\Delta x=2R/2^{12}$ , $\Delta t=0.1$ ) on $\Omega =({-R,R})$ for various values of $R$ , and the explicit steady state (4.3). As expected, the error decreases as $R$ tends to infinity, though the decay does not follow an obvious pattern.

We show the relative entropy of our numerical solution (for $\alpha = 1+\varepsilon$ ) with respect to the equilibrium $\rho _{\infty }$ in Figure 6. The results show good agreement with the exponential trend predicted by [24], using two most common entropy functions: $\Phi (x) = (x -1)^{2}$ and $\Phi (x) = x(\log x - 1) + 1$ .

Figure 5. Lévy–Fokker–Planck equation (1.2) in one dimension. $L^{1}({\Omega })$ distance between the numerical steady state with $\alpha =1+\varepsilon$ on $\Omega =({-R,R})$ and the explicit steady state with $\alpha =1$ . Scheme (3.1), $\Delta x=2R/2^{12}$ , $\Delta t=0.1$ . The mismatch decreases as $R$ increases.

Figure 6. Lévy–Fokker–Planck equation (1.2) in one dimension. Dissipation of the relative entropy with $\alpha =1 + \varepsilon$ with respect to the equilibrium $\rho _{\infty }$ . Scheme (3.1), $R=50$ , $\Delta x=0.05$ , $\Delta t=0.01$ . Black reference line has slope one.

A similar analysis can be performed when $\alpha/2=1$ . The Lévy–Fokker–Planck equation reduces to the classical Fokker–Planck equation

\begin{align*} \partial _t\rho = \Delta \rho + \nabla \cdot ({x\rho }), \end{align*}

whose unique, asymptotically stable steady state is

(4.3) \begin{align} \rho ^{\infty }(x) = \frac{1}{(2 \pi )^{\frac{d}{2}}} e^{-\frac{|x|^{2}}{2}} \end{align}

in dimension $d$ (for solutions with unit mass on the whole space). If we let $\frac{\alpha }{2} = 0.99$ , the steady state of our numerical scheme should be close to this one, and the agreement should improve as the domain grows.

Figure 7 shows the $L^{1}({\Omega })$ distance between the numerical steady state of the Lévy–Fokker–Planck equation ( $\alpha/2=0.99$ , $\Delta x=2R/2^{12}$ , $\Delta t=0.1$ ) on $\Omega =({-R,R})$ for various values of $R$ and the explicit steady state (4.4). Once again, the error decreases as $R$ tends to infinity, as expected.

Figure 8 shows the numerical steady states of the Lévy–Fokker–Planck equation (1.2) ( $R=50$ , $\Delta x = 0.05$ , $\Delta t = 0.01$ ) on $\Omega =({-}R,R)$ for various fractional orders $\alpha \in (1,2)$ . We recover symmetric distributions with algebraic tails that become thicker as $\alpha$ decreases. We compare the tails of our numerical results with the expected behaviour predicted in ref. [9], given by $\rho (t,x) \asymp \min \{ 1, |x|^{-\alpha -d} \}$ .

Figure 7. Lévy–Fokker–Planck equation (1.2) in one dimension. $L^{1}({\Omega })$ distance between the numerical steady state with $\alpha/2=0.99$ on $\Omega =({-R,R})$ and the explicit steady state with $\alpha =2$ . Scheme (3.1), $\Delta x=2R/2^{12}$ , $\Delta t=0.1$ . The distance decreases as $R$ increases.

Figure 8. Lévy–Fokker–Planck equation (1.2) in one dimension. Numerical steady state for varying fractional order $\alpha \in (1,2)$ . Scheme (3.1), $R=50$ , $\Delta x=0.05$ , $\Delta t=0.01$ . Top: profile at the centre of the domain. Bottom: detail of the algebraic tails, compared with the predicted asymptotic behaviour $|x|^{-\alpha -d}$ (dashed).

4.1.4. Convergence of steady states

To conclude, we verify the order of convergence of the scheme. We fix the domain and compute the steady state of the Lévy–Fokker–Planck equation (1.2) as in the previous section, for various values of $\alpha$ . We compute the steady states on a sequence of refining meshes, and study their convergence. Since the analytical steady state is not known explicitly, we shall monitor the error between numerical steady states and show that this decays with the mesh size.

Figure 9 shows the $L^{1}({\Omega })$ and $L^2({\Omega })$ distance between successive numerical steady states ( $R=50$ , $\Delta t=\Delta x=2R/ 2^{n}$ for $5\leq n\leq 10$ ) computed with scheme (3.1) on $\Omega =({-}R,R)$ . The scheme is first-order accurate for all fractional orders $\alpha \in ({1,2})$ .

Figure 9. Lévy–Fokker–Planck equation (1.2) in one dimension. Convergence of scheme (3.1) for varying fractional order $\alpha$ . $R=50$ , $\Delta x=2R/2^{n}$ for $5\leq n\leq 10$ , $\Delta t = \Delta x$ . Black reference line has slope one.

Figure 10 performs the same analysis on scheme (3.1) with the second-order flux (A.1) (viz. Remark 3.3), letting $\Delta t = \Delta x^2$ instead. The scheme is second-order accurate for all fractional orders $\alpha \in ({1,2})$ .

Figure 10. Lévy–Fokker–Planck equation (1.2) in one dimension. Convergence of scheme (3.1) with second-order flux (A.1) (viz. Remark 3.3) for varying fractional order $\alpha$ . $R=50$ , $\Delta x=2R/2^{n}$ for $5\leq n\leq 10$ , $\Delta t=\Delta x$ . Black reference line has slope two.

4.2. Two dimensions

4.2.1. Steady states as a function of domain size

We begin our two-dimensional experiments by verifying the behaviour of the dimensionally split scheme. Figure 11 shows the numerical steady states of the Lévy–Fokker–Planck equation (1.2) ( $R=20$ , $\Delta x = \Delta y = 0.15$ , $\Delta t = 0.2$ ) on $\Omega =({-}R,R)^2$ for various fractional orders. We recover radially symmetric distributions with algebraic tails that become thicker as $\alpha$ decreases. We compare the tails of our numerical results with the expected behaviour predicted in ref. [9], which is given by $\rho (t,x) \asymp \min \{ 1, |x|^{-\alpha -d} \}$ .

Figure 11. Lévy–Fokker–Planck equation (1.2) in two dimensions. Numerical steady state for varying fractional order $\alpha \in (0,2)$ . Scheme (3.4) with splitting (viz. subsection 3.2.1), $R=20$ , $\Delta x=\Delta y=0.15$ , $\Delta t=0.2$ . Top: central section. Bottom: detail of the algebraic tails, compared with the dotted lines of the predicted long time asymptotic behaviour $|x|^{-\alpha -d}$ .

As in the one-dimensional case, an explicit solution to the Lévy–Fokker–Planck equation on the whole space is known for $\alpha =1$ . The solution is found from the self-similar solution to the fractional heat equation [19] through the change of variables proposed in ref. [8], just as a solution to the classical Fokker–Planck equation can be derived from a solution to the heat equation. The solution in question is given by

\begin{equation*} \rho ^{\ast }(t, x, y) = \frac{1}{2\pi } \frac{ e^{2t}(e^{t}-1) }{{\left ({(1+x^2+y^2)e^{2t} - 2e^{t} + 1}\right )}^{\frac{3}{2}} }, \end{equation*}

which tends to the steady state

(4.4) \begin{equation} \rho _{\infty }(x, y) = \frac{1}{2\pi } \frac{1}{({1+x^2+y^2})^{\frac{3}{2}}}. \end{equation}

Figure 12 shows the numerical solution ( $\alpha =1$ , $R=20$ , $\Delta x=\Delta y=0.08$ , $\Delta t=0.1$ ) on $\Omega =({-}R,R)^2$ compared to the explicit steady state (4.4). The datum for the numerical solution is a uniform distribution with unit mass.

Figure 12. Lévy–Fokker–Planck equation (1.2) in two dimensions. Numerical steady state and explicit steady state. Scheme (3.4) with splitting (viz. subsection 3.2.1), $\alpha = 1+\varepsilon$ , $R=20$ , $\Delta x=\Delta y=0.08$ , $\Delta t=0.1$ .

We now study the convergence of the numerical steady state to the profile (4.4) as the size of the domain grows. Unlike the one-dimensional test, the two-dimensional analysis can be performed setting $\alpha =1$ exactly. Figure 13 shows the $L^{1}({\Omega })$ distance between the numerical steady state of the Lévy–Fokker–Planck equation ( $\alpha =1$ , $\Delta x=\Delta y=2R/2^{8}$ , $\Delta t=0.1$ ) on $\Omega =({-R,R})^2$ , for various values of $R$ , and the explicit steady state (4.4). As in the one-dimensional case, the error decreases as $R$ tends to infinity.

Figure 13. Lévy–Fokker–Planck equation (1.2) in two dimensions. $L^{1}({\Omega })$ distance between the numerical steady state on $\Omega =({-R,R})^2$ and the explicit steady state. Scheme (3.4) with splitting (viz. subsection 3.2.1), $\alpha =1$ , $\Delta x=2R/2^{8}$ , $\Delta t=0.1$ . The mismatch decreases as $R$ increases.

4.2.2. Convergence of steady states

We verify the order of convergence of the dimensionally split scheme. As in one dimension, we fix the domain size and compute the steady state of the Lévy–Fokker–Planck equation (1.2) for various values of $\alpha$ . Figure 14 shows the $L^{1}({\Omega })$ and $L^2({\Omega })$ distance between numerical steady states ( $R=20$ , $\Delta t=\Delta x=\Delta y=2R/ 2^{n}$ for $5\leq n\leq 8$ ) computed with scheme (3.4) on $\Omega =({-}R,R)^2$ as the mesh size is halved. The order of the scheme appears slightly less than one; this might be a consequence of the dimensional splitting. Noticeably, the convergence is initially very slow when the fractional order is close to zero.

Figure 14. Lévy–Fokker–Planck equation (1.2) in two dimensions. Convergence of scheme (3.4) with splitting (viz. subsection 3.2.1) for varying fractional orders $\alpha$ . $R=20$ , $\Delta x=\Delta y=2R/2^{n}$ for $n=5, \dots, 8$ , $\Delta t=\Delta x$ . Black reference line has slope one.

4.2.3. Long-time asymptotics

To conclude, we study the rate of convergence of the numerical solution of the Lévy–Fokker–Planck equation (1.2) to the corresponding steady states. Figure 15 shows the $L^{1}(\Omega )$ and $L^2(\Omega )$ distances of the numerical solution ( $R=20$ , $\Delta x = \Delta y = 0.15$ , $\Delta t = 0.08$ ) on $\Omega =({-}R,R)^2$ for various fractional orders to their asymptotic steady states as a function of time. Perhaps due to the highly symmetric initial data, the numerical solutions show an improved rate of convergence ( $e^{-2t}$ ) towards the steady state with respect to the result of [24] ( $e^{-\alpha t}$ ). This acceleration phenomena due to symmetry of the datum is well-documented, as it has been observed also in the classical Fokker–Planck setting [7], as well as in the porous medium equation [14].

Figure 15. Lévy–Fokker–Planck equation (1.2) in two dimensions. Convergence to steady state of numerical solutions for varying fractional order. Scheme (3.4) with splitting (viz. subsection 3.2.1), $R=20$ , $\Delta x = \Delta y = 0.15$ , $\Delta t = 0.08$ . Black reference line has slope two.