PaperlayerView exact revision on arXiv

Source-linked AI summary

Matlab code for Lyapunov exponents of fractional order systems

Marius-F. Danca, Nikolay Kuznetsov

arXiv:1804.01143v1physics.comp-phnlin.CD

TL;DR

The paper addresses numerical determination of Lyapunov exponents for fractional-order systems, where finite-time local exponents and their reliability require careful computation. It proposes a Benettin-Wolf-based Matlab program and reports representative exponent calculations, while identifying numerical error and scope constraints.

  • Problem

    Computing and interpreting Lyapunov exponents for fractional-order systems requires finite-time numerical procedures whose reliability depends on the computational setting.

  • Method

    The paper adapts the Benettin-Wolf algorithm to fractional-order systems, solving perturbed equations with periodic orthonormalization in a Matlab implementation.

  • Results

    1.E-2 is the reported order of relative numerical errors for the Benettin-Wolf algorithm in the considered examples.

  • Takeaways & Limitations

    The programs support numerical computation of finite-time local Lyapunov exponents for fractional-order systems, but results should be used with precaution.

  • Takeaways & Limitations

    Fractional-order initial-value systems cannot have non-constant periodic solutions, so periodic behavior is treated numerically through asymptotically or extremely-near periodic trajectories.

Abstract

from arXiv · show

In this paper the Benettin-Wolf algorithm to determine all Lyapunov exponents for a class of fractional-order systems modeled by Caputo's derivative and the corresponding Matlab code are presented. First it is proved that the considered class of fractional-order systems admits the necessary variational system necessary to find the Lyapunov exponents. The underlying numerical method to solve the extended system of fractional order, composed of the initial value problem and the variational system, is the predictor-corrector Adams-Bashforth-Moulton for fractional differential equations. The Matlab program prints and plots the Lyapunov exponents as function of time. Also, the programs to obtain Lyapunov exponents as function of the bifurcation parameter and as function of the fractional order are described. The Matlab program for Lyapunov exponents is developed from an existing Matlab program for Lyapunov exponents of integer order. To decrease the computing time, a fast Matlab program which implements the Adams-Bashforth-Moulton method, is utilized. Four representative examples are considered.

1. Introduction

The introduction motivates Lyapunov exponents for fractional-order systems while noting both their widespread computational use and theoretical concerns. It positions the paper’s Matlab implementation within this context of finite-time local exponent computation.

  • Lyapunov exponents measure average divergence or convergence rates of nearby orbits and help assess stability and sensitive dependence on initial conditions.
  • Although some researchers question their utility and rigorous interpretation, Lyapunov exponents remain widely used in nonlinear-dynamics research and software.
  • Two common definitions use exponential growth rates of fundamental-matrix columns or singular values, which can differ in general.
  • For numerical studies, the paper considers finite-time local Lyapunov exponents computed over grids covering attractors and limited reliable integration intervals.

2. Benettin-Wolf Algorithm for LEs of FO

The paper formulates Caputo fractional-order initial-value problems and their variational equations so the Benettin-Wolf algorithm can compute finite-time local Lyapunov exponents. Numerical integration uses the predictor-corrector Adams-Bashforth-Moulton method, with an optimized implementation for computation and Matlab programs for exponent analysis.

  • The Benettin-Wolf procedure for fractional-order systems requires integrating both the original differential equations and the associated variational system.
  • Under Lipschitz continuity of f, the Caputo fractional-order initial-value problem admits a unique solution.
  • The formulation adopts a dynamical-system definition for numerical purposes, while fractional systems lack the integer-order composition property and cannot have non-constant periodic solutions.
  • The predictor-corrector Adams-Bashforth-Moulton method integrates the fractional-order equations and is applicable to fully general systems without special assumptions.
  • The method has error roughly proportional to h^2, so an error near 1.0E-6 suggests a step size close to h = 1.0E-3.
  • Because fractional-order equations have long-memory processes, the implementation uses the faster optimized FDE12.m solver and accepts order, function, time span, initial conditions, and step size.

3. The Matlab code for LEs

The Matlab implementation computes Lyapunov exponents for commensurate autonomous fractional-order systems by integrating the original and variational systems, periodically orthonormalizing perturbations, and accumulating their growth rates.

  • Assumptions: The programs assume autonomous, commensurate fractional-order systems with smooth right-hand sides and typically choose orders near 1 for significant chaos.Incommensurate systems can be treated similarly, with a different fractional-equation numerical method.
  • Algorithm: Algorithm 1 integrates the extended fractional-order system, applies Gram-Schmidt orthonormalization every normalization interval, and accumulates logarithmic vector magnitudes.The extended system combines the original equations with the variational equations.
  • Implementation: FO_Lyapunov.m adapts existing integer-order Benettin-Wolf implementations for fractional differential equations.The paper identifies Matlab, BASIC, and Fortran implementations as source variants for the adaptation.
  • Interface: The Matlab function FO_Lyapunov accepts the equation count, extended-system function, time span, normalization step, initial state, solver step, fractional order, and output interval.The extended-system file must include both the right-hand side f and the Jacobian J, and the program plots exponent evolution over time.

4. Numerical tests

The numerical tests apply the implementation to a representative fractional-order system over finite time, using an extended state that includes both the original variables and variational variables.

  • Simulation scope: The numerical simulations generally use finite intervals, with t ∈[0, 300] for the reported tests and t ∈[0, 500] for the Lorenz system.The paper notes that numerical simulations can be reliable only over relative finite-time intervals.
  • RF system: For the three-dimensional RF system, the extended calculation uses 12 variables: 3 original states and 9 variational variables.The variational matrix contributes ne×ne variables, giving ne(ne+1) total variables.
  • RF system: The RF-system example calls FO_Lyapunov with ne=3, t ∈[0, 300], initial state [0.1;0.1;0.1], and fractional order q=0.999.Intermediate exponent values are printed every 1000 normalization steps.
  • RF system: For t ∈[0, 300], the RF system has LE = (0.0749, 0.0018, -2.0850).This is the last line reported in Table 1.

2. If one considers the Lorenz system

The paper demonstrates exponent calculations across systems and parameter sweeps, including a four-dimensional example with three positive exponents and surfaces showing how exponents vary with parameters and fractional order.

  • Numerical examples: For the four-dimensional system with a=1, b=0.5, and q=0.98, the computed exponents are LE=(0.1262, 0.0846, 0.0778, -1.5244).Because three exponents are positive, the paper classifies this example as hyperchaotic.
  • Numerical examples: The four-dimensional example requires 20 variables in the extended system and a modified output format for printing four exponents.Its discontinuous non-smooth dynamics require a prior smooth approximation for integration and Jacobian construction.
  • Parameter sweeps: The programs plot Lyapunov exponents as functions of the bifurcation parameter p and fractional order q.The parameter-sweep programs modify FO_Lyapunov while retaining the extended-system function for the q sweep.
  • Computational limits: The parameter-sweep programs can be substantially accelerated with a while loop, but shared constants across parameter steps still slow computation.The paper recommends supplementary optimization for these programs.
  • Parameter sweeps: For the Chen-system surfaces over q∈[0.9,1] and p∈[20,30], S1 is positive only for some fractional orders, while S2 can approach the zero plane at relatively large q and p.The surfaces use a=35 and b=3, with q and p treated as variables.

5. Conclusion and discussion

The paper proposes a Matlab implementation of the Benettin-Wolf algorithm for finite-time local Lyapunov exponents in fractional-order systems. It emphasizes that numerical results require caution because errors depend strongly on setup choices and can reach order 1.E-2.

  • The paper proposes a Matlab program based on the Benettin-Wolf algorithm to compute finite-time local Lyapunov exponents for fractional-order systems.
  • Numerical results depend strongly on initial conditions, the integration interval, and the renormalization step size h norm.
  • 1.E-2 is the reported order of the Benettin-Wolf algorithm’s relative numerical errors in the considered examples.
  • The algorithm can produce relatively large errors for numerically periodic cycles, where the maximum Lyapunov exponent should be zero.

A. Program for LEs of FO

The fractional-order Lyapunov-exponent program integrates the extended system, repeatedly orthonormalizes variational vectors, accumulates logarithmic growth factors, and prints and plots the resulting exponents.

  • The program uses a fast predictor-corrector Adams-Bashforth-Moulton implementation, FDE12.m, to solve the fractional differential equations.
  • The main loop integrates the extended system from t to t + h_norm and updates the state with the solution at the new time.
  • The variational vectors are orthonormalized with Gram-Schmidt normalization after each integration step.
  • Running exponent sums accumulate log(zn(k)) before normalization by elapsed time.
  • The program periodically prints exponent values and plots the Lyapunov exponents against time.

B. Function LE RF.m

The LE_RF.m function organizes the variational variables as a matrix and computes their evolution through the system Jacobian.

  • The function stores variational-equation variables in a matrix X assembled from entries of the extended state vector.
  • The right-hand side of the variational equations is computed as f(4:12)=J*X.
  • The displayed implementation is intended for systems with ne≤3.

C. Program for LEs as function on p

The parameter-sweep program evaluates Lyapunov exponents across a specified range of a parameter using repeated runs of the fractional-order computation.

  • run_Lyapunov_p accepts parameter bounds p_min and p_max and a requested number n of sweep intervals.
  • The program computes the parameter increment as p_step=(p_max-p_min)/n and iterates from p_min while p<=p_max.

D. Function LE RF p.m

The LE_RF_p function implements the RF system and its variational dynamics for Lyapunov-exponent computation, with p supplied as a parameter.

  • LE_RF_p defines the RF system’s state equations using the parameter p.The function initializes the derivative vector and assigns p as the bifurcation parameter.
  • The function constructs the 3×3 variational matrix X from the extended state variables.X uses x(4) through x(12) to represent the variational matrix.
  • The Jacobian matrix J is computed explicitly from the RF state variables and p.The Jacobian supplies the linearized dynamics needed by the variational system.
  • The variational equations are evaluated as J*X and stored in the remaining components of f.This couples the RF equations with their variational system in the extended integration function.
Loading 1804.01143v1…