Using Lagrangian descriptors to reveal the phase space structure of dynamical systems described by fractional differential equations: Application to the Duffing oscillator

TL;DR

This work extends Lagrangian descriptors (LDs) to fractional-order dynamical systems by applying them to the unforced undamped Duffing oscillator described with $D^{\alpha}$, and by comparing backward-time schemes for different fractional derivatives. It analyzes two backward-time strategies—a nonlocal implicit inverse for Grünwald-Letnikov and a time-reversing inverse for Caputo—across orders $\alpha$ near 1 and varying integration times $\tau$, revealing how phase-space structures evolve with memory effects. The central finding is that for $\alpha$ close to 1 the LDs reproduce the classical infinity-shaped NHIM geometry, while smaller $\alpha$ yields $ au$- and method-dependent distortions; the nonlocal implicit inverse generally provides the most consistent backward dynamics, albeit with convergence challenges at very small $\alpha$. The study demonstrates the utility of LDs for qualitative geometric interpretation of fractional-order systems and points to future work on more robust backward integration and broader $\alpha$ ranges, including exploration of right-handed derivatives.

Abstract

We showcase the utility of the Lagrangian descriptors method in qualitatively understanding the underlying dynamical behavior of dynamical systems governed by fractional-order differential equations. In particular, we use the Lagrangian descriptors method to study the phase space structure of the unforced and undamped Duffing oscillator when its time evolution is governed by fractional-order differential equations. In our study, we implement two types of fractional derivatives, namely the standard Grünwald-Letnikov method, which is a finite difference approximation of the Riemann-Liouville fractional derivative, and a Grünwald-Letnikov method with a correction term that approximates the Caputo fractional derivative. While there is no issue with forward-time integrations needed for the evaluation of Lagrangian descriptors, we discuss in detail ways to perform the non-trivial task of backward-time integrations and implement two methods for this purpose: a `nonlocal implicit inverse' technique and a `time-reverse inverse' approach. We analyze the differences in the Lagrangian descriptors results due to the two backward-time integration approaches, discuss the physical significance of these differences, and eventually argue that the nonlocal implicit inverse implementation of the Grünwald-Letnikov fractional derivative manages to reveal the phase space structure of fractional-order dynamical systems correctly.

Figures (4)

• Figure 1: Initial conditions colored according to their $p$-'norm' Lagrangian descriptor value \ref{['p-norm LD definition']}, by using the color bar above each panel, of the classical Duffing oscillator, i.e., system \ref{['fractional differential equation duffing oscillator']} for $\alpha=1$. The three panels are constructed for a grid of $1000\times1000$ equidistant ICs over the intervals $x\in[-1.5,1.5]$ and $y\in[-1,1]$, and show results for an integration time of (a) $\tau=5$, (b) $\tau=10$, and (c) $\tau=20$.
• Figure 2: The phase space of the Duffing oscillator governed by the fractional differential equations \ref{['fractional differential equation duffing oscillator']} for different fractional orders $\alpha$ (indicated on the left of each row) and integration times $\tau$ (indicated on the top of each column), where each IC is colored according to its Lagrangian descriptor value (using the color bar above each panel). Each plot is created using $100\times100$ equidistant ICs over the intervals $x\in[-1.5,1.5]$ and $y\in[-1,1]$. The first column of results [panels (a), (e), and (i)] reveals the phase space structure using the Grünwald-Letnikov method \ref{['grunwald']} following the 'nonlocal implicit inverse' approach (see text for details). The remaining panels in the second, third, and fourth columns correspond to results obtained by the 'time-reversing inverse' method using the Caputo derivative. The first and second columns [panels (a), (b), (e), (f), (i), and (j)] show phase space features for integration time $\tau=5$, while the results of the third [plots (c), (g), and (k)] and the fourth column [panels (d), (h), and (l)] were obtained for $\tau=10$ and $\tau=20$ respectively. Each row corresponds to the Duffing oscillator \ref{['fractional differential equation duffing oscillator']} for different fractional orders: $\alpha=0.9999$ [top row, panels (a)--(d)], $\alpha=0.99$ [middle row, panels (e)--(h)], and $\alpha=0.98$ [bottom row, panels (i)--(l)].
• Figure 3: Plots similar to those seen in \ref{['LDs for Duffing Oscillator with fractional differential equations of order close to one']}, but for $\alpha=0.95$ [first row, panels (a)--(d)] and $\alpha=0.9$ [second row, panels (e)--(h)]. White-colored points in panels (a) and (e) correspond to ICs whose backward time evolution was not computed due to the numerical instabilities of the implemented algorithm.
• Figure 4: Initial conditions of the Duffing oscillator governed by the fractional differential equations \ref{['fractional differential equation duffing oscillator']} with $\alpha=0.95$, colored according to their Lagrangian descriptor value (using the color scale above each panel). The plots are created using an equidistant grid of $100\times100$ points over the intervals $x\in[-1.5,1.5]$ and $y\in[-1,1]$, with an integration time of $\tau=5$. The Grünwalled-Letnikov derivative \ref{['grunwald']} is used for the forward time evolution of orbits to create a panel (a), and its 'nonlocal implicit inverse' version to obtain panel (b) for the backward time evolution of ICs. The results shown in panels (c) and (d) are obtained by the forward time implementation of the Caputo derivative \ref{['eq:Caputo']} and the backward time application of its 'time-reversing inverse' approach, respectively. We superimpose the forward [panels (a) and (c)] and backward [panels (b) and (d)] time evolution of three orbits (shown as cyan curves) with ICs $(x_0,y_0)=(-0.351,0.116)$, $(x_0,y_0)=(0.525,0.343)$, and $(x_0,y_0)=(0.186,-0.599)$ (indicated by cyan points) using the respective integration methods.