Weighted average finite difference methods for fractional diffusion equations

TL;DR

This work develops weighted-average finite-difference schemes for fractional diffusion equations with a Riemann–Liouville derivative, combining standard time/space discretization with RL discretizations via weight sequences whose generating functions recover BDF1/2/3 and NG$p$ forms. A simple von Neumann–type stability analysis yields the bound $\dfrac{1}{S} \ge \dfrac{1}{S_\times}$ with $\dfrac{1}{S_\times}=2(2\lambda-1)\omega(-1,1-\gamma)$, implying unconditional stability for $\lambda\le 1/2$ and conditional stability otherwise, with explicit and Crank–Nicolson schemes examined. The paper validates the stability bound through numerical experiments against exact solutions expressed with Mittag-Leffler functions and discusses truncation errors, showing Crank–Nicolson as particularly robust (second-order) and suggesting a practical explicit-then-implicit integration strategy. The results offer a straightforward, stable framework for simulating fractional diffusion processes and can extend to higher dimensions and fractional wave equations.

Abstract

Weighted averaged finite difference methods for solving fractional diffusion equations are discussed and different formulae of the discretization of the Riemann-Liouville derivative are considered. The stability analysis of the different numerical schemes is carried out by means of a procedure close to the well-known von Neumann method of ordinary diffusion equations. The stability bounds are easily found and checked in some representative examples.

Figures (7)

• Figure 1: Stability phase diagram for the weighted average methods versus the weight factor $\lambda$. The line corresponds to Eq. \ref{['1StimesMain']}. The square corresponds to the case shown in Fig. \ref{['fig5']} and the star corresponds to the cases shown in Figs. \ref{['fig6']} and \ref{['fig7']}.
• Figure 2: Stability bound for the explicit method versus the anomalous exponent $\gamma$. The lines BDF1, BDF2, BDF3 and NG2 correspond to the theoretical stability bounds of the explicit method when the BDF1, BDF2, BDF3 and NG2 discretization formulae of the Riemann-Liouville derivative are used. The circles are numerical estimates of the BDF1 stability bound YusteAcedoFracExpliYusteAcedoFracExpli2. The squares correspond to the cases shown in Fig. \ref{['fig3']} and the star is the case of Fig. \ref{['fig4']}.
• Figure 3: Solution of the subdiffusion equation \ref{['subdeq']} with absorbing boundary conditions $u(0,t)=u(1,t)=0$ and initial condition $u(x,0)=x(1-x)$. The symbols correspond to the BDF1-explicit numerical solution and the lines correspond to the exact analytical solution. The solution is shown for the time $t=0.5$ for $\gamma=0.5$, $S=0.33$, $\Delta x=1/10$ (triangles), $\gamma=0.75$, $S=0.4$, $\Delta x=1/20$ (squares) and $\gamma=1$, $S=0.5$, $\Delta x=1/50$ (circles). These cases are marked by squares in Fig. \ref{['fig2']}
• Figure 4: BDF1-explicit numerical solution $u(x,t)$ for the same problem as in Fig. \ref{['fig3']} but with $\gamma=1/2$, $K_\gamma=1$, $S=0.37$, and $\Delta x=1/20$ after 150 time steps (squares) and 200 time steps (circles). The lines are plotted as a visual aide. This case corresponds to the point marked by a star in Fig. \ref{['fig2']}.
• Figure 5: Numerical solution $u(x,t)$ obtained by means of the weighted average method with weight factor $\lambda=0.8$ for the same problem as in Fig. \ref{['fig3']} but with $\gamma=1/2$, $K_\gamma=1$, $S=0.55$, and $\Delta x=1/20$ after 500 time steps (circles). The line is the exact analytical result. This case corresponds to the point marked by a square in Fig. \ref{['fig1']}.