---
title: Fractional HBVMs for FDEs
url: https://www.emergentmind.com/topics/fractional-hbvms-fhbvms
type: topic
---

# Fractional HBVMs for FDEs

Fractional Hamiltonian Boundary Value Methods (FHBVMs) are a class of spectrally accurate, step-by-step numerical integrators for initial value problems (IVPs) of fractional differential equations (FDEs), generalizing the Hamiltonian Boundary Value Methods (HBVMs) to the fractional-order setting. FHBVMs admit an efficient Runge–Kutta-type structure leveraging Jacobi polynomial expansions and are implemented in MATLAB codes such as `fhbvm` and `fhbvm2`. These methods efficiently solve both single-order and multi-order FDEs, with rigorous stability and convergence properties and demonstrated superior performance over classical fractional solvers [2403.04916, 2503.17719, 2512.04737, 2407.11460].

## 1. Problem Class, Fractional Calculus, and Integral Reformulation

FHBVMs target systems of Caputo-type FDEs of the general form
\[
D_t^\alpha y(t) = f(y(t)), \quad t \in [0, T],\quad y^{(j)}(0) = y_0^j,\quad j = 0, \ldots, \ell-1,
\]
where $D_t^\alpha$ denotes the Caputo derivative of order $\alpha \in (\ell-1, \ell)$, $y \in \mathbb{R}^m$, and $f: \mathbb{R}^m \to \mathbb{R}^m$ is sufficiently smooth. The Caputo derivative admits the integral representation
\[
D_t^\alpha y(t) = \frac{1}{\Gamma(\ell - \alpha)} \int_0^t (t-x)^{\ell-1-\alpha} y^{(\ell)}(x)\, dx,
\]
with the equivalent Volterra–Riemann–Liouville integral form for the solution,
\[
y(t) = \sum_{j=0}^{\ell-1} \frac{t^j}{j!} y_0^j + \frac{1}{\Gamma(\alpha)} \int_0^t (t-x)^{\alpha-1} f(y(x))\, dx.
\]
This integral form is critical for the polynomial-based approximation central to FHBVMs [2403.04916, 2503.17719, 2512.04737, 2407.11460].

## 2. Method Derivation: Jacobi Expansions and Stepwise Discretization

On each interval $[t_{n-1}, t_n]$, the right-hand side $f(y(t))$ is projected onto an orthonormal Jacobi basis $\{P_j(c)\}$ for $c = (t - t_{n-1})/h_n \in [0, 1]$, with weight $\omega(c) = (1-c)^{\alpha-1}$:
\[
f(y(t_{n-1} + c h_n)) \approx \sum_{j=0}^{s-1} P_j(c) \gamma_j^{(n)}, \quad \gamma_j^{(n)} = \int_0^1 \omega(\tau) P_j(\tau) f(\sigma_n(\tau h_n))\, d\tau.
\]
The approximate solution $\sigma_n$ then satisfies, after inverting the fractional integral,
\[
\sigma_n(c h_n) = T_\ell(c h_n) + h_n^\alpha \sum_{j=0}^{s-1} I^\alpha P_j(c)\, \gamma_j^{(n)} + \text{history terms},
\]
where $I^\alpha P_j(c) = (1/\Gamma(\alpha)) \int_0^c (c-\tau)^{\alpha-1} P_j(\tau) d\tau$ and $T_\ell$ encodes the initial data [2403.04916, 2503.17719, 2512.04737].

To close the system, the integrals $\gamma^{(n)}_j$ are computed by Gauss–Jacobi quadrature of order $2k$ at nodes $\{c_i\}$:
\[
\gamma^{(n)}_j \approx \sum_{i=1}^k b_i P_j(c_i) f(\sigma_n(c_i h_n)),\quad j=0, \dots, s-1.
\]
The FHBVM scheme then yields a nonlinear system for $\gamma^{(n)} = [\gamma^{(n)}_0, \ldots, \gamma^{(n)}_{s-1}]$ and, consequently, for the stage values $Y_i = \sigma_n(c_i h_n)$.

This construction accommodates the memory effect inherently in FDEs by maintaining history terms, and can be naturally extended to graded meshes for handling initial weak singularity, or mixed graded–uniform meshes for solutions exhibiting both early-time singularity and long-time oscillations [2403.04916, 2503.17719, 2407.11460].

## 3. Runge–Kutta Formulation and Theoretical Properties

The FHBVM(k, s) scheme admits a Runge–Kutta-type structure:
\[
Y = \phi^{n-1} + h_n^\alpha A F(Y),\quad y_n = y_{n-1} + h_n^\alpha b^T F(Y),
\]
where $Y = [Y_1, \dots, Y_k]^T$, $F(Y) = [f(Y_1), \dots, f(Y_k)]^T$, $A$ involves the integrals of Jacobi polynomials at Gauss–Jacobi nodes, and $b = [b_1, \ldots, b_k]^T$ are the quadrature weights. The spectral expansion enables spectral (in $s$) and geometric (in $k$) accuracy.

Key properties include:
- **Convergence**: For sufficiently regular $f$, global error on a uniform mesh is $O(h^{s+\alpha-1})$, while, for graded meshes $h_i = r^{i-1} h_1$, $O(h_1^2 + h_\nu^{s+\alpha})$, with $h_\nu$ the last graded step [2503.17719].
- **Stability**: The linear stability function,
  \[
  R^{(\alpha)}_s(q) = 1 + \frac{q}{\Gamma(\alpha+1)} e^T (I - qA)^{-1} e,
  \]
  with $q = \lambda h^\alpha$, admits A-stability; the stability region $S^{(\alpha)}_s$ contains the Mittag–Leffler sector $\Lambda_\alpha = \{z: |\arg(z)| > \pi \alpha/2\}$, and as $s \to \infty$, the discrete stability region approaches the continuous one [2503.17719, 2512.04737].
- **Energy Conservation**: For Hamiltonian FDEs $y^{(\alpha)} = J \nabla H(y)$ with skew-symmetric $J$, FHBVMs inherit discrete Hamiltonian conservation for polynomial $H$ of degree $\leq s$ up to spectral error [2403.04916, 2512.04737].
- **Short Memory**: Old history vectors $\gamma^{(\nu)}$ can be discarded once their effect decays, reducing memory cost [2407.11460].

## 4. Mesh Strategies and Implementation: MATLAB Codes fhbvm, fhbvm2

Implementation is provided in MATLAB codes `fhbvm` [2403.04916, 2407.11460] and `fhbvm2` [2503.17719], both available at https://people.dimai.unifi.it/brugnano/fhbvm/. The codes perform the following steps:
- Automatically select mesh: uniform, graded ($h_n = r h_{n-1}$, $r>1$ for $\alpha<1.8$), or mixed graded–uniform (graded for early times, uniform elsewhere).
- Precompute Gauss–Jacobi nodes and weights, Jacobi polynomials, and quadrature integrals $I^\alpha P_j(c_i)$ and history integrals $J_j(x)$.
- For each step, assemble memory terms, solve the nonlinear system for $\gamma^n$ using a blended Newton iteration (or simplified Newton/traditional fixed point in multi-order case), recover stage values, and update the solution.

The codes support error estimation via mesh-doubling, memory-efficient storage, and allow the user to specify polynomial degree $s$ and quadrature order $k$ [$k \geq s$], with practical choices $s=10 \ldots 22$, $k = s$ or $s+2$ [2403.04916, 2407.11460, 2503.17719].

## 5. Multi-Order Extension: Common Nodes via Multiple Orthogonal Polynomials

For systems with components having different fractional orders $\boldsymbol{\alpha} = [\alpha_1, \ldots, \alpha_\nu]$ (multi-order FDEs), FHBVMs construct a single mesh with common quadrature nodes suitable for all orders. This is achieved using multiple orthogonal polynomials (MOPs): the monic degree-$k$ MOP $\pi_k$ satisfies vanishing moments up to $q-1$ for each Jacobi weight associated with each $\alpha_i$. The nodes (zeros of $\pi_k$) and corresponding barycentric weights yield quadrature formulas exact for all polynomials up to degree $2s-1$ in each component.

Implementation solves a system for all unknown Fourier coefficients $\gamma_{ij}^n$ via simplified Newton (for $\nu > 1$) or blended iteration ($\nu=1$), with precomputation of all relevant matrices (Jacobi and their fractional integrals) [2512.04737].

## 6. Numerical Performance and Benchmarks

Benchmarking against the MATLAB BDF2-based `flmm2` and multi-order predictor–corrector and product-integration schemes, FHBVMs (typical configuration: FHBVM(22,22) or FHBVM(22,20)) yield:
- **Spectral accuracy:** Up to 16 digits for analytic solutions and smooth $f$, typically within $0.05 – 0.3$ s for scalar/small system problems, even in stiff nonlinear and multi-order settings.
- **Reduced nodal mesh size:** Orders of magnitude fewer points compared to standard methods for comparable accuracy.
- **Robust long-time simulation:** Accurate phase portraits for nonlinear, oscillatory, or limit-cycle problems (e.g., fractional Brusselator) and stiff linear systems.
- **Superiority in efficiency:** Full machine accuracy in fractions of a second versus several seconds or minutes for flmm2 and other classical solvers, especially under mesh grading and for high-order or multi-order systems [2403.04916, 2512.04737, 2407.11460].

### Representative Timings and Accuracies

| Problem Type                | FHBVM (CPU / Accuracy) | flmm2/Other (CPU / Accuracy)   |
|-----------------------------|------------------------|--------------------------------|
| Nonlinear, low $\alpha$     | 0.1s / 16 digits       | 6–30s / 8–12 digits           |
| Stiff linear, small $\alpha$| 0.8s / 10+ digits      | 60-85s / 3–5 digits           |
| Multi-order Predator–Prey   | 22s / 12 digits        | 70–600s / 4–5 digits          |

## 7. Practical Guidelines, Extensions, and Outlook

- **Parameter Selection:** $s = 10 \ldots 22$ for spectral or near-spectral accuracy; $k \geq s$ ensures quadrature accuracy [2407.11460].
- **Mesh Policy:** Use graded mesh ($r=1.1\ldots1.5$) in presence of endpoint singularity; mixed mesh for long-time simulations with oscillations. Allow automatic selection by the code [2403.04916, 2503.17719].
- **Solver Tolerance:** Blended Newton/modified Newton provides robust convergence in $2$–$5$ iterations/step with default tolerance $10^{-12}$ to $10^{-14}$ [2403.04916].
- **Short Memory Implementations:** Discard very old memory coefficients in large-scale, long-time simulations [2407.11460].
- **Extension to multi-order and large-scale systems:** Multiple orthogonal polynomial quadrature and block-parallel solvers are available [2512.04737].

**Future Directions** include adaptive step-size and order selection, fast-memory (e.g., oblivious convolution) techniques, and further extensions to fractional-in-time PDEs [2503.17719].

FHBVMs and their variants thus provide a spectrally-accurate, memory-aware, and Hamiltonian-structure preserving framework for the numerical solution of both single-order and multi-order FDEs, with rigorous theoretical guarantees and demonstrated computational superiority across canonical and challenging fractional problems. MATLAB codes are available at https://people.dimai.unifi.it/brugnano/fhbvm/ [2403.04916, 2503.17719, 2512.04737, 2407.11460].

Source: https://www.emergentmind.com/topics/fractional-hbvms-fhbvms