On Runge-Kutta convolution quadrature based fractional variational integrators

Technical Summary: Runge-Kutta Convolution Quadrature Based Fractional Variational Integrators

Problem Statement
Fractional differential equations are essential for modeling systems with long-term memory, such as viscoelastic materials and anomalous diffusion. While fractional calculus extends variational principles to include dissipative systems (fractional damping), constructing high-order numerical integrators for these systems remains challenging. Previous work by [23] combined backward-differentiation convolution quadrature (BDFCQ) with Galerkin variational methods for the conservative part of the Lagrangian. However, this approach suffers from a "saturation" effect: despite BDFCQ's potential for higher-order accuracy, the resulting fractional variational integrators (FVIs) are limited to second-order accuracy. This limitation arises because the multistep structure of BDFCQ operates only on main nodal values, failing to utilize the internal stage information inherent in higher-order Galerkin discretizations of the conservative action. Consequently, there is a need for fractional variational integrators where the discretization of the fractional part is structurally compatible with the internal-stage structure of the conservative discretization.

Methodology
The authors propose a new class of Fractional Variational Integrators (FVIs) that combine Runge-Kutta Convolution Quadrature (RKCQ) with higher-order Galerkin variational methods. The methodology rests on three pillars:

1. Augmented Variational Framework: Following the restricted Hamilton principle [22], the system is modeled by doubling the state variables ($x$ and $y$) to incorporate fractional damping. The action functional is split into a conservative part ($L_{con}$) and a fractional part ($L_{fra}$).
2. Discretization Strategy:
   • Conservative Part: Approximated using higher-order Galerkin variational methods with polynomial interpolation of degree $s$ and Lobatto quadrature.
   • Fractional Part: Approximated using RKCQ based on Lobatto IIIC methods. Unlike BDFCQ, RKCQ evaluates the convolution integral at the internal Runge-Kutta stages ($t_k + c_i h$). This allows the fractional derivative approximation to align perfectly with the internal nodes used in the Galerkin discretization of the conservative action.
3. Derivation of Equations: The authors derive the Fractional Discrete Euler-Lagrange (FDEL) equations by applying the discrete restricted Hamilton principle. Key to this derivation are new discrete identities established for RKCQ:
   • A discrete semigroup property for the convolution weights (Lemma 4.1).
   • Asymmetric discrete integration-by-parts formulas (Lemmas 4.2 and 4.3), which are necessary to handle the non-local nature of the fractional derivatives within the variational framework.

Key Contributions
The paper makes four primary contributions:

1. Construction of Compatible Integrators: The authors construct FVIs where the stage structure of the RKCQ (specifically Lobatto IIIC) matches the internal nodes of the Galerkin discretization. This structural compatibility avoids the order saturation observed in previous BDFCQ-based approaches.
2. Theoretical Foundations: The paper establishes the necessary discrete identities (semigroup property and integration by parts) required to derive the FDEL equations within the RKCQ setting.
3. Convergence Analysis:
   • Theoretical: A rigorous proof of second-order convergence is provided for the 2-stage Lobatto IIIC scheme applied to a quadratic Lagrangian with $\beta = 2\alpha \in (0, 1]$.
   • Numerical: Numerical experiments demonstrate that the method achieves higher-order accuracy (4th and 6th order) for 3- and 4-stage Lobatto IIIC schemes, respectively, when applied to smooth problems.
4. Midpoint Scheme: A brief discussion and implementation of a midpoint-based fractional integrator (MIDCQ) are included, which also achieves second-order accuracy.

Results
The proposed schemes were tested on two models: a two-dimensional forced harmonic oscillator and the fractional Bagley-Torvik equation.

• Accuracy: For the 2-stage Lobatto IIIC scheme, the numerical results confirmed the theoretical second-order convergence. For 3- and 4-stage schemes, the observed convergence orders were 4 and 6, respectively, matching the theoretical expectations for the underlying Galerkin methods.
• Energy Decay: The integrators successfully reproduced the expected energy decay behavior of fractionally damped systems. In long-time simulations, the proposed FVI maintained better energy decay characteristics compared to classical RK4 methods.
• Regularity Dependence: The study highlighted that high-order convergence is contingent on the regularity of the solution. For problems with discontinuous forcing (e.g., step functions), the convergence order saturated at first order, consistent with the regularity requirements of the underlying RKCQ theory. However, for smooth forcing, the full higher-order accuracy was realized.
• Comparison with BDFCQ: The results confirm that the RKCQ approach overcomes the order limitation of the BDFCQ approach by fully utilizing internal stage information.

Significance and Claims
The paper claims that the primary significance of this work lies in the structural compatibility between the fractional derivative approximation and the conservative discretization. By utilizing RKCQ based on Lobatto IIIC, the method naturally incorporates the same stage nodes used by the Galerkin discrete Lagrangian. This allows both the conservative and fractional terms to be evaluated on a common grid, thereby avoiding the order saturation that plagues multistep-based approaches like BDFCQ.

The authors modestly note that while numerical evidence supports 4th and 6th-order convergence for higher-stage schemes, a rigorous theoretical proof for these higher orders remains an open problem. Additionally, the computational cost of the current implementation is $O(N^2)$ due to the history dependence of convolution weights, though the authors acknowledge that fast convolution techniques could potentially reduce this to $O(N \log^2 N)$ in future work. The paper concludes that this approach provides a robust framework for constructing high-order, structure-preserving integrators for fractionally damped Lagrangian systems.