Fast evaluation of the Caputo fractional derivative and its applications to fractional diffusion equations

Abstract

We present an efficient algorithm for the evaluation of the Caputo fractional derivative $_0^C\!D_t^αf(t)$ of order $α\in (0,1)$, which can be expressed as a convolution of $f'(t)$ with the kernel $t^{-α}$. The algorithm is based on an efficient sum-of-exponentials approximation for the kernel $t^{-1-α}$ on the interval $[Δt, T]$ with a uniform absolute error $\varepsilon$, where the number of exponentials $N_{\text{exp}}$ needed is of the order $O\left(\log\frac{1}{\varepsilon}\left( \log\log\frac{1}{\varepsilon}+\log\frac{T}{Δt}\right) +\log\frac{1}{Δt}\left( \log\log\frac{1}{\varepsilon}+\log\frac{1}{Δt}\right) \right)$. As compared with the direct method, the resulting algorithm reduces the storage requirement from $O(N_T)$ to $O(N_{\text{exp}})$ and the overall computational cost from $O(N_T^2)$ to $O(N_TN_{\text{exp}})$ with $N_T$ the total number of time steps. Furthermore, when the fast evaluation scheme of the Caputo derivative is applied to solve the fractional diffusion equations, the resulting algorithm requires only $O(N_SN_{\text{exp}})$ storage and $O(N_SN_TN_{\text{exp}})$ work with $N_S$ the total number of points in space; whereas the direct methods require $O(N_SN_T$) storage and $O(N_SN_T^2)$ work. The complexity of both algorithms is nearly optimal since $N_{\text{exp}}$ is of the order $O(\log N_T)$ for $T\gg 1$ or $O(\log^2N_T)$ for $T\approx 1$ for fixed accuracy $\varepsilon$. We also present a detailed stability and error analysis of the new scheme for solving linear fractional diffusion equations. The performance of the new algorithm is illustrated via several numerical examples. Finally, the algorithm can be parallelized in a straightforward manner.