On an explicit finite difference method for fractional diffusion equations

TL;DR

This paper develops an explicit FTCS-type scheme for the fractional diffusion equation $\frac{\partial}{\partial t}u(x,t)=K_\gamma\,{}_0D_t^{1-\gamma}\frac{\partial^2}{\partial x^2}u(x,t)$ by discretizing the Grünwald-Letnikov fractional derivative. The resulting update $U_j^{(m+1)}=U_j^{(m)}+S_\gamma\sum_{k=0}^{m}\omega_k^{(1-\gamma)}[U_{j-1}^{(m-k)}-2U_j^{(m-k)}+U_{j+1}^{(m-k)}]$ with $S_\gamma=K_\gamma(\Delta t)^\gamma/(\Delta x)^2$ is non-Markovian and memory-intensive but highly explicit. A von Neumann stability analysis yields a bound $S_\gamma \le S_\gamma^{\times}$, with $S_\gamma^{\times}=1/2^{2-\gamma}$ for first-order and $S_\gamma^{\times}=1/4^{3/2-\gamma}$ for second-order approximations, consistent with numerical tests; truncation error shows the scheme is consistent with $T(x,t)=O(h^p)+O(\Delta t)+O((\Delta x)^2)$. Numerical experiments against exact subdiffusive solutions and absorbing-boundary problems demonstrate high accuracy near the stability bound, and the study discusses memory considerations and the potential for short-memory approaches or implicit schemes for long-time computations. The method generalizes to higher dimensions and other fractional dynamics problems, offering a transparent explicit alternative to implicit schemes in many scenarios.

Abstract

A numerical method to solve the fractional diffusion equation, which could also be easily extended to many other fractional dynamics equations, is considered. These fractional equations have been proposed in order to describe anomalous transport characterized by non-Markovian kinetics and the breakdown of Fick's law. In this paper we combine the forward time centered space (FTCS) method, well known for the numerical integration of ordinary diffusion equations, with the Grunwald-Letnikov definition of the fractional derivative operator to obtain an explicit fractional FTCS scheme for solving the fractional diffusion equation. The resulting method is amenable to a stability analysis a la von Neumann. We show that the analytical stability bounds are in excellent agreement with numerical tests. Comparison between exact analytical solutions and numerical predictions are made.

Figures (6)

• Figure 1: First values of $S_{\gamma,m}$ versus $m$ for $\gamma=1/2$ when the first-order coefficients (circles) and second-order coefficients (squares) are used. The lines mark the corresponding limit values $S_{\gamma}^\times$ given by Eqs. \ref{['Sfirst']} and \ref{['Ssecond']}
• Figure 2: The difference $\Delta S_{\gamma}=S_{\gamma,2}^{\times}-S_{\gamma,1}^{\times}$ versus $\gamma$ when the first-order coefficients for $\omega_k^{(1-\gamma)}$ [c.f. Eq. \ref{['coefO1']}] (solid line) and second-order coefficients [c.f. Eq. \ref{['G2']}] (dotted line) are used.
• Figure 3: Comparison between the exact subdiffusion propagator (lines) and the numerical integration results for $\gamma=1/4$ (squares), $\gamma=1/2$ (circles), $\gamma=3/4$ (triangles) and $\gamma=1$ (crosses) and $t=10$.
• Figure 4: Numerical solution of the subdiffusion equation for the problem with absorbing boundary conditions, $u(0,t)=u(1,t)=0$, and initial condition $u(x,0)=x(1-x)$ versus the exact analytical result (lines) for $t=0.5$. The solution $u(x,t)$ is shown for $\gamma=0.5$ (triangles), $\gamma=0.75$ (squares) and $\gamma=1$ (circles).
• Figure 5: Values of $S_\gamma^\times$ corresponding to the onset of instability versus the subdiffusion exponent $\gamma$. The solid line is the prediction of the Fourier--von Neumann analysis and the symbols denote the results of the numerical tests with the criterion in Eq. (\ref{['scriter']}): stars, triangles and squares for the absorbing boundary problem with $u(x,0)=x(1-x)$ with $M=50$, $100$ and $1000$, respectively, and circles for the propagator with $M=1000$.