A second-order accurate numerical scheme for a time-fractional Fokker-Planck equation

A time-stepping $L1$ scheme for solving a time fractional Fokker-Planck equation of order $\alpha \in (0, 1)$, with a general driving force, is investigated. A stability bound for the semi-discrete solution is obtained for $\alpha\in(1/2,1)$ {via a novel and concise approach.} Our stability estimate is $\alpha$-robust in the sense that it remains valid in the limiting case where $\alpha$ approaches $1$ (when the model reduces to the classical Fokker-Planck equation), a limit that presents practical importance. Concerning the error analysis, we obtain an optimal second-order accurate estimate for $\alpha\in(1/2,1)$. A time-graded mesh is used to compensate for the singular behavior of the continuous solution near the origin. The $L1$ scheme is associated with a standard spatial Galerkin finite element discretization to numerically support our theoretical contributions. We employ the resulting fully-discrete computable numerical scheme to perform some numerical tests. These tests suggest that the imposed time-graded meshes assumption could be further relaxed, and we observe second-order accuracy even for the case $\alpha\in(0,1/2]$, that is, outside the range covered by the theory.