Source-linked AI summary
Matrix approach to discrete fractional calculus II: partial fractional differential equations
Igor Podlubny, Aleksei V. Chechkin, Tomas Skovranek, YangQuan Chen, Blas M. Vinagre Jara
TL;DR
Partial fractional differential equations and delays require a numerical discretization method that handles arbitrary real-order derivatives across different time and space formulations. The paper develops a triangular-matrix approach that unifies integer and fractional operators, demonstrates it on multiple fractional diffusion equations, and reports agreement with prior numerical results.
Problem
The expanding range of fractional kinetic equations calls for a numerical method that is relatively simple while handling different time and space fractional formulations.
Method
The paper uses matrix representations of discretized fractional operators, including triangular strip matrices, to unify differentiation of arbitrary real order and n-fold integration.
Results
The method is applied to integer, fractional, and delayed diffusion equations, with the time-fractional solution in perfect agreement with a prior approach and the classical example matching prior values.
Takeaways & Limitations
The suggested method provides a unifying approach to numerical solution of partial differential equations of integer and non-integer order, including equations with delays.
Takeaways & Limitations
The paper considers one spatial variable and linear problems, although it states that extensions to multiple spatial variables and nonlinear problems are possible.
Abstract
from arXiv · showhide
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.
1. Introduction
Fractional diffusion models extend classical transport equations to systems with nonstandard temporal or spatial behavior, motivating numerical methods that are both general and simple. The paper proposes a matrix-based approach that unifies integer- and fractional-order discretization, including delayed equations.
- Motivation: Fractional kinetic equations model transport whose temporal evolution deviates from exponential Debye or Gaussian laws and fast-decaying correlations.Applications include disordered media, turbulence, underground pollution, protein dynamics, and complex reactions.
- Fractional diffusion models: Time fractionalization replaces the first-order time derivative with a Caputo derivative of order α less than 1.Setting α = 1 recovers the classical diffusion equation; Caputo derivatives allow initial conditions stated using function values and integer-order derivatives.
- Fractional diffusion models: Space fractionalization replaces the second-order spatial derivative with a derivative of order β between 1 and 2.The symmetric Riesz derivative corresponds to symmetric jump probabilities, whereas asymmetric derivatives represent preferable jump directions.
- Motivation: The growing range of fractional kinetic models makes a numerical method that is both relatively simple and sufficiently general desirable.The paper presents its approach as a response to this need.
- Related numerical methods: Existing numerical methods use substantially different discretizations for normal and fractional derivatives, including L2, L1, Grünwald-Letnikov, finite-element, spectral, Laplace-inversion, and random-walk approaches.These methods target different derivative types, domains, coefficients, and equation classes.
- Contribution: The paper proposes a matrix representation of discretized fractional operators that unifies differentiation of arbitrary, including integer, order with n-fold integration.The approach is intended to simplify numerical solution across different partial fractional equations, including equations with delays.
2. The idea of the suggested method
The method discretizes the full space-time domain simultaneously and represents fractional differentiation and integration with triangular strip matrices. Boundary and initial conditions then enter a single algebraic system assembled from the interior-node equations.
- Core representation: The method uses triangular strip matrices to discretize differentiation and integration operators of arbitrary real order.These matrices provide the basic representation underlying the suggested approach.
- Global discretization: Rather than advancing between time layers, the method considers the whole time interval at once and forms a two-dimensional space-time node net.Unknown values lie at interior nodes, while boundary values are known.
- Algebraic system: Approximating the equation at all interior nodes simultaneously produces the left-hand side of the algebraic system.Initial and boundary conditions supply terms on the right-hand side.
- Node organization: Node numbering proceeds right to left within each time level and bottom to top across time levels, although standard numberings also work.The paper adopts this ordering for clarity of presentation.
- Construction tools: The construction combines triangular strip matrices with the Kronecker product, eliminators, and shifters to approximate partial derivatives and assemble the resulting system.The subsequent sections introduce these tools before deriving the matrix equations.
- Triangular strip matrices: A triangular strip matrix is determined by its first column or first row, which can be generated from a truncated power series.Addition, subtraction, multiplication, and inversion can be expressed through operations on generating series.
4. Kronecker matrix product
The Kronecker product represents matrices as block arrays and preserves useful structural properties. Identity-matrix products place repeated copies or diagonal blocks of a matrix into larger systems.
- Definition: The Kronecker product A⊗B of an n×m matrix and a p×q matrix is an np×mq block matrix whose blocks are scalar multiples of B.The block layout repeats each entry of A as the coefficient of a corresponding B block.
- Properties: If A and B are band matrices, then A⊗B is also a band matrix.
- Properties: If A and B are lower or upper triangular, then A⊗B has the same triangular orientation.
- Special products: Left multiplication by an identity Kronecker product creates block-diagonal repetition, while right multiplication creates sparse matrices of diagonal blocks.These products support the placement of one-dimensional operators within multidimensional discretizations.
5. Eliminators
Eliminators select subsets of rows or columns from square matrices, while shifters encode discrete shifts such as delays. Together they support operator approximations and delayed differential equations.
- Eliminators: Eliminators are formed by deleting selected rows from the identity matrix.The notation S_r1,r2,...,rk records the omitted row indices.
- Eliminators: Left multiplication by an eliminator retains only matrix rows whose indices were not omitted.Right multiplication by its transpose analogously retains the corresponding columns.
- Shifters: Shifters represent discrete shifts, including delays, in approximations of differential operators and delayed differential equations.
- Shifters: The method uses upper and lower shifters with ones on diagonals above or below the main diagonal and zeroes elsewhere.
- Shifters: Shifters and eliminators shift coefficients in triangular strip matrices in opposing matrix directions.The paper gives constructions for south-west and north-east shifts.
7. Discretization of ordinary fractional derivatives
The method approximates ordinary fractional derivatives across all discretization nodes using triangular strip matrices and alternative constructions for symmetric Riesz derivatives.
- Descending node numbering swaps the original roles of B^(α)_n and F^(α)_n while preserving their backward- and forward-difference notation.
- Lower and upper triangular strip matrices simultaneously approximate left- and right-sided fractional derivatives on equidistant discretization nodes.The matrices support fractional derivatives of arbitrary order through fractional difference approximations.
- The symmetric Riesz derivative combines left- and right-sided approximations or uses Ortigueira’s centred fractional-difference approximation.
- Both symmetric Riesz approximations produce practically the same numerical results and yield a well-posed algebraic-system matrix for partial fractional differential equations.
8. Discretization of partial derivatives in time and space
Partial fractional derivatives are discretized simultaneously over a space-time grid by arranging unknowns as a column vector and applying Kronecker-product matrix operators.
- Integer-order stencils use two time layers, whereas fractional time derivatives involve all time levels from the beginning.
- With fractional derivatives in both variables, all spatial nodes at the considered time layer contribute to the approximation.
- The unknown values are stacked by descending time layers so that all grid values form one column vector for simultaneous differentiation.
- The time derivative operator is formed as the Kronecker product of the fractional matrix B^(α)_n and spatial identity matrix E_m.
- The spatial fractional derivative operator uses the Kronecker product of the time identity E_n and the symmetric Riesz matrix R^(β)_m.
- Replacing partial derivatives with their discrete analogs directly discretizes the general fractional diffusion equation and produces the matrix structure shown in Figure 7.
9. Initial and boundary conditions
The discretization assumes zero initial and boundary conditions; nonzero conditions are handled by introducing an auxiliary function.
- Nonzero initial and boundary conditions are transferred to the equation’s right-hand side through an auxiliary function satisfying homogeneous conditions.
10. Implementation in MATLAB
MATLAB routines implement the matrix approach, its fractional-derivative operators, and the demonstrations accompanying the numerical examples.
- BCRECUR computes fractional-difference coefficients, while BAN and FAN construct left- and right-sided fractional-derivative matrices.
- RANSYM and RANORT construct matrices for symmetric Riesz derivatives using two alternative formulas.
- ELIMINATOR and SHIFT implement the eliminator and shifter operations used by the matrix method.
- Demo functions implement Examples 1–5, including the delayed fractional differential-equation example.
11. Examples
The examples apply the matrix discretization method to classical, time-fractional, space-fractional, mixed fractional, and delayed diffusion equations. Across these cases, the method reproduces reference solutions and handles varying derivative orders and delays within the same framework.
- Example 1: Classical diffusion equation: The method first reproduces the classical integer-order diffusion solution in agreement with analytical and numerical results from.For h = 0.1, τ = h2/6, and n = 37, the computed u(x, t) values agree with.
- Example 2: Time-fractional diffusion equation: The time-fractional diffusion solution agrees with a numerical solution obtained by a different approach.The comparison is reported for the Caputo time-fractional formulation.
- Example 5: Fractional diffusion equation with delay: The framework extends equally simply to fractional differential equations with delays by choosing the delay as an integer multiple of the time step.With δ = kτ, k represents the number of time steps corresponding to the delay.
- Example 2: Time-fractional diffusion equation: Varying the time-fractional order α connects the fractional model to the classical case and causes y(x, t) to approach u(x, 0) as α approaches zero.At α = 1, Example 2 becomes the classical case; as α goes to zero, y(x, t) tends toward 4x(1−x) for all t.
- Example 3: Space-fractional diffusion equation: The method also discretizes diffusion equations with a Riesz spatial fractional derivative for orders 1 < β ≤ 2.Computations are shown for four different β values.
- Example 4: Fractional diffusion equation: A general diffusion equation with both time and spatial derivatives of non-integer order is handled using the same discretization approach.The examples vary non-integer α and β, with results shown for several combinations.
12. Conclusion
The suggested matrix method unifies numerical treatment of integer- and fractional-order partial differential equations, including delays. The paper demonstrates broader applicability to multidimensional, Riesz, variable-order, distributed-order, and delayed equations, while noting that its examples are linear.
- The suggested method provides a unifying numerical approach for partial differential equations with integer or non-integer orders, including delays.
- The one-dimensional formulation can be extended to two and three spatial variables using triangular strip matrices and Kronecker products.
- The presented problems are linear, although the approach can be extended to nonlinear problems.
- Example 5 presents solutions for α = 0.9, γ = 0.8, and β = 1.9 with delays δ_k = kτ for k = 6, 12, 24, 36.
- The method also applies to partial fractional equations derived from the Laplace equation by replacing second-order spatial derivatives with fractional Riesz derivatives.
- It can handle partial fractional equations with variable or distributed orders and equations with delays.
Appendix: sample evaluation of the symmetric Riesz fractional derivative
The appendix evaluates the symmetric Riesz fractional derivative of φ(x) = x(1 − x) by combining its left- and right-sided Riemann–Liouville derivatives. It presents the derivative calculations for orders 1 < β < 2.
- For φ(x) = x(1 − x) and 1 < β < 2, the appendix evaluates the left-sided Riemann–Liouville fractional derivative.
- It similarly evaluates the right-sided Riemann–Liouville fractional derivative of the same function and order range.
- The symmetric Riesz fractional derivative is then obtained from the left- and right-sided derivative results.