---
title: Fractional Differentiation Matrices
url: https://www.emergentmind.com/topics/fractional-differentiation-matrices
type: topic
---

# Fractional Differentiation Matrices

Fractional differentiation matrices are finite-dimensional operators designed to discretize and approximate fractional derivatives—operators of non-integer order—on grids, polynomial bases, or function spaces. These matrices underpin a wide range of numerical methods for fractional ordinary and partial differential equations, including approaches based on finite differences, Toeplitz structure, spectral collocation, and operational matrices. In recent developments, they also serve as key algorithmic tools in neural network optimization with fractional-order derivatives and as discrete analogues for more general operator powers, including matrix-order differentiation.

## 1. Classical Discretizations: Binomial, Grünwald-Letnikov, and Toeplitz Structures

Foundational approaches to fractional differentiation matrices originate from discrete analogues of the Grünwald-Letnikov formula and the properties of Toeplitz matrices.

**Discrete Binomial Approximations:**  
Let α > 0 be the fractional order and t ∈ {0,1,…,N}. Discrete Grünwald-type definitions result in lower- and upper-triangular (N+1) × (N+1) matrices with binomial weights:

- Left fractional difference (delta or nabla):
  $$
  D_{i,j}^{(\alpha)} = (-1)^{i-j}\binom{\alpha}{i-j}\mathbf{1}_{j\le i}
  $$
- Right fractional difference:
  $$
  D_{i,j}^{(\alpha)} = (-1)^{j-i}\binom{\alpha}{j-i}\mathbf{1}_{j\ge i}
  $$
where $\binom{\alpha}{k} = \alpha(\alpha-1)\cdots(\alpha-k+1)/k!$ for $k \ge 0$.

For a grid function $\vec{f} = (f_0, ..., f_N)^T$, the product $D^{(\alpha)}\vec{f}$ computes the fractional difference at each node. These matrices, and their definitions, extend to both delta and nabla operators, with practical boundary conventions (e.g., $f_j=0$ outside the domain) and direct implementation pathways [1210.0445].

**Toeplitz Matrix Discretizations:**  
Fractional differentiation and integration operators are also realized through Toeplitz matrix constructions. Given a symbol $\varphi_{α,R}(θ) = (1-Re^{iθ})^α(1+Re^{-iθ})^α$, the Toeplitz matrix $T_N(\varphi_{α,R})$ with Fourier-coefficient entries approximates (on the uniform grid) the action of the Marchaud/Riemann–Liouville derivative. The explicit asymptotic decay of the symbol’s coefficient $\varphi_{α}(n) \sim 2^α/\Gamma(-\alpha)(-1)^n|n|^{-\alpha-1}$ ensures the convergence of $N^\alpha T_N(\varphi_\alpha)\vec{f}$ to the continuous fractional derivative. Inversion of this Toeplitz matrix produces a discrete fractional integral operator corresponding to the Riemann–Liouville kernel, and all classical analytic relations (such as those for Green’s kernels and integer-order reduction) are recovered in the matrix framework [1803.10832].

## 2. Fractional Differentiation Matrices via Orthogonal Polynomials and Spectral Collocation

High-order discretizations leverage the spectral properties of orthogonal polynomials, especially Jacobi-type polynomials on Jacobi–Gauss–Lobatto (JGL) grids.

Given a polynomial interpolation basis $F_j(x)$ (Lagrange at JGL points $x_k$), any function $u(x)$ is expanded as $I_N u(x) = \sum_{j=0}^N u(x_j) F_j(x) = \sum_{k=0}^N p_k P_k^{a,b}(x)$. The construction of the fractional Caputo or Riemann-Liouville differentiation matrix involves several recurrence and expansion steps:

- Compute recurrence relations for fractional integrals of Jacobi polynomials,
- Form integral/differentiation matrices (e.g., $_C\widehat D_L^{(a,b,\alpha)}$ for left Caputo),
- Apply conversion matrices from nodal values to the polynomial basis (matrix $M$),
- Combine to yield $_C D_L^{(a,b,\alpha)} = _C\widehat D_L^{(a,b,\alpha)} M$.

These matrices, when applied, realize high-accuracy spectral collocation methods. Spectral-radius analysis shows that the interior differentiation matrices $M_i^{(\alpha,a,b)}$ scale as $N^{2\alpha}$ in spectral radius with uniform $O(1)$ proportionality constants, indicating stable and spectrally accurate discretizations. Empirical tests on Baglay–Torvic and fractional diffusion equations demonstrate spectral (exponential) convergence and surpass the accuracy of alternative operational-matrix and finite-difference schemes for moderate polynomial degree $N$ [1404.4429].

## 3. Matrix Formulations for PDEs: Banded-Strip and Operational Matrices

For time- and space-fractional PDEs, matrix discretizations generalize to Kronecker-product structures, banded matrices, and precomputable operational matrices.

**Grünwald-Type Banded Strip Matrices:**  
On a uniform time grid $t_j = j\tau$, the backward-difference approximation for a Caputo or Riemann–Liouville fractional derivative yields a banded strip matrix:
$$
B_n^{(\alpha)} = \tau^{-\alpha} \begin{pmatrix}
\omega_0 & \omega_1 & \cdots & \omega_n\\
0 & \omega_0 & \cdots & \omega_{n-1}\\
\vdots & \ddots & \ddots & \vdots\\
0 & \cdots & 0 & \omega_0\\
\end{pmatrix}
$$
with $\omega_k^{(\alpha)} = (-1)^k \binom{\alpha}{k}$. The discrete operator for a space-fractional Riesz derivative is built as a symmetric banded matrix by taking specific linear combinations or direct weights (e.g., Ortigueira’s formula) [0811.1355].

**Kronecker Product Construction:**  
To handle grids in multi-dimensions (e.g., for $u(x_i,t_j)$), the full operator is assembled as $A = B_n^{(\alpha)} \otimes I_{m+1} - \chi I_{n+1} \otimes R_m^{(\beta)}$, enabling direct encoding of time and space derivatives, with boundary and delay conditions accommodated via eliminators and shifters.

**Caputo Operational Matrices for PINNs:**  
The Caputo derivative can be discretized on arbitrary grids using telescoping auxiliary weights, producing a lower-triangular operational matrix $A$ such that for vectorized function values $\vec{f}$, the approximation $\vec{f}^{(\alpha)} \approx A\vec{f}$ holds. This supports non-uniform domains and is particularly suitable for embedding in physics-informed neural networks (PINNs), enabling matrix-vector multiplications in place of convolution quadrature steps, and admits integration with polynomial neural blocks (e.g., Legendre Block) [2401.14081].

## 4. Extensions: Matrix-Order Fractional Differentiation and Jordan Functional Calculus

Fractional calculus admits further extension where the order of differentiation itself is an $n \times n$ complex matrix $A$, with differentiations defined via matrix-valued gamma and power functions.

Consider the Riemann–Liouville matrix-order derivative for a vector function $F(t)$:
$$
D_{0+}^A F(t) = \frac{1}{\Gamma(-A)} \int_0^t (t-\tau)^{-A-I} F(\tau) d\tau,
$$
with Gamma and power functions applied via the Jordan decomposition of $A$. Discrete analogues leverage block Toeplitz matrices assembled from the matrix-valued binomial coefficients:
$$
\binom{A}{k} = \frac{A(A-I)\cdots (A-(k-1)I)}{k!}\,,
$$
enabling the construction of matrix-difference operators $D_h^A$ approximating $D_{0+}^A$. Applications involve systems of viscoelasticity, coupled anomalous diffusions, and vector-valued control and econometric models, with all analytical expressions (integral representations, inverse relations, semigroup laws) extending naturally to the matrix-order case [2005.00091].

## 5. Algebraic and Implementation Properties

**Triangular Structures and Spectral Stability:**  
- Left-fractional matrices (lower-triangular) and right-fractional matrices (upper-triangular) are well-posed for forward and backward time-stepping.
- For Toeplitz and banded matrices, spectral analysis and factorization via Hankel–Toeplitz techniques guarantee invertibility and correct analytic behavior as grid size increases [1803.10832].

**Convergence and Accuracy:**  
- Grünwald-type time discretizations yield $O(\tau)$ accuracy for $0<\alpha<1$; Riesz spatial approximations provide $O(h^2)$ accuracy.
- Spectral collocation with fractional differentiation matrices achieves exponential convergence in $N$ for smooth problems and outperforms low-order schemes [1404.4429].

**Numerical Recipes:**  
- For uniform or non-uniform grids, weights and matrix entries are constructed via binomial coefficients, spectral polynomial recurrences, or integral formulas.
- Implementation details (pseudocode) are available for a variety of scenarios—see [0811.1355], [2512.10330], [2401.14081].

## 6. Applications and Computational Impact

Fractional differentiation matrices are employed in:

| Domain                    | Matrix Approach                                  | Reference         |
|---------------------------|--------------------------------------------------|-------------------|
| Fractional ODE/PDE Solvers| Banded, Toeplitz, Spectral Collocation, Kronecker| [0811.1355], [1404.4429], [1803.10832] |
| PINNs for FDEs            | Caputo operational matrices in neural blocks     | [2401.14081]      |
| Matrix-order FDEs         | Block Toeplitz and Jordan-calculus based matrices| [2005.00091]      |
| Operator Powers           | Diagonalization, Balakrishnan’s formulae         | [2512.10330]      |
| Autograd in DNN training  | Fractional-order Jacobians (block matrices)      | [2506.07408]      |

In PINNs and deep-learning contexts, fractional-order differentiation matrices provide efficient means for integrating non-integer order constraints directly into the training loop, lead to intrinsic regularization and memory properties, and are shown empirically to accelerate or improve test set metrics for certain time-series tasks [2506.07408], [2401.14081].

## 7. Synthesis and Theoretical Unification

Fractional differentiation matrices unify and discretize all classical continuous fractional operators (Riemann–Liouville, Caputo, Weyl, Grünwald–Letnikov, Marchaud), as well as their integer-order limits, within a single matrix-based numerical and analytic framework. This extends to matrix-order and operator-function settings, spectral methods, neural operator learning architectures, and iterative solvers leveraging Kronecker and Toeplitz structures. Their properties, stability, and convergence have been analyzed across finite-difference, spectral, and functional-analytic perspectives, demonstrating the adaptability and efficacy of the matrix approach for fractional calculus in numerical and computational analysis [1803.10832], [0811.1355], [1404.4429], [2005.00091], [2401.14081], [2506.07408], [2512.10330], [1210.0445].

Source: https://www.emergentmind.com/topics/fractional-differentiation-matrices