Fractional HBVMs for FDEs

• FHBVMs are spectrally accurate numerical integrators for fractional differential equations, using polynomial-based approximations and robust memory management.
• They employ Jacobi polynomial expansions and a Runge–Kutta-type formulation to achieve rigorous stability, convergence, and energy conservation.
• MATLAB implementations (fhbvm and fhbvm2) demonstrate superior performance over classical solvers by reducing computational cost and handling multi-order systems efficiently.

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 (Brugnano et al., 2024, Brugnano et al., 22 Mar 2025, Brugnano et al., 4 Dec 2025, Brugnano et al., 2024).

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 (Brugnano et al., 2024, Brugnano et al., 22 Mar 2025, Brugnano et al., 4 Dec 2025, Brugnano et al., 2024).

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 (Brugnano et al., 2024, Brugnano et al., 22 Mar 2025, Brugnano et al., 4 Dec 2025).

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 (Brugnano et al., 2024, Brugnano et al., 22 Mar 2025, Brugnano et al., 2024).

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 (Brugnano et al., 22 Mar 2025).
• 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 (Brugnano et al., 22 Mar 2025, Brugnano et al., 4 Dec 2025).

• 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 (Brugnano et al., 2024, Brugnano et al., 4 Dec 2025).
• Short Memory: Old history vectors $\gamma^{(\nu)}$ can be discarded once their effect decays, reducing memory cost (Brugnano et al., 2024).

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

Implementation is provided in MATLAB codes fhbvm (Brugnano et al., 2024, Brugnano et al., 2024) and fhbvm2 (Brugnano et al., 22 Mar 2025), 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$ (Brugnano et al., 2024, Brugnano et al., 2024, Brugnano et al., 22 Mar 2025).

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) (Brugnano et al., 4 Dec 2025).

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 (Brugnano et al., 2024, Brugnano et al., 4 Dec 2025, Brugnano et al., 2024).

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 (Brugnano et al., 2024).
• 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 (Brugnano et al., 2024, Brugnano et al., 22 Mar 2025).
• Solver Tolerance: Blended Newton/modified Newton provides robust convergence in $2$–$5$ iterations/step with default tolerance $10^{-12}$ to $10^{-14}$ (Brugnano et al., 2024).
• Short Memory Implementations: Discard very old memory coefficients in large-scale, long-time simulations (Brugnano et al., 2024).
• Extension to multi-order and large-scale systems: Multiple orthogonal polynomial quadrature and block-parallel solvers are available (Brugnano et al., 4 Dec 2025).

Future Directions include adaptive step-size and order selection, fast-memory (e.g., oblivious convolution) techniques, and further extensions to fractional-in-time PDEs (Brugnano et al., 22 Mar 2025).

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/ (Brugnano et al., 2024, Brugnano et al., 22 Mar 2025, Brugnano et al., 4 Dec 2025, Brugnano et al., 2024).