A Fast Time Stepping Method for Evaluating Fractional Integrals

We evaluate the fractional integral $I^\alpha[f](t)=\frac{1}{\Gamma(\alpha)}\int_0^t(t-\tau)^{\alpha-1}\,f(\tau)\,d\tau$, $0<\alpha<1$, at time steps $t=\Delta t,2\Delta t,\dots,N\Delta t$ by making use of the integral representation of the convolution kernel $t^{\alpha-1}=\frac{1}{\Gamma(1-\alpha)}\int_0^{\infty}e^{-\xi\,t}\,\xi^{-\alpha}\,d\xi$. We construct an efficient $Q$-point quadrature of this integral representation and use it as a part of a fast time stepping method. The new method has algorithmic complexity $O(NQ)$ and storage requirement $O(Q)$. The number of quadrature nodes $Q$ is independent of $N$ and grows like $O\bigl(\bigl(-\log\epsilon-\log\Delta t\bigr)^2\bigr)$, where $\epsilon$ is the quadrature error tolerance and $\Delta t$ is the size of the time step. The (possible) singularity of $f$ near $\tau=0$ is taken into account. This new method is particularly well-suited for long time simulations.