🔢 mathematics

Chaotic attractors and Lyapunov exponents of fractional-order systems: a tutorial on the augmented-system algorithm, with a complete Python program

This tutorial presents a complete Python implementation of the augmented-system algorithm based on diffusive representation to simulate chaotic fractional-order systems and compute their Lyapunov exponents, validating the method on linear, Lorenz, and Chen systems while explaining the underlying mathematical principles.

Original authors: Stefano Fazzino

Published 2026-10-20

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Stefano Fazzino

Original paper licensed under CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

In the study of complex systems, scientists often look for patterns in chaos. They examine how tiny differences in a starting point can grow into vastly different outcomes, a phenomenon known as the butterfly effect. To measure this sensitivity, researchers use a tool called a Lyapunov exponent. In simple terms, this number tells us how fast two nearly identical paths diverge from one another. If the number is positive, the system is chaotic and unpredictable in the long run. If it is negative, the system settles down into a steady state. For decades, this method worked perfectly for systems that change based only on their current state, like a swinging pendulum or a flowing river.

However, many real-world processes do not work that way. Materials like rubber or biological tissues often remember their past. Their current behavior depends not just on where they are right now, but on everything that happened to them before, with older events fading away slowly over time. Scientists describe this using fractional-order systems, where the rate of change is weighted by a long history. The problem is that the standard tools for measuring chaos break down here. Because these systems have no true "current state" independent of their past, the usual way of calculating how fast paths separate becomes undefined. For twenty years, researchers have simulated these systems, but they have struggled to measure their chaos accurately, often relying on shortcuts that might be measuring the wrong thing entirely.

A new tutorial by Stefano Fazzino offers a clear path forward, providing both a theoretical explanation and a complete, ready-to-use computer program to solve this problem. The core idea is to stop trying to force a fractional system into a standard box and instead treat the system's memory as its actual state. The author shows that the complex, fading memory of these systems can be broken down into a collection of simple, decaying signals. By adding these signals as extra variables to the equations, the difficult fractional system transforms into a standard, ordinary system that computers can handle easily. This new approach, called the augmented-system algorithm, allows researchers to track the system's evolution and calculate its Lyapunov exponents correctly for the first time.

The program, written in Python and designed to run on a standard laptop, simulates famous chaotic systems like the Lorenz and Chen models but with fractional-order rules. When the author tested the code on a simple linear system where the answer was already known, the program matched the exact theoretical value to within a tiny fraction of a percent. This validation proved that the method works. When applied to the more complex chaotic systems, the results revealed something surprising and fundamental about how these systems behave. In ordinary chaotic systems, the Lyapunov exponents usually include negative numbers, representing directions where paths converge and stabilize. In these fractional systems, however, the program found no negative exponents at all.

Instead of converging to negative numbers, the exponents that would normally be negative drifted slowly toward zero. This confirms a deep theoretical insight: in a system with long-term memory, disturbances do not die out exponentially fast as they do in standard systems. Instead, they fade away algebraically, following a power law that is much slower. The program also showed that if a researcher tries to ignore the distant past by cutting off the memory, the results become wildly inaccurate. The simulation demonstrated that removing even the slowest, most distant parts of the memory changes the system entirely, inflating or shrinking the calculated chaos by large margins.

The tutorial includes a complete, one-hundred-and-fifty-line program that anyone can run to see these effects. It generates visualizations of the strange, twisting shapes these systems trace out in space, known as attractors, and calculates how chaotic they are. The author warns that reading the results requires a new perspective. Because the non-chaotic directions in these systems do not settle into negative numbers, a result that looks like zero might not mean the system is stable; it might just mean the system is in a temporary, periodic phase. The program provides the full trajectory data so users can look at the actual path and distinguish between a true equilibrium and a slow, repeating cycle.

This work does more than just provide a calculator; it changes how we understand the nature of chaos in memory-dependent systems. It proves that the state of such a system is not a single point in time, but the entire history of its movements. By converting this infinite history into a manageable set of variables, the author has given scientists a reliable way to measure the unpredictability of processes ranging from viscoelastic materials to the adaptation of neurons. The findings suggest that the rules of chaos are different when the past is never truly forgotten, and the tools to measure them must be built to respect that memory.

Technical Summary: Chaotic Attractors and Lyapunov Exponents of Fractional-Order Systems

1. The Problem

Fractional-order dynamical systems, defined by the Caputo derivative of non-integer order $\nu \in (0, 1)$ , model processes where the current rate of change depends on the entire past history with a power-law weight. While these systems (such as fractional Lorenz, Rössler, and Chen systems) are known to exhibit chaos, measuring this chaos via Lyapunov exponents has remained problematic.

Standard Lyapunov algorithms rely on the concept of a "state" that determines the future evolution of a system. However, a fractional system has no such state; its future depends on its entire history. Consequently, the system lacks a flow in the traditional sense. Existing methods either truncate the memory (introducing significant errors) or apply standard algorithms to truncated schemes, computing exponents for a different system entirely. Furthermore, the mathematical definition of Lyapunov exponents for fractional systems is non-trivial because the memory of the system ages, causing the norm of the memory to grow without bound even on bounded attractors.

2. Methodology: The Augmented-System Algorithm

The paper proposes a practical solution based on the "diffusive representation," which transforms the fractional system into an infinite-dimensional ordinary differential equation (ODE) system. The core methodology rests on three pillars:

A. Diffusive Representation

The power-law memory kernel $k_\nu(\tau) = \tau^{\nu-1}/\Gamma(\nu)$ is expressed as a superposition of decaying exponentials:
$k_\nu(\tau) = \int_0^\infty \mu_\nu(\omega) e^{-\omega\tau} d\omega$
where $\mu_\nu(\omega) \propto \omega^{-\nu}$ . This allows the fractional system to be rewritten as an "augmented system" where the state is the function $\omega \mapsto z(\omega, t)$ , representing memory modes. Each mode satisfies a first-order ODE:
$\partial_t z(\omega, t) = -\omega z(\omega, t) + f(x(t))$
The original variable $x(t)$ is recovered via a weighted sum of these modes. This formulation converts the fractional system into an autonomous ODE system with an infinite-dimensional state, making it amenable to standard dynamical systems tools.

B. Discretization and Time Stepping

To implement this numerically, the continuum of modes is discretized using a quadrature rule that approximates the integral. The paper employs a specific strategy to handle the singularity of $\mu_\nu(\omega)$ at $\omega=0$ :

Singular Panel: A Gauss-Jacobi rule is used for the interval $[0, \omega_{min}]$ to capture the slow, aging modes essential for long-time behavior.
Geometric Panels: Gauss-Legendre rules in the variable $\ln(\omega)$ are used for subsequent decades up to $\omega_{max}$ .
Exponential Time Differencing (ETD): Because the system is stiff (rates $\omega$ range from $10^{-7}$ to $10^7$ ), the linear part of the mode equations is integrated analytically, while the nonlinear coupling is handled numerically using a second-order Runge-Kutta scheme.

C. Tangent System and Norm

To compute Lyapunov exponents, the linearized dynamics of the augmented system are integrated alongside the trajectory. A critical theoretical insight dictates the choice of norm for orthonormalization (via QR decomposition):

The linearized dynamics consist of a pure relaxation (norm 1) and a compact coupling.
The theory predicts that fractional systems have no negative Lyapunov exponents. Perturbations that decay do so algebraically (as a power of $t$ ), not exponentially.
Therefore, the tangent vectors must be orthonormalized in a specific weighted norm: $\|( \delta a, \delta z )\|^2 = |\delta a|^2 + \sum w_i(1+\omega_i)|\delta z_i|^2$ . Using the standard Euclidean norm would yield incorrect results.

3. Key Contributions

The paper provides a complete, self-contained Python program (approx. 150 lines) that implements this algorithm. Its specific contributions include:

A Working Implementation: A script that simulates fractional chaotic systems and computes their Lyapunov exponents, validated on linear systems with known exact solutions.
Theoretical Justification for Code Design: The paper explains why specific design choices (such as the inclusion of the singular panel and the specific weighted norm) are necessary, linking them directly to the mathematical properties of the diffusive representation.
Clarification of Output Interpretation: The paper establishes how to read the results of a fractional Lyapunov calculation:
- Positive exponents are well-defined and discrete.
- Non-positive exponents do not converge to negative numbers but drift slowly toward zero.
- A vanishing exponent ( $\lambda \approx 0$ ) requires inspection of the trajectory to distinguish between an equilibrium, a periodic window, or a transient.

4. Results

The program was validated and applied to several systems:

Linear Validation: On a linear system $CD^\nu x = \text{diag}(0.05, -1, -2)x$ with $\nu=0.8$ , the program computed the largest exponent as $0.023636$, matching the theoretical value of $0.05^{1/0.8} \approx 0.023644$ within $0.04\%$ . The other exponents correctly drifted toward zero rather than converging to $-1$ and $-2$.
Fractional Lorenz System ( $\nu=0.995$ ): The program computed a largest exponent $\lambda_1 \approx 0.856$ (compared to $0.906$ for the integer-order system). It correctly captured the trajectory's lingering near equilibria before chaos sets in.
Fractional Chen System ( $\nu=0.9$ ): The program found $\lambda_1 \approx 2.33$ .
Periodic Windows: The algorithm successfully identified a periodic window in the Chen system at $\nu=0.99$ , where the exponent drifted to zero and the trajectory became periodic after a long transient.

5. Significance and Claims

The paper positions itself as a practical companion to a recent theoretical article, aiming to resolve the ambiguity in measuring chaos in fractional systems. Its significance lies in:

Correcting Misconceptions: It demonstrates that standard Lyapunov algorithms applied to truncated memory schemes measure the wrong system.
Defining the State: It operationalizes the concept that "the state of a fractional system is its memory," allowing the application of standard dynamical systems theory.
Avoiding Truncation Errors: It explicitly warns against truncating the memory (removing slow modes), showing via kernel error analysis that such truncation fundamentally alters the system and leads to large errors in exponents.
Modest Scope: The author does not claim to solve all fractional dynamics problems but provide a robust tool for simulating and analyzing chaotic attractors and their exponents for systems of the form $CD^\nu x = f(x)$ . The paper emphasizes that the reader must understand the specific behavior of fractional exponents (the band at zero) to interpret the output correctly.