Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

A fast and memoryless numerical method for solving fractional differential equations

Nicola Guglielmi, Ernst Hairer

TL;DR

The paper tackles the computational challenge of fractional differential equations by replacing the memory-heavy nonlocal kernel with a sum of exponentials, turning the problem into a large, stiff system of ODEs. This memoryless augmentation enables the direct use of established stiff solvers such as Radau5, aided by fast structured linear algebra that exploits the resulting block structure. The authors detail the kernel-approximation method (Beylkin–Monzón), its impact on accuracy, and techniques for handling cases with α>1, including kernel splitting and extended chain tricks. Numerical experiments on scalar ODEs, the Brusselator, and 1D fractional PDEs validate the approach, showing high accuracy with reduced computational cost and memory usage, and they provide public driver codes for replication.

Abstract

The numerical solution of implicit and stiff differential equations by implicit numerical integrators has been largely investigated and there exist many excellent efficient codes available in the scientific community, as Radau5 (based on a Runge-Kutta collocation method at Radau points) and Dassl, based on backward differentiation formulas, among the others. When solving fractional ordinary differential equations (ODEs), the derivative operator is replaced by a non-local one and the fractional ODE is reformulated as a Volterra integral equation, to which these codes cannot be directly applied. This article is a follow-up of the article by the authors (Guglielmi and Hairer, SISC, 2025) for differential equations with distributed delays. The main idea is to approximate the fractional kernel $t^{α-1}/ Γ(α)$ ($α>0$) by a sum of exponential functions or by a sum of exponential functions multiplied by a monomial, and then to transform the fractional integral (of convolution type) into a set of ordinary differential equations. The augmented system is typically stiff and thus requires the use of an implicit method. It can have a very large dimension and requires a special treatment of the arising linear systems. The present work presents an algorithm for the construction of an approximation of the fractional kernel by a sum of exponential functions, and it shows how the arising linear systems in a stiff time integrator can be solved efficiently. It is explained how the code Radau5 can be used for solving fractional differential equations. Numerical experiments illustrate the accuracy and the efficiency of the proposed method. Driver examples are publicly available from the homepages of the authors.

A fast and memoryless numerical method for solving fractional differential equations

TL;DR

The paper tackles the computational challenge of fractional differential equations by replacing the memory-heavy nonlocal kernel with a sum of exponentials, turning the problem into a large, stiff system of ODEs. This memoryless augmentation enables the direct use of established stiff solvers such as Radau5, aided by fast structured linear algebra that exploits the resulting block structure. The authors detail the kernel-approximation method (Beylkin–Monzón), its impact on accuracy, and techniques for handling cases with α>1, including kernel splitting and extended chain tricks. Numerical experiments on scalar ODEs, the Brusselator, and 1D fractional PDEs validate the approach, showing high accuracy with reduced computational cost and memory usage, and they provide public driver codes for replication.

Abstract

The numerical solution of implicit and stiff differential equations by implicit numerical integrators has been largely investigated and there exist many excellent efficient codes available in the scientific community, as Radau5 (based on a Runge-Kutta collocation method at Radau points) and Dassl, based on backward differentiation formulas, among the others. When solving fractional ordinary differential equations (ODEs), the derivative operator is replaced by a non-local one and the fractional ODE is reformulated as a Volterra integral equation, to which these codes cannot be directly applied. This article is a follow-up of the article by the authors (Guglielmi and Hairer, SISC, 2025) for differential equations with distributed delays. The main idea is to approximate the fractional kernel () by a sum of exponential functions or by a sum of exponential functions multiplied by a monomial, and then to transform the fractional integral (of convolution type) into a set of ordinary differential equations. The augmented system is typically stiff and thus requires the use of an implicit method. It can have a very large dimension and requires a special treatment of the arising linear systems. The present work presents an algorithm for the construction of an approximation of the fractional kernel by a sum of exponential functions, and it shows how the arising linear systems in a stiff time integrator can be solved efficiently. It is explained how the code Radau5 can be used for solving fractional differential equations. Numerical experiments illustrate the accuracy and the efficiency of the proposed method. Driver examples are publicly available from the homepages of the authors.

Paper Structure

This paper contains 27 sections, 2 theorems, 58 equations, 1 figure, 7 tables, 1 algorithm.

Table of Contents

  1. Introduction
  2. Fractional differential equations
  3. Integral and differential fractional operators
  4. Fractional differential equations
  5. General formulation as integro-differential equations
  6. Ordinary differential equations approximating fractional problems
  7. Linear chain trick
  8. Structure of the Jacobian matrix
  9. Solving fractional differential equations by the code RADAU5
  10. Applying the code RADAU5
  11. Fast linear algebra
  12. Use of the fast linear algebra (DC_SUMEXP)
  13. Approximation of $t^{\alpha -1}$ by a sum of exponential functions
  14. Approach of Beylkin and Monzón
  15. Effect of approximating the kernel
  16. ...and 12 more sections

Key Result

Theorem 1

\newlabelthm:approx Consider the problem Volt1-m1, and assume that the approximation $\tilde{k}(t)$ of the kernel $k(t) = t^{\alpha -1}/\Gamma (\alpha )$ (with $0 < \alpha < 1$) satisfies Then, we have where $u(t)$ is solution of the scalar equation

Figures (1)

  • Figure 7.1: Solution of \ref{['example2']} with parameter $\alpha_1=1.3, \alpha_2=0.8$ in the time interval $[0,220]$ ($y_1(t)$ is in blue, $y_2(t)$ is in red) \newlabelfig:Brus

Theorems & Definitions (6)

  • Remark 1
  • Theorem 1
  • proof
  • Remark 2
  • Lemma 8.1
  • proof