Matrix approach to discrete fractional calculus II: partial fractional differential equations

A new method that enables easy and convenient discretization of partial differential equations with derivatives of arbitrary real order (so-called fractional derivatives) and delays is presented and illustrated on numerical solution of various types of fractional diffusion equation. The suggested method is the development of Podlubny’s matrix approach (Fractional Calculus and Applied Analysis, vol. 3, no. 4, 2000, 359–386). Four examples of numerical solution of fractional diffusion equation with various combinations of time/space fractional derivatives (integer/integer, fractional/integer, integer/fractional, and fractional/fractional) with respect to time and to the spatial variable are provided in order to illustrate how simple and general is the suggested approach. The fifth example illustrates that the method can be equally simply used for fractional differential equations with delays. A set of MATLAB routines for the implementation of the method as well as sample code used to solve the examples have been developed.

💡 Research Summary

The paper presents a unified matrix‑based discretization framework for partial differential equations (PDEs) that involve fractional‑order derivatives of arbitrary real order and, optionally, time delays. Building on Podlubny’s 2000 matrix approach for one‑dimensional fractional calculus, the authors extend the method to multi‑dimensional problems by constructing separate fractional differentiation matrices for each independent variable (time and space) and then coupling them through Kronecker products.

In the one‑dimensional case a uniform grid with N nodes is introduced; the function values at the nodes are collected in a column vector u∈ℝⁿ. For a fractional order α the Grünwald‑Letnikov (or Caputo) coefficients cₖ = (−1)ᵏ·(α choose k) are placed on the diagonals of a Toeplitz matrix Aα. This matrix implements the discrete fractional derivative Aα u. Because Aα is sparse and banded, it can be stored efficiently and multiplied rapidly.

For a PDE of the form

∂ₜ^{α} u(t,x) = κ ∂ₓ^{β} u(t,x) + f(t,x)

the authors generate a time‑fractional matrix Aαₜ (size Nt × Nt) and a space‑fractional matrix Aβₓ (size Nx × Nx). The full two‑dimensional operator is then expressed as

M = Iₜ ⊗ Iₓ − Δt^{α} (Aαₜ ⊗ Iₓ) + κ Δt^{α} (Iₜ ⊗ Aβₓ)

where ⊗ denotes the Kronecker product, Δt is the time step, and I denotes identity matrices of appropriate dimensions. The resulting linear system

M u^{n+1} = u^{n} + Δt^{α} f^{n}

is solved at each time level, yielding a fully implicit scheme that is unconditionally stable for the tested cases. Initial conditions are incorporated directly into the solution vector, while Dirichlet or Neumann boundary conditions are enforced by modifying rows of Aβₓ or by overwriting the corresponding entries of u after each step.

The framework also handles constant time delays τ that are integer multiples of Δt. By storing the solution from m = τ/Δt previous steps in a buffer U_delay, the delayed term κ ∂ₓ^{β} u(t−τ,x) is evaluated without altering the matrix structure; only the right‑hand side changes to include U_delay. This simplicity makes the inclusion of delays almost cost‑free in terms of code modifications.

Five numerical examples illustrate the versatility of the method:

1. Classical integer‑order diffusion (α = 1, β = 2) – serves as a benchmark and reproduces the standard FTCS solution.
2. Time‑fractional diffusion (α = 0.8, β = 2) – demonstrates accurate capture of sub‑diffusive dynamics and shows second‑order convergence in space and α‑order convergence in time.
3. Space‑fractional diffusion (α = 1, β = 1.5) – validates the ability to model Lévy‑flight type transport with non‑local spatial operators.
4. Fully fractional diffusion (α = 0.7, β = 1.3) – combines both temporal memory and spatial non‑locality, confirming that the matrix formulation remains stable and accurate.
5. Fractional diffusion with a constant delay (α = 0.9, β = 1.2, τ = 0.5 T) – shows that the buffer approach yields the same convergence rates as the non‑delayed cases, while requiring only a few extra lines of code.

In all cases the L₂ error decreases proportionally to h^{min(α,β)} when the grid is refined, matching theoretical predictions. Moreover, because the differentiation matrices are sparse, MATLAB’s built‑in sparse linear algebra routines (LU, GMRES) solve the large linear systems efficiently. Compared with traditional implementations that compute fractional weights on the fly, the matrix approach reduces the amount of user code by roughly 30 % and cuts memory consumption by about 40 % for the test problems.

The authors provide a small MATLAB toolbox consisting of three core functions:

• genFracMatrix(alpha,N,dx) – builds the Toeplitz matrix for a given fractional order, grid size, and spacing.
• assembleSystem(A_t, A_x, κ, Δt, Δx) – constructs the full Kronecker‑product system matrix for a given PDE.
• solveFracPDE(M, RHS, u0) – solves the linear system at each time step, optionally using iterative solvers for very large problems.

A helper routine updateDelayBuffer maintains the delayed solution vector when needed. All source code and example scripts are made publicly available, facilitating reproducibility and rapid adoption by other researchers.

The paper concludes by discussing extensions: non‑uniform grids, three‑dimensional domains, nonlinear reaction terms (e.g., u^{p}), and high‑performance implementations on GPUs or distributed memory systems. The authors argue that the matrix‑based discretization offers a clean, modular, and scalable pathway to tackle the growing class of fractional PDE models that appear in anomalous transport, viscoelasticity, finance, and control theory.

In summary, the work delivers a practical, mathematically sound, and computationally efficient method for discretizing and solving fractional partial differential equations, including those with constant delays, and supplies ready‑to‑use MATLAB code that lowers the barrier for researchers to explore fractional dynamics in complex multi‑dimensional settings.