Improved Matlab code for Lyapunov exponents of fractional order systems

This paper introduces FO_LE, an enhanced and efficient Matlab routine that combines QR-based reorthonormalization with a quadratic LIL predictor-corrector scheme to robustly compute Lyapunov exponents for both commensurate and non-commensurate fractional-order systems modeled by Caputo derivatives.

The Gist

Imagine you are watching a pot of water boil. If you drop a single leaf in, it swirls around. If you drop a second leaf right next to it, they might swirl together for a while, but eventually, the chaotic currents will push them apart.

Lyapunov Exponents are basically a way to measure how fast those two leaves drift apart.

• If they drift apart slowly, the system is predictable (like a calm river).
• If they fly apart instantly, the system is chaotic (like a violent storm). You can't predict where the leaf will be in five minutes because a tiny change in where you dropped it makes a huge difference.

This paper is about building a better, faster, and more reliable calculator to measure that "drifting apart" speed, specifically for a special type of math called Fractional-Order Systems.

The Problem: The "Memory" of the System

Most standard math models assume the system has no memory. It only cares about what is happening right now.

• Analogy: Think of a standard car. If you take your foot off the gas, it stops accelerating immediately. It doesn't remember that you were speeding up 10 seconds ago.

Fractional-Order Systems are different. They have memory.

• Analogy: Think of a heavy piece of jelly or a memory foam mattress. If you push it, it doesn't just react to your current push; it remembers how hard you pushed it 5 seconds ago, 10 seconds ago, and so on. The past influences the present. This is common in real-world things like blood flow, battery charging, or how materials stretch.

Calculating the "drifting apart" speed (Lyapunov Exponents) for these "memory-keeping" systems is incredibly hard. The old calculators (Matlab codes) were like trying to count every single grain of sand in a beach while walking backward; they were slow and sometimes made mistakes.

The Solution: A New Toolkit

The author, Marius-F. Danca, has created a new, improved toolkit called FO_LE. Here is how it works, using simple metaphors:

1. The "Memory-Safe" Engine (The LIL Solver)

To calculate the path of the leaves, you need a powerful engine to simulate the water.

• The Old Way: The previous engines used a method called "Adams-Bashforth-Moulton." It was like a hiker taking small, careful steps but getting tired quickly because they had to look back at every single step they ever took to decide where to go next.
• The New Way (LIL): The new engine uses a technique called LIL (Lagrange Interpolation at the Last step).
  • Analogy: Imagine the hiker now has a smart GPS. Instead of looking at every single step, the GPS looks at the last few steps and predicts the curve of the path perfectly. It's faster, smoother, and handles the "memory" of the system much better. It's like upgrading from a bicycle to a high-speed train that still respects the winding roads.

2. The "Tidy-Up" Crew (QR Decomposition vs. Gram-Schmidt)

As the simulation runs, the "leaves" (mathematical vectors) get stretched and squashed. If you don't fix them, they get so long they break the computer's memory, or they get so squashed they become useless. You have to "re-normalize" them (reset them to a standard size) every few seconds.

• The Old Way (Gram-Schmidt): This was like a clumsy janitor trying to straighten a pile of tangled ropes. He would pull one rope, then the next, but often he'd accidentally mess up the first one he fixed. It was messy and prone to error.
• The New Way (QR Decomposition): This is like a professional rope-tying expert. They take the whole pile of ropes and use a special technique (QR factorization) to instantly straighten them all out perfectly without tangling the others. It's cleaner, faster, and much more reliable.

Why Does This Matter?

The author tested this new toolkit on two things:

1. A Test Problem: A math problem where we already know the exact answer. The new toolkit got the answer much closer to the truth than the old methods, and it did it faster.
2. The Rabinovich-Fabrikant System: A complex, chaotic system (like a weird, swirling fluid). The new toolkit successfully identified when the system was chaotic (leaves flying apart) and when it was stable (leaves staying together), even when the "memory" settings were changed.

The Bottom Line

This paper is essentially a software upgrade.

• Before: Calculating chaos in "memory-based" systems was slow, clunky, and sometimes inaccurate.
• Now: The author gives us a FO_LE code that is:
  • Faster: It uses a smarter engine (LIL).
  • Cleaner: It uses a better way to tidy up the math (QR).
  • Robust: It works for both simple "memory" systems and complex ones with different memory settings.

It's like giving scientists a new pair of high-definition glasses to see the chaotic behavior of the universe more clearly, helping them understand everything from how diseases spread to how new materials behave.

Technical Summary

1. Problem Statement

The paper addresses the numerical computation of Lyapunov Exponents (LEs) for fractional-order (FO) dynamical systems modeled by Caputo derivatives.

• Context: LEs are crucial for characterizing chaos, stability, and predictability in dynamical systems. While efficient algorithms exist for integer-order systems, their application to FO systems—particularly non-commensurate systems (where fractional orders differ across equations)—remains challenging.
• Limitations of Existing Methods: Previous implementations (e.g., FO_Lyapunov and FO_NC_Lyapunov) relied on:
  1. The Gram–Schmidt (GS) orthogonalization procedure, which can suffer from numerical instability over long integration times.
  2. The Adams–Bashforth–Moulton (ABM) predictor-corrector method (specifically FFT-accelerated variants like fde12) for integration, which, while fast, may have lower convergence orders compared to newer high-order schemes.
• Goal: To develop a robust, efficient, and compact Matlab routine (FO_LE) that improves numerical stability and accuracy for both commensurate and non-commensurate FO systems while preserving the "full memory" structure inherent to Caputo derivatives.

2. Methodology

The proposed approach integrates three core components into a unified algorithm:

A. Variational System Formulation

The method solves an extended system consisting of the original FO equations and their linearized variational equations.

• System: $C D_t^{\alpha_i} x_i(t) = f_i(t, x(t))$
• Variational Equations: $C D_t^{\alpha_i} \phi_{ij}(t) = \sum_{k} \frac{\partial f_i}{\partial x_k} \phi_{kj}(t)$
• This creates a system of $n_e + n_e^2$ equations, where $n_e$ is the system dimension.

B. QR-Based Reorthonormalization

Instead of the classical Gram–Schmidt procedure, the paper replaces the orthogonalization step with QR decomposition.

• Process: After integrating the extended system over a renormalization interval $h_{norm}$, the perturbation matrix $\Phi$ is decomposed as $\Phi = QR$.
• Advantages:
  • Numerical Stability: QR decomposition is generally more robust than explicit GS recursion, especially when perturbation vectors become nearly linearly dependent.
  • Efficiency: It leads to a more compact Matlab implementation.
  • Stretching Factors: The diagonal entries of the upper triangular matrix $R$ provide the stretching factors ($z_k$) used to update the cumulative LE sums ($s_k \leftarrow s_k + \log |R_{kk}|$).

C. LIL Predictor-Corrector Solver

The extended system is integrated using the LIL (Lagrange Interpolation at the Last step) scheme.

• Type: A quadratic, two-step backward interpolation method (implicit predictor-corrector).
• Properties:
  • It is a full-memory method, essential for accurately modeling the non-local nature of Caputo derivatives.
  • Under smoothness assumptions, it achieves third-order accuracy, theoretically outperforming standard ABM schemes (which are typically second-order).
  • The paper utilizes a fast implementation (LIL_nc) available for both commensurate and non-commensurate orders.

D. Impulsive Algorithm Framework

The routine FO_LE is formulated as an impulsive algorithm. It alternates between continuous-time integration (using LIL_nc) and discrete orthonormalization actions (QR) at intervals $h_{norm}$. Crucially, this framework preserves the full memory structure of the fractional system, ensuring that the reorthonormalization does not destroy the intrinsic memory effects.

3. Key Contributions

1. New Matlab Routine (FO_LE): A comprehensive code that computes finite-time LEs for both commensurate and non-commensurate FO systems.
2. Algorithmic Improvement: The substitution of Gram–Schmidt with QR decomposition significantly enhances numerical robustness.
3. High-Order Integration: The integration of the LIL scheme offers higher convergence orders and better accuracy compared to previous ABM-based implementations.
4. Unified Framework: The code handles both commensurate and non-commensurate cases within the same structure, simplifying the workflow for researchers.
5. Open Source Availability: The paper provides the code for the LIL solver and the FO_LE routine, with a fast variant (LIL_nc) available on MathWorks File Exchange.

4. Results and Validation

The paper validates the method through two primary benchmarks:

A. Benchmark with Exact Solution (Non-Commutative System)

A 2D non-commensurate system with a known exact solution (involving Mittag-Leffler functions) was used to compare LIL_nc against fde12_nc2 (FFT-accelerated ABM) and classical ABM_nc.

• Accuracy: LIL_nc produced errors of the same order of magnitude as the highly optimized fde12_nc2, though fde12_nc2 was slightly more accurate in absolute terms.
• Convergence: LIL_nc demonstrated a stable experimental convergence order of ~1.61, whereas fde12_nc2 remained below 1.57. This confirms the theoretical advantage of the LIL scheme.
• Efficiency: LIL_nc was significantly faster than the non-accelerated ABM_nc and competitive with fde12_nc2 in terms of CPU time.

B. Rabinovich–Fabrikant (RF) System

The method was applied to the fractional RF system to analyze different dynamical regimes:

• Chaotic Regime: For commensurate orders ($\alpha \approx 0.999$), the system exhibited chaos, confirmed by a positive largest LE ($\approx 0.10$).
• Stable Equilibrium: For non-commensurate orders ($\alpha = [0.6, 0.8, 0.7]$), the system converged to a stable equilibrium, confirmed by all negative LEs.
• Stability Analysis: The paper successfully applied generalized stability criteria (Matignon's criterion extended to non-commensurate systems) to verify the stability of equilibria, matching the numerical LE results.

5. Significance

• Robustness: The combination of QR decomposition and the LIL solver provides a more stable tool for studying long-term dynamics in fractional systems, where numerical errors can accumulate rapidly due to the "long memory" of the derivatives.
• Accessibility: By providing a ready-to-use Matlab routine, the paper lowers the barrier for researchers to investigate stability, bifurcations, and chaos in complex fractional-order models.
• Theoretical Consistency: The method respects the mathematical structure of Caputo derivatives (full memory), ensuring that the computed LEs are physically meaningful approximations of the asymptotic spectrum.
• Future Applications: The tool is positioned as a foundation for further studies on hidden attractors, bifurcation analysis, and control in fractional-order systems.

In conclusion, the paper presents a significant advancement in the numerical analysis of fractional-order chaos, offering a tool that is simultaneously more accurate, more stable, and computationally efficient than previous state-of-the-art implementations.