Source-linked AI summary
Trapezoidal methods for fractional differential equations: theoretical and computational aspects
Roberto Garrappa
TL;DR
The paper addresses how to generalize the implicit trapezoidal method to fractional differential equations while accounting for theoretical differences, solution irregularity, and computational efficiency. It compares product integration and fractional linear multistep methods, finding method-dependent performance across fractional-order regimes and an overall efficiency advantage for fractional linear multistep methods. The analysis also identifies stability and computational limitations that constrain applicability.
Problem
Methods generalizing the same ordinary differential equation scheme can have distinct theoretical and computational behavior for fractional differential equations, requiring comparison of their strengths, weaknesses, and preferred circumstances.
Method
The paper analyzes product integration and fractional linear multistep generalizations of the implicit trapezoidal rule, including their theory, stability, implementations, and numerical behavior.
Results
The fractional linear multistep approach stabilizes discontinuous-right-hand-side solutions without unwanted spurious oscillations, while FT has the smallest error for 0 < α < 1 and NG performs best for 1 < α < 2; overall, FLMM efficiency is generally superior to PI.
Takeaways & Limitations
The preferred trapezoidal generalization depends on the fractional-order regime and method properties, with FT described as the most natural generalization based on error and stability characteristics.
Takeaways & Limitations
Despite their implicit nature, these methods can become numerically unstable for insufficiently small step sizes and therefore seem unsuited to stiff problems when 1 < α < 2.
Abstract
from arXiv · showhide
The paper describes different approaches to generalize the trapezoidal method to fractional differential equations. We analyze the main theoretical properties and we discuss computational aspects to implement efficient algorithms. Numerical experiments are provided to illustrate potential and limitations of the different methods under investigation.
1 Introduction
The paper compares product integration and fractional linear multistep approaches to generalizing the implicit trapezoidal rule for fractional differential equations. It focuses on theoretical properties, computational efficiency, and when each method is preferable, especially given stability and regularity challenges.
- Fractional differential equations model real-life phenomena across fields including bioengineering, chemistry, control, mechanics, physics, and signal processing.
- Product integration and fractional linear multistep methods derive from different approaches but coincide as the fractional order approaches the nearest integer.
- In the pure fractional case, the two approaches can have distinct theoretical and computational behavior, including different responses to discontinuous right-hand sides.
- The paper aims to describe and compare methods, identify their strengths and weaknesses, and determine circumstances favoring one method over others.
- The analysis restricts attention to product integration rules and fractional linear multistep methods generalizing the implicit second-order trapezoidal rule.
- Implicit methods receive particular attention because they have been seldom investigated, while explicit methods are often preferred despite stability issues.
2 Fractional differential equations
Fractional differential equations are formulated as initial-value problems using Caputo derivatives and equivalent integral equations. Their solutions may lack smoothness near the initial point, with mixed-power expansions and unbounded derivatives affecting numerical convergence.
- The paper considers an initial-value problem for a system of fractional differential equations with initial integer-order derivatives.
- The order parameter determines m as the smallest integer greater than α, while f is a sufficiently smooth vector field and y is the unknown solution.
- The derivative operator is introduced according to Caputo’s definition, with the formulation also connected to Riemann–Liouville integrals.
- The fractional differential equation can be equivalently written as an integral equation involving the Taylor expansion of y centered at the initial point.
- The true solution has an expansion in mixed integer and real powers, causing a lack of smoothness over the integration interval.
- The derivatives of y(t) are unbounded at the initial point, creating a major difficulty for numerical treatment.
3 Product integration rules
Product integration trapezoidal rules approximate the integral formulation by polynomial interpolation of the vector field. Uniform grids enable convolution and FFT acceleration, whereas graded grids address endpoint singularities at additional computational cost.
- 3 Product integration rules: Product integration rules evaluate the fractional integral by approximating the vector field with suitable polynomials.
- 3 Product integration rules: The product integration trapezoidal rule replaces f on each grid subinterval with a first-degree polynomial interpolant and integrates the resulting expression exactly.
- 3 Product integration rules: Under sufficient smoothness, the scheme is second-order convergent, but unbounded initial derivatives reduce the actual convergence rate below order 2.
- 3 Product integration rules: For 0 < α < 1, the error away from the initial point decreases proportionally to h1+α under ordinary smoothness assumptions, while α > 1 yields order 2 convergence.
- 3 Product integration rules: On uniform grids, convolution weights depend only on n − j, making evaluation less expensive and enabling FFT-based implementation.
- 3.2 Graded meshes: Graded grids cluster points near the left endpoint to counteract solution irregularity and can restore second-order convergence for suitable linear problems.
- 3.2 Graded meshes: Nonuniform product integration is more demanding because weights must be recomputed and the lack of convolution prevents fast algorithms.
4 Fractional linear multistep methods
Fractional linear multistep methods (FLMMs) extend classical ODE multistep formulas by generating fractional convolution weights from powers of their generating functions. The section develops trapezoidal generalizations, convergence results, and computational procedures for evaluating their weights.
- FLMM framework: FLMMs approximate the Riemann–Liouville integral using convolution quadrature on uniform grids, with convolution and starting weights independent of the step size.Starting weights address the potentially singular integrand near the initial time.
- FLMM framework: The convolution weights characterize each FLMM and are obtained as coefficients of a fractional power of the generating function of an underlying ODE linear multistep method.The generating function is constructed from the characteristic polynomials of the underlying method.
- Convergence: A stable, consistent implicit linear multistep method of order p yields a convergent FLMM of order p under the stated root condition.The zeros of the second characteristic polynomial must have absolute value at most 1.
- Computational aspects: The Miller formula evaluates the first N fractional-series coefficients in O(N^2) operations generally, while special factors can reduce the work to 2N multiplications and N additions.For the fractional trapezoidal weights, factoring the generating function and using FFT multiplication requires O(N log_2 N) work.
- Trapezoidal generalizations: The fractional trapezoidal rule is obtained from the classical trapezoidal generating function with ρ(z)=z−1 and σ(z)=(z+1)/2 by taking its α-power.The trapezoidal rule is a one-step Adams–Moulton method and belongs to the broader Adams–Moulton family.
- Trapezoidal generalizations: Newton–Gregory methods reverse the order of fractional powering and series truncation, producing a distinct family from the fractional trapezoidal construction.The interchange between truncation and α-power is explicitly identified as the source of the difference.
- Additional second-order method: The second-order BDF method is included because it has the same convergence order as trapezoidal methods and offers strong stability properties, although it is not derived from trapezoidal integration.Its fractional weights can still be computed with an overall cost proportional to N.
5 Linear stability analysis
Linear stability is characterized through the fractional stability sector Σα and the corresponding numerical stability regions. The analysis shows that stability behavior depends strongly on the fractional order and on the chosen trapezoidal generalization.
- Stability criterion: For the scalar fractional test equation, the steady state y=0 is stable if and only if λ belongs to the sector Σα.This sector generalizes the classical ODE stability condition.
- Stability criterion: When α=1, Σα coincides with the left complex half-plane C−, recovering the classical ODE stability sector.Figure 1 illustrates Σα for representative fractional orders below and above 1.
- Numerical stability: For convolution quadrature, the stability region Sα is determined under the sequence assumption stated in Theorem 6.The paper verifies that the investigated PI rule and FLMMs satisfy the theorem’s assumptions.
- Numerical stability: A(απ/2)-stability requires the entire sector Σα to lie within the numerical stability region, making it the relevant fractional analogue of A-stability.Classical A-stability requires inclusion of the whole complex left half-plane.
- Orders 0<α<1: For 0<α<1, FT preserves the exact-sector property, while NG and PI have slightly larger regions and FBDF has the largest stability region.The FT region coincides with Σα, which corresponds to convergence toward the steady state in the classical case.
- Orders 1<α<2: For 1<α<2, FT still coincides with Σα and FBDF remains strongest, but NG and PI lose A(απ/2)-stability.The loss means these implicit methods can become unstable at insufficiently small step sizes and appear unsuitable for stiff problems in this range.
6 Implementation details and Matlab codes
Efficient implementation of the fractional methods requires managing long convolution lags, starting weights, nonlinear solves, and initialization. FFT acceleration reduces the dominant convolution cost, while small fractional orders create conditioning, memory, and step-size constraints.
- Efficient convolution: FFT implementation reduces convolution effort from proportional to N^2 to N log(N)^2, avoiding impracticable growth as grid points increase.The lag term is a major bottleneck because direct integration over N steps can require O(N^2) operations.
- Starting weights: Starting weights enforce exactness for prescribed power functions to address solution singularities near the integration interval's left endpoint.For order p = 2, the weights are evaluated by solving s + 1 linear equations at each step.
- Starting weights: The starting-weight coefficient matrix is mildly ill-conditioned, but its small size usually permits analytical inversion; very small α is the main exception.Its size is s = ⌈1/α⌉ + 1, typically 3 or 4, while very small α can make the system substantially larger.
- Initialization: Initialization computes y1 through ys together in one implicit nonlinear system of size s · q, which can cause memory problems for small α and large systems.Only y0 is supplied by the problem, so the remaining initial approximations must be generated computationally.
- Nonlinear solves: Newton iterations preserve the implicit methods' stability properties, but convergence can depend on h^α, requiring small h when α is low.This dependence is harmless for 1 < α < 2 but can create unexpected difficulties for small fractional orders.
7 Numerical experiments
The experiments compare accuracy and efficiency across PI rules and FLMM trapezoidal generalizations on linear and nonlinear fractional problems. They confirm order reduction for uniform-grid PI when α < 1, show a roles swap between FT and NG across α = 1, and identify method- and grid-dependent efficiency trade-offs.
- Linear FDE experiments: For α = 0.5, uniform-grid PI converges with order 1 + α rather than 2, whereas graded-grid PI achieves full order 2.These results confirm the theoretical order predictions for the linear test problem.
- Linear FDE experiments: FT produces the smallest error for 0 < α < 1, while NG produces the lowest error for 1 < α < 2, indicating a roles swap across α = 1.The stability and accuracy trade-off helps explain why the preferred method changes with the fractional order.
- Grid effects: Graded grids perform very poorly for low α because the exponent r = 2α makes them coarse near the right endpoint, and they add no improvement for α > 1.The reported utility of graded grids is therefore restricted to the intermediate range indicated by the experiments.
- Efficiency comparisons: For 0 < α < 1, FT, NG, and FBDF are more efficient than PI; PI_U loses accuracy, while PI_G lacks convolution structure for FFT acceleration.FBDF is less efficient than FT and NG because its weights are more complex to evaluate.
- Efficiency comparisons: For 1 < α < 2, uniform-grid PI becomes more efficient without order reduction but remains less efficient than NG, which also has lower errors.The comparison concerns the error-versus-execution-time trade-off.
- Brusselator experiment: In the nonlinear Brusselator test with α = 0.8, PI_U gives good results, but FT still has the smallest error and PI_U is no more efficient than FT or NG.This provides a nonlinear example where favorable PI_U accuracy does not translate into superior efficiency.
8 Concluding remarks
The paper investigates implicit second-order fractional methods and finds that stability and efficiency depend strongly on the fractional order. Overall, FLMMs are more efficient than PI rules, while stability advantages vary across methods.
- The study discusses implicit second-order methods for fractional differential equations and their efficient implementation.
- The BDF-derived method has the largest stability regions, as in the corresponding ODE case.
- FT gives the smallest error for 0 < α < 1, whereas NG performs best for 1 < α < 2.
- Overall, FLMMs are more efficient than PI rules, despite PI methods receiving greater attention in the literature.