A semidiscrete finite element approximation of a time-fractional Fokker-Planck equation with nonsmooth initial data

TL;DR

This paper develops a new stability and convergence analysis for a finite element method applied to a time-fractional Fokker-Planck equation with nonsmooth initial data, supported by numerical experiments.

Contribution

It introduces a novel energy-based analysis for spatial discretization of the equation with low-regularity initial data, enhancing understanding of numerical stability and convergence.

Findings

The analysis confirms stability and convergence for the spatial discretization.

Numerical experiments align with theoretical convergence predictions.

The method effectively handles nonsmooth initial data.

Abstract

We present a new stability and convergence analysis for the spatial discretization of a time-fractional Fokker--Planck equation in a convex polyhedral domain, using continuous, piecewise-linear, finite elements. The forcing may depend on time as well as on the spatial variables, and the initial data may have low regularity. Our analysis uses a novel sequence of energy arguments in combination with a generalized Gronwall inequality. Although this theory covers only the spatial discretization, we present numerical experiments with a fully discrete scheme employing a very small time step, and observe results consistent with the predicted convergence behavior.

Tables

Table 1. Table 1: Weighted errors ( 35 ) and convergence rates ( 36 ) for different α 𝛼 \alpha , when u 0 ​ h = P h ​ u 0 subscript 𝑢 0 ℎ subscript 𝑃 ℎ subscript 𝑢 0 u_{0h}=P_{h}u_{0} .

$Q_h$	$\alpha =0.25$		$\alpha =0.50$		$\alpha =0.75$
7	7.98e-03		7.77e-03		7.84e-03
15	1.96e-03	2.024	1.91e-03	2.024	1.94e-03	2.017
31	4.88e-04	2.008	4.75e-04	2.008	4.82e-04	2.007
63	1.21e-04	2.014	1.18e-04	2.014	1.19e-04	2.015

Table 2. Table 2: Weighted errors ( 35 ) and convergence rates ( 36 ) for different α 𝛼 \alpha , when u 0 ​ h = I h ​ u 0 subscript 𝑢 0 ℎ subscript 𝐼 ℎ subscript 𝑢 0 u_{0h}=I_{h}u_{0} .

$Q_h$	$\alpha =0.25$		$\alpha =0.50$		$\alpha =0.75$
7	7.79e-02		7.46e-02		7.27e-02
15	4.04e-02	0.948	3.86e-02	0.950	3.76e-02	0.952
31	2.06e-02	0.973	1.97e-02	0.973	1.91e-02	0.974
63	1.04e-02	0.987	9.93e-03	0.987	9.65e-03	0.987

Equations

$$\partial_{t}u-\nabla\cdot\bigl{(}\partial_{t}^{1-\alpha}\kappa_{\alpha}\nabla u-\vec{F}\partial_{t}^{1-\alpha}u\bigr{)}={0}\quad\text{for $\vec{x}\in\Omega$ and $0<t<T$,}$$

$${\mathcal{I}}^{\alpha }u(t)={\omega }_{\alpha }\ast u(t)={\int }_0^t{\omega }_{\alpha }(t-s)u(s)ds,{\omega }_{\alpha }(t)=\frac{t^{\alpha -1}}{\Gamma (\alpha )}.$$

$$u(\vec{x},t)=0\text{for x∈∂Ω and 0<t<T,}$$

$$\mathcal{I}^{1-\alpha}(\partial_{t}u)-\nabla\cdot\bigl{(}\kappa_{\alpha}\nabla u-\vec{F}(\vec{x})u\bigr{)}={0},$$

$$⟨{\partial }_{t}u,v⟩+⟨{\partial }_t^{1-\alpha }\nabla u,\nabla v⟩-⟨\vec{F}{\partial }_t^{1-\alpha }u,\nabla v⟩=0$$

$$⟨{\partial }_t{u}_h,\chi ⟩+⟨{\partial }_t^{1-\alpha }\nabla u_h,\nabla \chi ⟩-⟨\vec{F}{\partial }_t^{1-\alpha }{u}_h,\nabla \chi ⟩=0\text{for all χ∈Sh,}$$

$$\|u_{h}(t)-u(t)\|\leq Ch^{2}\biggl{(}\|u_{0}\|_{2}^{2}+\int_{0}^{t}\|u^{\prime}(s)\|_{2}^{2}\,ds\biggr{)}^{1/2}.$$

$$\|v\|_{r}=\|(-\nabla^{2})^{r/2}v\|=\biggl{(}\sum_{m=1}^{\infty}\lambda_{m}^{r}\langle v,\varphi_{m}\rangle^{2}\biggr{)}^{1/2}\quad\text{for $r\geq 0$,}$$

$$-{\nabla }^2{\phi }_m={\lambda }_m{\phi }_m\text{in Ω,}\text{with φm=0 on ∂Ω.}$$

$$\Vert u_h(t)\Vert \le C\Vert u_{0h}\Vert \text{for 0≤t≤T.}$$

$$\Vert u(t){\Vert }_2+t\Vert u^{\prime }(t){\Vert }_2\le t^{\alpha (r-2)/2}{K}_r\text{for 0<t≤T.}$$

$$\Vert u_h(t)-u(t)\Vert \le C{t}^{\alpha (r-2)/2}{h}^2{K}_r\text{for 0≤t≤T and 0≤r≤2.}$$

$${\partial }_t^{1-\alpha }{\kappa }_{\alpha }\frac{\partial u}{\partial n}-(\vec{F}\cdot \vec{n}){\partial }_t^{1-\alpha }u=0\text{for x∈∂Ω and 0<t<T,}$$

$${\int }_0^t\Vert {\mathcal{I}}^{\nu }\varphi {\Vert }^{2}ds\le C{t}^{2(\nu -\mu )}{\int }_0^t\Vert {\mathcal{I}}^{\mu }\varphi {\Vert }^{2}ds.$$

$${\int }_0^T\Vert {\mathcal{I}}^{\alpha }\psi (t){\Vert }^{2}dt\le {\omega }_{\alpha +1}(T){\int }_0^T{\omega }_{\alpha }(T-t){\int }_0^t\Vert \psi (s){\Vert }^{2}dsdt,$$

$$\biggl{|}\int_{0}^{t}\langle\phi,\mathcal{I}^{\alpha}\psi\rangle\,ds\biggr{|}\leq\frac{1}{4\epsilon(1-\alpha)^{2}}\int_{0}^{t}\langle\mathcal{I}^{\alpha}\phi,\phi\rangle\,ds+\epsilon\int_{0}^{t}\langle\mathcal{I}^{\alpha}\psi,\psi\rangle\,ds,$$

$${\int }_0^t\Vert {\mathcal{I}}^{\alpha }\varphi {\Vert }^{2}ds\le \frac{C{t}^{\alpha }}{1-\alpha }{\int }_0^t⟨{\mathcal{I}}^{\alpha }\varphi ,\varphi ⟩ds,$$

$${\int }_0^t⟨\varphi ,{\mathcal{I}}^{\alpha }\varphi ⟩ds\le C{t}^{\alpha }{\int }_0^t\Vert \varphi {\Vert }^{2}ds.$$

$$\biggl{|}\int_{0}^{t}\langle\phi,\mathcal{I}^{\alpha}\psi\rangle\,ds\biggr{|}\leq\frac{1}{\cos(\alpha\pi/2)}\biggl{(}\int_{0}^{t}\langle\phi,\mathcal{I}^{\alpha}\phi\rangle\,ds\biggr{)}^{1/2}\biggl{(}\int_{0}^{t}\langle\psi,\mathcal{I}^{\alpha}\psi\rangle\,ds\biggr{)}^{1/2},$$

$${\int }_0^t⟨\varphi ,{\mathcal{I}}^{\alpha }\varphi ⟩ds\ge \cos (\pi \alpha /2){\int }_0^t\Vert {\mathcal{I}}^{\alpha /2}\varphi {\Vert }^{2}ds,$$

$${\int }_0^t\Vert {\mathcal{I}}^{\alpha }\varphi {\Vert }^{2}ds\le C{t}^{\alpha }{\int }_0^t\Vert {\mathcal{I}}^{\alpha /2}\varphi {\Vert }^{2}ds,$$

$${\int }_0^t\Vert {\mathcal{I}}^{\alpha }\varphi {\Vert }^{2}ds\le C{t}^{2\alpha }{\int }_0^t\Vert \varphi {\Vert }^{2}ds,$$

$$\int_{0}^{t}\langle\phi,\mathcal{I}^{\alpha}\phi\rangle\,ds\leq\biggl{(}\int_{0}^{t}\|\phi\|^{2}\,ds\biggr{)}^{1/2}\biggl{(}\int_{0}^{t}\|\mathcal{I}^{\alpha}\phi\|^{2}\,ds\biggr{)}^{1/2}\leq Ct^{\alpha}\int_{0}^{t}\|\phi\|^{2}\,ds,$$

$${\int }_0^t\Vert {\mathcal{I}}^{\alpha }\varphi {\Vert }^2\le \frac{C{t}^{\alpha /2}}{1-\alpha }{\int }_0^t{\omega }_{\alpha /2}(t-s)y(s)ds\text{for}y(t)={\int }_0^t⟨\varphi ,{\mathcal{I}}^{\alpha }\varphi ⟩ds.$$

$${\int }_0^T\Vert {\mathcal{I}}^{\nu }\psi (t){\Vert }^{2}dt\le {\omega }_{\nu +1}(T){\int }_0^T{\omega }_{\nu }(T-t){\int }_0^t\Vert \psi (s){\Vert }^{2}dsdt,$$

$$\Vert \varphi (t)-\varphi (0){\Vert }^2\le \frac{t^{1-\alpha }}{(1-\alpha {)}^2}{\int }_0^t⟨{\varphi }^{\prime }(s),({\mathcal{I}}^{\alpha }{\varphi }^{\prime })(s)⟩ds.$$

$$\Vert \varphi (t)-\varphi (0){\Vert }^2\le \frac{t^{1-\alpha }}{1-\alpha }{\int }_0^t\Vert {\mathcal{I}}^{\alpha /2}{\varphi }^{\prime }(s){\Vert }^{2}ds,$$

$$\mathcal{M}\varphi (t)$$

$${\vec{B}}_2(\varphi )$$

$$\mathcal{M}{\mathcal{I}}^{\alpha }-{\mathcal{I}}^{\alpha }\mathcal{M}=\alpha {\mathcal{I}}^{\alpha +1}$$

$$({\partial }_t^{1-\alpha }\varphi )(t)=({\mathcal{I}}^{\alpha }\varphi {)}^{\prime }=\varphi (0){\omega }_{\alpha }(t)+({\mathcal{I}}^{\alpha }{\varphi }^{\prime })(t).$$

Full text

1 Introduction

We consider the spatial discretisation via Galerkin finite elements of a time-fractional Fokker–Planck equation AngstmannEtAl2015 ; HenryLanglandsStraka2010 ,

[TABLE]

with initial condition $u(\vec{x},0)=u_{0}(\vec{x}),$ where $\partial_{t}=\partial/\partial t$ and $\Omega$ is a polyhedral domain in $\mathbb{R}^{d}$ ($d\geq 1$). The fractional exponent is restricted to the range $0<\alpha<1$, $\kappa_{\alpha}>0$ is the diffusivity coefficient. In our analysis, we put $\kappa_{\alpha}=1$ for convenience, but it is straight forward to extend our methods to allow for a spatially-varying diffusivity. The fractional derivative is taken in the Riemann–Liouville sense, that is, $\partial_{t}^{1-\alpha}u=\partial_{t}\mathcal{I}^{\alpha}u$, where the fractional integration operator $\mathcal{I}^{\alpha}$ is defined by

[TABLE]

Though we impose a homogeneous Dirichlet boundary condition,

[TABLE]

the proposed stability and errors analysis remain valid for zero-flux boundary condition, see Remark 1.

The time-space dependent driving force $\vec{F}$ and it time partial derivative, $\partial_{t}\vec{F}$, are assumed to be in  L_{\infty}\bigl{(}\Omega\times(0,T),\mathbb{R}^{d}\bigr{)}. When $\vec{F}$ is independent of $t$, the model problem (1) can be rewritten in the form

[TABLE]

where the first term is just the Caputo fractional derivative of order $\alpha$. For a one- or two-dimensional spatial domain $\Omega$, numerical methods applicable to (3) have been widely studied CaoFuHuang2012 ; ChenLiuZhuangAnh2009 ; Cui2015 ; Deng2007 ; FairweatherZhangYangXu2015 ; GaoSun2015a ; GaoSun2015b ; GraciaORiordanStynes2015 ; Jiang2015 ; SaadatmandiDehghanAzizi2012 ; VongWang2015 ; WeiZhangHe2013 ; ZhengLiuTurnerAnh2015 . In all of these works, the solution $u$ was assumed to be sufficiently regular, including at $t=0$. Although (3) is in many respects more convenient for constructing and analyzing the accuracy of numerical schemes, only (1) is physically valid for a time-dependent forcing $\vec{F}$ HeinsaluPatriarcaGoychukHanggi2007 .

Our earlier paper KimMcLeanMustapha2016 presented an analysis of the semidiscrete finite element solution of (1) that is limited to cases in which

the solution $u$ is sufficiently regular, 2. 2.

the spatial domain $\Omega$ is an interval on the real line (that is, $d=1$), 3. 3.

the fractional exponent is in the range $1/2<\alpha<1$, 4. 4.

the boundary condition is of homogeneous Dirichlet type (2).

By employing a different approach that based on novel energy arguments, we are able to relax significantly the regularity requirements on $u$, in addition to permitting $d\geq 1$, $0<\alpha<1$, and zero-flux (10) as well as Dirichlet boundary conditions. This new approach leads to an error bound of optimal order in $L_{2}(\Omega)$ at each fixed $t>0$, even for non-smooth initial data $u_{0}$. We consider only continuous piecewise linear elements and (unlike our earlier paper KimMcLeanMustapha2016 ) do not analyse any time discretisation.

In Section 2, we define the semidiscrete finite element scheme and outline our main results in the context of our previous work KimMcLeanMustapha2016 . Section 3 gathers together some technical estimates involving fractional integrals. Section 4 presents the new stability result (Theorem 4.1) and Section 5 the new error bound (Theorem 5.1). Finally, in Section 6, we discuss two numerical examples. The first confirms both the convergence rate and the dependence on $t$ predicted by our theory. The second looks briefly at how the method behaves when $u_{0}$ is a point mass, and therefore does not even belong to $L_{2}(\Omega)$.

2 The finite element solution

The continuous solution $u:(0,T]\to H^{1}_{0}(\Omega)$ of problem (1) subject to the homogeneous Dirichlet boundary condition (2), satisfies the weak form,

[TABLE]

for all $v\in H^{1}_{0}(\Omega)$, where $\langle u,v\rangle=\int_{\Omega}uv$ and $\langle\vec{u},\vec{v}\rangle=\int_{\Omega}\vec{u}\cdot\vec{v}$. Let $h$ denote the maximum element diameter from a shape-regular triangulation of $\Omega$, and let $\mathbb{S}_{h}\subseteq H^{1}_{0}(\Omega)$ denote the usual space of continuous, piecewise-linear functions that vanish on $\partial\Omega$. The semidiscrete finite element solution $u_{h}:[0,T]\to\mathbb{S}_{h}$ is then defined by

[TABLE]

together with the initial condition $u_{h}(0)=u_{0h}$, where $u_{0h}\in\mathbb{S}_{h}$ is a suitable approximation to $u_{0}$.

Previously, for $0\leq t\leq T$, we showed (KimMcLeanMustapha2016, , Theorems 3.3 and 3.4) that, $\|u_{h}(t)\|\leq C\|u_{0h}\|_{1}$ and, provided $u_{0h}$ is chosen to be the Ritz projection of $u_{0}$ onto $\mathbb{S}_{h}$,

[TABLE]

Here, $\|v\|=\sqrt{\langle v,v\rangle}$ denotes the norm in $L_{2}(\Omega)$, $u^{\prime}(t)=\partial_{t}u$,

[TABLE]

and $\varphi_{1}$, $\varphi_{2}$, $\varphi_{3}$, …is a complete orthonormal system in $L_{2}(\Omega)$ consisting of Dirichlet eigenfunctions of the Laplacian: $\langle\varphi_{m},\varphi_{k}\rangle=\delta_{mk}$ and

[TABLE]

The associated function space $\dot{H}^{r}(\Omega)=\{\,v\in L_{2}(\Omega):\|v\|_{r}<\infty\,\}$ is a subspace of the usual Sobolev space $H^{r}(\Omega)$ for $0\leq r\leq 1$; in particular, $\dot{H}^{0}(\Omega)=L_{2}(\Omega)$ and $\dot{H}^{1}(\Omega)=H^{1}_{0}(\Omega)$. Also, $\dot{H}^{2}(\Omega)=H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$ provided $\Omega$ is convex (so the Poisson problem is $H^{2}$-regular).

We prove in Theorem 4.1 a stronger stability estimate,

[TABLE]

Also, whereas the previous error bound (6) is meaningful only if $u_{0}\in\dot{H}^{2}(\Omega)$ and u^{\prime}\in L_{2}\bigl{(}(0,T),\dot{H}^{2}(\Omega)\bigr{)}, our new error analysis makes a much weaker regularity assumption: for some $r$ in the range $0\leq r\leq 2$ there is a constant $K_{r}$ such that

[TABLE]

When $\vec{F}\equiv\vec{0}$ and the domain $\Omega$ is convex, it is known (McLean2010, , Theorem 4.4) that such an estimate holds with $K_{r}=C\|u_{0}\|_{r}$ in the case of Dirichlet boundary conditions (2). Since the term of (1) involving $\vec{F}$ is of lower order in the spatial variables, we conjecture that the same is true for a nonzero (but sufficiently regular) forcing $\vec{F}$. In Theorem 5.1, we show that if $u_{0h}$ is chosen to be the $L_{2}$-projection of $u_{0}$ onto $\mathbb{S}_{h}$, then

[TABLE]

For instance, in the worst case when $r=0$, the error is $O(t^{-\alpha}h^{2})$.

Remark 1

If we impose a zero-flux boundary condition,

[TABLE]

where $\vec{n}$ denotes the outward unit normal to $\Omega$, then $u:(0,T]\to H^{1}(\Omega)$ satisfies (4) for all $v\in H^{1}(\Omega)$. Likewise, $u_{h}$ is defined as in (5) but the finite element space $\mathbb{S}_{h}\subseteq H^{1}(\Omega)$ now consists of all continuous piecewise-linear functions (that is, the elements of $\mathbb{S}_{h}$ need not vanish on $\partial\Omega$). The stability estimate (7) remains valid, and the error bound (9) holds assuming $u$ satisfies (8), where $\|\cdot\|_{2}$ is now the norm in $H^{2}(\Omega)$ rather than $\dot{H}^{2}(\Omega)$. Note that for either choice of boundary condition, the variational equation (5) is equivalent to a system of Volterra integral equations (KimMcLeanMustapha2016, , Theorem 3.1) that admits a unique continuous solution $u_{h}:[0,T]\to\mathbb{S}_{h}$. Moreover, the methods of Miller and Feldstein (MillerFeldstein1971, , Theorem 1) show that $u_{h}$ is continuously differentiable on $(0,T]$. Finally, notice that in the case of the zero-flux boundary condition (10), the total mass $\int_{\Omega}u(\cdot,t)$ within $\Omega$ is conserved.

3 Fractional integrals

In this section only, $C$ is an absolute constant. Our analysis of the semidiscrete finite element solution $u_{h}$ will rely on the following technical lemmas, in which $\phi$ and $\psi$ are suitably regular functions of $t>0$ taking values in a Hilbert space.

Lemma 1

If  $0\leq\mu\leq\nu\leq 1$, then

[TABLE]

Proof

If $\mu=\nu$ then there is nothing to prove, so assume $\mu<\nu$. In a previous paper (KimMcLeanMustapha2016, , Lemma 2.3), we showed that for $0<\alpha\leq 1$,

[TABLE]

and the right-hand side is bounded by $\omega_{\alpha+1}(T)^{2}\int_{0}^{T}\|\psi(s)\|^{2}\,ds$. Putting $\psi=\mathcal{I}^{\mu}\phi$ and $\alpha=\nu-\mu$, it follows that $\mathcal{I}^{\alpha}\psi=\mathcal{I}^{\nu}\phi$ and $\omega_{\alpha+1}(T)\leq CT^{\alpha}=CT^{\nu-\mu}$. ∎

Lemma 2

If  $0<\alpha<1$ and $\epsilon>0$, then

[TABLE]

[TABLE]

[TABLE]

Proof

From a result of Mustapha and Schötzau (MustaphaSchoetzau2014, , Lemma 3.1(iii)),

[TABLE]

so (11) follows because $\cos(\alpha\pi/2)\geq 1-\alpha$. The same paper (MustaphaSchoetzau2014, , Lemma 3.1(ii)) showed that

[TABLE]

and by choosing $\nu=\alpha$ and $\mu=\alpha/2$ in Lemma 1 have

[TABLE]

proving (12). Instead choosing $\nu=\alpha$ and $\mu=0$ in Lemma 1 gives

[TABLE]

so

[TABLE]

proving (13).∎

Lemma 3

If  $0<\alpha<1$, then

[TABLE]

Proof

From our earlier paper (KimMcLeanMustapha2016, , Lemma 2.3),

[TABLE]

so the result follows by letting $\nu=\alpha/2$ and $\psi=\mathcal{I}^{\alpha/2}\phi$, and then using (14). ∎

Lemma 4

If  $0<\alpha<1$, then

[TABLE]

Proof

We showed previously (KimMcLeanMustapha2016, , Lemma 2.1) that

[TABLE]

so the desired estimate follows from (14) and the inequality $\cos(\alpha\pi/2)\geq 1-\alpha$.∎

4 Stability

We seek to estimate the finite element solution $u_{h}(t)$ in terms of the initial data $u_{0h}$. Throughout, the generic constant $C$ may depend on $\alpha$, $T$ and the vector norms of $\vec{F}$ and $\vec{F}^{\prime}=\partial_{t}\vec{F}$ in $L_{\infty}(\Omega\times(0,T))$.

It will be convenient to define

[TABLE]

and we will use the elementary identities

[TABLE]

and

[TABLE]

Lemma 5

For $0\leq t\leq T$,

[TABLE]

Proof

Integration by parts (in time) shows that

[TABLE]

and our assumptions on $\vec{F}$ imply

[TABLE]

so the first estimate follows at once. The second estimate follows immediately from the first one and the inequality

[TABLE]

With the help of the identities (18) and (16), we find that

[TABLE]

so

[TABLE]

Thus,

[TABLE]

which implies the third estimate.∎

In the next two lemmas, we prove preliminary stability estimates for $u_{h}$ and $\mathcal{M}u_{h}$.

Lemma 6

The finite element solution satisfies, for $0\leq t\leq T$,

[TABLE]

and

[TABLE]

Proof

We integrate (5) in time to obtain

[TABLE]

and then choose $\chi=\mathcal{I}^{\alpha}u_{h}(t)$ so that

[TABLE]

Therefore, after cancelling the term $\tfrac{1}{2}\|\mathcal{I}^{\alpha}\nabla u_{h}\|^{2}$, integrating in time and applying Lemma 5, we deduce that

[TABLE]

From (11) with $\phi=u_{0h}$ and $\psi=u_{h}$,

[TABLE]

so if we define

[TABLE]

then

[TABLE]

Noting that $(\mathcal{I}^{\alpha}u_{0h})(t)=u_{0h}\omega_{\alpha+1}(t)$, and applying Lemma 3 with $\phi=u_{h}$, it follows that

[TABLE]

where

[TABLE]

Let $E_{\beta}(z)=\sum_{n=0}^{\infty}z^{n}/\Gamma(1+n\beta)$ denote the Mittag–Leffler function. A generalised Gronwall inequality of Dixon and McKee (DixonMcKee1986, , Theorem 3.1) (also stated in our earlier paper (KimMcLeanMustapha2016, , Lemma 2.6)) then yields

[TABLE]

The first estimate of the lemma follows at once, and the second is then a consequence of (12).∎

Lemma 7

For $0\leq t\leq T$,

[TABLE]

and

[TABLE]

Proof

We multiply both sides of (20) by $t$, and then use (16), to obtain

[TABLE]

By integrating (20) in time, we find that

[TABLE]

and so, noting that $\mathcal{I}^{1}u_{0h}=\mathcal{M}u_{0h}$,

[TABLE]

Now choose $\chi=\mathcal{I}^{\alpha}\mathcal{M}u_{h}$, cancel the term $\tfrac{1}{2}\|\nabla\chi\|^{2}$ and integrate in time to arrive at the estimate

[TABLE]

Using (11) twice, with $\epsilon=1/4$, we see that the second term on the right-hand side is bounded by

[TABLE]

so

[TABLE]

Since $\mathcal{I}^{\alpha}\mathcal{M}u_{0h}=u_{0h}\mathcal{I}^{\alpha}\omega_{2}=u_{0h}\omega_{\alpha+2}$, we have

[TABLE]

and, using (13) followed by Lemma 1 with $\nu=1$ and $\mu=\alpha$,

[TABLE]

Thus, by Lemma 5,

[TABLE]

which, when combined with the second estimate from Lemma 6, proves the first claim. The second follows at once thanks to (12).∎

Next, we show that $u_{h}$ may be replaced with $(\mathcal{M}u_{h})^{\prime}$ in the first estimate of Lemma 6.

Lemma 8

For $0\leq t\leq T$,

[TABLE]

Proof

Differentiate (24) to obtain

[TABLE]

and note that

[TABLE]

We choose $\chi=\partial_{t}^{1-\alpha}\mathcal{M}u_{h}=(\mathcal{I}^{\alpha}\mathcal{M}u_{h})^{\prime}$, and observe that $(\mathcal{M}u_{h})(0)=0$ so (17) implies $\chi=\mathcal{I}^{\alpha}(\mathcal{M}u_{h})^{\prime}$. Thus,

[TABLE]

By (11),

[TABLE]

so by Lemma 5,

[TABLE]

The first integral on the right-hand side is bounded by $Ct^{1+\alpha}\|u_{0h}\|^{2}$, and so is the second via Lemmas 6 and 7. It follows using Lemma 3 that $y(t)$ satisfies an inequality of the form (22) with $a(t)=Ct^{1+\alpha}\|u_{0h}\|^{2}$ and $b(t)=Ct^{\alpha/2}$, so (23) holds, proving the result. ∎

The stability of $u_{h}(t)$ in $L_{2}(\Omega)$ now follows.

Theorem 4.1

There is a constant $C$, depending on $\alpha$, $T$ and $\vec{F}$, such that

[TABLE]

Proof

Using Lemma 4 with $\phi=\mathcal{M}u_{h}$, followed by Lemma 8, we obtain

[TABLE]

∎

Because some of the estimates of Section 3 break down as $\alpha\to 1$, the same is true of the stability result above. That is, the proof of Theorem 4.1 yields a constant $C$ that tends to infinity as $\alpha\to 1$. However, we can easily prove stability in the limiting case when $\alpha=1$, that is, when (1) reduces to the classical Fokker–Planck equation,

[TABLE]

and the finite element equation (5) to

[TABLE]

5 Error estimate

We now seek to estimate the accuracy of the semidiscrete finite element solution $u_{h}$. Recall that the Ritz projection $R_{h}v\in\mathbb{S}_{h}$ of a function $v\in H^{1}(\Omega)$ is defined by

[TABLE]

here, the lower-order terms are included to allow for a zero-flux boundary condition (10), in which case the functions in $\mathbb{S}_{h}$ do not have to vanish on $\partial\Omega$ and so the Poincaré inequality is not applicable. Since the Galerkin finite element method is quasi-optimal in $H^{1}(\Omega)$, we know that $\|v-R_{h}v\|_{1}\leq Ch\|v\|_{2}$ for $v\in H^{2}(\Omega)$. Assuming that $\Omega$ is convex, so that the Poisson problem is $H^{2}$-regular, the usual duality argument implies that

[TABLE]

We now decompose the error into

[TABLE]

and deduce from (4) and (5) that

[TABLE]

With this equation, we can use the techniques of Section 4 to estimate $\theta_{h}$ in terms of $\rho_{h}$. The next lemma provides our basic estimate for the latter.

Lemma 9

Let $\beta\geq 0$ and $0\leq r\leq 2$. If $u$ has the regularity property (8), then

[TABLE]

Proof

For the case $\beta=0$, we see from (25) that

[TABLE]

whereas for $\beta>0$,

[TABLE]

and the result follows after making the substitution $s=ty$ for $0\leq y\leq 1$. ∎

The proofs of Lemmas 10 and 11 below parallel those of Lemmas 6 and 7 from Section 4. We let $P_{h}$ denote $L_{2}$-projector onto the finite element subspace $\mathbb{S}_{h}$, that is, for any $v\in L_{2}(\Omega)$ we define $P_{h}v\in\mathbb{S}_{h}$ by $\langle P_{h}v,\chi\rangle=\langle v,\chi\rangle$ for all $\chi\in\mathbb{S}_{h}$.

Lemma 10

If $u_{0h}=P_{h}u_{0}$ then, for $0\leq t\leq T$ and $0\leq r\leq 2$,

[TABLE]

and

[TABLE]

Proof

We integrate (27) in time to obtain

[TABLE]

where $\tilde{\rho}_{h}=\rho_{h}-\mathcal{I}^{\alpha}\rho_{h}$. Our choice of $u_{0h}$ means that $\langle e_{h}(0),\chi\rangle=0$, so by letting $\chi=\mathcal{I}^{\alpha}\theta_{h}$ and recalling the definitions (15), we see that

[TABLE]

Thus, by Lemma 5,

[TABLE]

After applying (11) with $\phi=\tilde{\rho}_{h}$ and $\psi=\theta_{h}$, followed by Lemma 3 with $\phi=\theta_{h}$, we see that the function

[TABLE]

satisfies an inequality of the form (22) with

[TABLE]

For brevity, put $\eta=h^{2}K_{r}$. By Lemma 9,

[TABLE]

and \|\mathcal{I}^{\alpha}\rho_{h}\|^{2}\leq C\bigl{(}\eta t^{\alpha+\alpha(r-2)/2}\bigr{)}^{2}=C\eta^{2}t^{\alpha r}, so $a(t)\leq C\eta^{2}t^{\alpha(r-1)+1}$. Thus, the two estimates follow from (23) followed by (12).∎

Lemma 11

If $u_{0h}=P_{h}u_{0}$ then, for $0\leq t\leq T$ and $0\leq r\leq 2$,

[TABLE]

and

[TABLE]

Proof

We multiply both sides of (28) by $t$, remembering that $\langle e_{h}(0),\chi\rangle=0$, and then use (16) to obtain

[TABLE]

By integrating (28), we find that

[TABLE]

and hence, with $\vec{B}_{2}(\phi)$ defined as before in (15),

[TABLE]

Now choose $\chi=\mathcal{I}^{\alpha}\mathcal{M}\theta_{h}$ so that, after cancelling a term $\tfrac{1}{2}\|\nabla\chi\|^{2}$ and integrating,

[TABLE]

Using (11) with $\epsilon=1/4$, $\phi=(\mathcal{M}-\alpha\mathcal{I}^{1})\tilde{\rho}_{h}$ and $\psi=\mathcal{M}\theta_{h}$, and a second time with $\phi=\alpha\mathcal{I}^{1}\theta_{h}$, we see that

[TABLE]

Lemma 5 implies that

[TABLE]

and, putting $\eta=h^{2}K_{r}$ as before, we find with the help of Lemma 9 that

[TABLE]

Using (13), followed by Lemma 1 with $\nu=1$ and $\mu=\alpha$,

[TABLE]

so, recalling that $\|\mathcal{I}^{\alpha}\rho_{h}\|^{2}\leq C\eta^{2}t^{\alpha r}$, the first estimate follows by Lemma 10. The second is then an immediate consequence of (12).

Techniques like those of Lemma 8 and Theorem 4.1 now yield our error bound.

Theorem 5.1

If $\Omega$ is convex and the solution of the fractional Fokker–Planck equation (1) has the regularity property (8), then the finite element solution, given by (4), satisfies

[TABLE]

for $0<t\leq T$ and $0\leq r\leq 2$. The constant $C$ may depend on $\alpha$, $T$ and $\vec{F}$.

Proof

Suppose in the first instance that $u_{0h}=P_{h}u_{0}$, as required for Lemmas 10 and 11. Differentiate (29) to obtain

[TABLE]

where $\vec{B}_{3}(\phi)$ is again defined as in (15). Noting that

[TABLE]

we choose $\chi=\partial_{t}^{1-\alpha}\mathcal{M}\theta_{h}=(\mathcal{I}^{\alpha}\mathcal{M}\theta_{h})^{\prime}$, and observe that $(\mathcal{M}\theta_{h})(0)=0$ so (17) implies $\chi=\mathcal{I}^{\alpha}(\mathcal{M}\theta_{h})^{\prime}$. Thus, after cancelling $\|\nabla\chi\|^{2}$,

[TABLE]

Integrating in time, and then applying (11) to the first term on the right hand side, with $\epsilon=1/2$, $\phi=(\mathcal{M}\tilde{\rho}_{h})^{\prime}$ and $\psi=(\mathcal{M}\theta_{h})^{\prime}$, it follows that

[TABLE]

Since, using (16),

[TABLE]

we see from (25), (8) and Lemma 9 that $\|(\mathcal{M}\tilde{\rho}_{h})^{\prime}\|\leq C\eta t^{\alpha(r-2)/2}(1+t^{\alpha})$ where, as before, $\eta=h^{2}K_{r}$. Consequently,

[TABLE]

and by Lemma 5,

[TABLE]

showing that

[TABLE]

Using Lemmas 10 and 11, we find that the second term on the right is bounded by $Ct^{1+\alpha(r-1)}\eta^{2}$. It follows using Lemma 3 that the function

[TABLE]

satisfies an inequality of the form (22) with $a(t)=Ct^{1+\alpha(r-1)}\eta^{2}$ and $b(t)=Ct^{\alpha/2}$. Therefore, using Lemma 4 with $\phi=\mathcal{M}\theta_{h}$, followed by (23), we have

[TABLE]

which is equivalent to the estimate $\|\theta_{h}\|\leq Ct^{\alpha(r-2)/2}h^{2}K_{r}$. Recalling (26), the desired error bound in the case $u_{0h}=P_{h}u_{0}$ follows by the triangle inequality and the case $\beta=0$ of Lemma 9.

The error bound for general $u_{0h}$ now follows from the stability result of Theorem 4.1. In fact, if $u_{h}^{*}$ and $u_{h}$ denote the finite element solutions satisfying $u_{h}^{*}(0)=P_{h}u_{0}$ and $u_{h}(0)=u_{0h}$, then the difference $u_{h}-u^{*}_{h}$ is the finite element solution with initial value $u_{0h}-P_{h}u_{0}$ so

[TABLE]

We obtain the desired estimate for $\|u_{h}(t)-u(t)\|$ after applying the triangle inequality, noting that $\|u_{h}^{*}(t)-u(t)\|\leq Ct^{\alpha(r-2)/2}h^{2}K_{r}$.∎

If $r<2$, then the error estimate in the theorem becomes unbounded as $t\to 0$, but the stability result of Theorem 4.1 shows that the error must in fact remain bounded.

6 Numerical examples

We discuss experiments with two problems, using a fully-discrete scheme of implicit Euler type. For time levels $0=t_{0}<t_{1}<t_{2}<\cdots<t_{N}=T$, we denote the $n$th step size by $k_{n}=t_{n}-t_{n-1}$ and the associated subinterval by $I_{n}=(t_{n-1},t_{n})$, for $1\leq n\leq N$. The maximum step size $k=\max_{1\leq n\leq N}k_{n}$ is sometimes used to label quantities that depend on the mesh. With any sequence of values $V^{1}$, $V^{2}$, …, $V^{N}$ we associate the piecewise-constant function $\check{V}$ defined by

[TABLE]

Integrating the finite element equation (5) over the $n$th time interval $I_{n}$ gives

[TABLE]

for all $\chi\in\mathbb{S}_{h}$, and we approximate $u_{h}(t_{n})$ by $U^{n}_{h}\in\mathbb{S}_{h}$ satisfying

[TABLE]

for all $\chi\in\mathbb{S}_{h}$ and for $1\leq n\leq N$, with $U^{0}_{h}=u_{0h}$. For $1\leq p\leq Q_{h}:=\dim\mathbb{S}_{h}$, let $\vec{x}_{p}$ denote the $p$th free node of the spatial mesh, and let $\phi_{p}\in\mathbb{S}_{h}$ denote the $p$th nodal basis function, so that $\phi_{p}(\vec{x}_{q})=\delta_{pq}$ and

[TABLE]

We define $Q_{h}\times Q_{h}$ matrices $\vec{M}$ and $\vec{G}^{n}$ with entries

[TABLE]

where $\vec{F}^{n}(\vec{x})=\vec{F}(\vec{x},t_{n})$, and the $Q_{h}$-dimensional column vector $\vec{U}^{n}$ with components $U^{n}_{p}$. It follows from (30) that

[TABLE]

with weights $\omega_{nj}=\int_{I_{j}}\omega_{\alpha}(t_{n}-s)\,ds$ for $1\leq j\leq n\leq N$. Thus, at the $n$th time step we must solve the linear system

[TABLE]

Although this fully-discrete scheme lacks a theoretical error analysis, we observed numerically that first-order accuracy in time is achieved, for $t$ bounded away from zero, if we use a graded mesh of the form

[TABLE]

Our earlier paper (KimMcLeanMustapha2016, , Table 5.3) includes computations with smooth initial data, in which we observed that the $L_{2}$ error is $O(h^{2})$ uniformly for $0\leq t\leq T$, consistent with Theorem 5.1 when $r=2$. Here, we instead focus on the case of non-smooth initial data.

6.1 Dirichlet boundary condition

In our first example, $F(x,t)=-x+\sin t$, $T=1$ and $\Omega=(0,\pi)$, with homogeneous Dirichlet boundary conditions $u(0,t)=0=u(\pi,t)$ and discontinuous initial data given by

[TABLE]

Figure 1 shows $u_{0}$ and its $L_{2}$-projection $P_{h}u_{0}$, as well as the nodal interpolant $I_{h}u_{0}\in\mathbb{S}_{h}$ defined by

[TABLE]

The Dirichlet eigenvalues and orthonormal eigenfunctions of $-\nabla^{2}=-\partial_{x}^{2}$ are

[TABLE]

so for $0\leq r<1/2$ we have

[TABLE]

If our conjecture that $K_{r}=C\|u_{0}\|_{r}$ in (8) is valid, then applying Theorem 5.1 with $r=\tfrac{1}{2}-\epsilon$ and $\epsilon^{-1}=\log(e^{2}+t^{-1})$, so that $t^{-\epsilon}\leq e$ and $0<\epsilon<1/2$, gives

[TABLE]

In our computations, we employed nonuniform time levels given by (31), but a uniform spatial mesh with $h=1/(Q_{h}+1)$. In all cases, $Q_{h}+1$ was divisible by $4$ so that the points $\pi/4$ and $3\pi/4$ (where $u_{0}$ is discontinuous) coincided with two of the nodes. We first computed a reference solution $U_{\text{ref}}^{n}=U^{n}_{h}$ using a fine mesh with $N=10,000$ and $Q_{h}=511$. We then computed $U^{n}_{h}$ for $Q_{h}\in\{7,15,31,63\}$, again with $N=10,000$. The initial data was chosen as $u_{0h}=P_{h}u_{0}$ in each case. With such a small $k$, the error,

[TABLE]

was dominated by the influence of the spatial discretisation, and we sought to estimate the convergence rates $\sigma_{h,k}$ such that

[TABLE]

from the relation

[TABLE]

Table 1 shows the values of $E^{*}_{h,k}$ and $\sigma_{h,k}$ for three different values of $\alpha$. The computed values of $\sigma_{h,k}$ are close to $2$, as expected from Theorem 5.1. Figure 2 shows how the $L_{2}$-error $E_{h,k}^{n}$ varies with $t_{n}$ for different $h$ when $\alpha=0.75$, again keeping $N=10,000$. Due to the log-log scale, the graph of a function proportional to $t^{-3\alpha/4}$ appears as a straight line with gradient $-3\alpha/4$, indicated by the small triangle, and we observe exactly this behaviour of the error for $t$ close—but not too close—to zero.

Physically, the solution $u$ must be non-negative, but the oscillations in the discrete initial data $P_{h}u_{0}$ mean that $U^{n}_{h}(x)$ was negative for some values of $(x,t_{n})$ near the points of discontinuity $(\pi/4,0)$ and $(3\pi/4,0)$. It is tempting to choose as the discrete initial data $u_{0h}=I_{h}u_{0}$, the nodal interpolant (33). In this way, $U^{0}_{h}=u_{0h}(x)\geq 0$ for all $x$. However, since

[TABLE]

by choosing $\chi=u_{0h}-P_{h}u_{0}$ we see that

[TABLE]

Thus, Theorem 5.1 now yields an error bound of order $h+t^{-3\alpha/4}h^{2}$ (ignoring the log factor), and Table 2 indeed shows only first-order convergence for this choice of initial data.

At the end of Section 4, we remarked that in our stability estimate the constant tends to infinity as $\alpha$ approaches $1$. Since the finite element method is stable in the classical case $\alpha=1$, we suspect that the dependence of the stability constant on $\alpha<1$ is an artefact of the method of proof. To investigate this question numerically, we computed $\|u_{h}(t)\|$ for random initial data, that is, when the value of $u_{0h}$ at each node was a random number from a uniform distribution in $[0,1]$. In practice, we did not observe any deterioration in the stability of the method for $\alpha$ close to $1$.

6.2 Zero-flux boundary condition

In our second example,

[TABLE]

where $V$ is a double-well potential perturbed by an oscillation in time,

[TABLE]

Gammaitoni et al. GammaitoniEtAl1998 used this potential for the classical Fokker–Planck equation ($\alpha=1$) in their study of stochastic resonance. We imposed the zero-flux boundary condition (10) and chose as the initial data $u_{0}(x)=\delta(x)$. The solution $u$ then gives the probability distribution for a single diffusing particle initially located at $x=0$. Since the Dirac delta functional does not belong to $L_{2}(\Omega)$, our stability result (Theorem 4.1) does not apply, and $P_{h}u_{0}$ is not defined. Nevertheless, the functions in $\mathbb{S}_{h}$ are continuous, so by extending the $L_{2}$ inner product to a dual pairing we can define the discrete initial data $u_{0h}\in\mathbb{S}_{h}$ by

[TABLE]

Figure 3 shows a surface plot of the numerical solution using $N=\text{4,096}$ time steps, now with a stronger mesh grading $\gamma=2$ in (31), and $Q_{h}=65$ spatial degrees of freedom. (Thus the delta function is centred on the node $\vec{x}_{33}=0$). We cut off the initial part of the plot where $t<0.005$ to avoid the oscillations, shown separately in Figure 4, which are much larger than was the case for our first example. The total mass should be constant and we observed in practice that $\int_{\Omega}U^{n}_{h}=1$ to ten significant figures, for $0\leq n\leq N$.