arxiv: 2604.09334 · v1 · submitted 2026-04-10 · 🌊 nlin.CD

Recognition: unknown

Improved Matlab code for Lyapunov exponents of fractional order systems

Authors on Pith no claims yet

Pith reviewed 2026-05-10 16:44 UTC · model grok-4.3

classification 🌊 nlin.CD

keywords Lyapunov exponentsfractional-order systemsMatlab codeCaputo derivativenumerical integrationchaosstability analysis

0 comments

The pith

An improved Matlab code computes Lyapunov exponents for fractional-order systems using QR reorthonormalization and a quadratic LIL integrator.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

This paper presents an enhanced Matlab routine called FO_LE for calculating Lyapunov exponents in systems defined by Caputo fractional derivatives. It upgrades earlier codes by adopting QR-based reorthonormalization in place of Gram-Schmidt and integrating the variational equations with the quadratic LIL predictor-corrector method. The routine works for both commensurate and non-commensurate orders while preserving the full memory structure required by the underlying model. The code is supplied directly, along with tests on a benchmark problem and the Rabinovich-Fabrikant system.

Core claim

Replacing the Gram-Schmidt procedure with QR reorthonormalization and switching to the quadratic LIL predictor-corrector scheme produces a single compact Matlab function that computes Lyapunov exponents for both commensurate and non-commensurate fractional-order systems while retaining the full-memory Caputo structure.

What carries the argument

The quadratic LIL predictor-corrector scheme together with QR reorthonormalization applied to the extended variational system in full-memory form.

If this is right

The same code handles both commensurate and non-commensurate fractional orders without separate implementations.
Higher-order integration from the quadratic LIL scheme improves accuracy relative to earlier predictor-corrector versions.
Full-memory retention ensures consistency with the original Caputo model for long-time simulations.
Benchmark comparisons confirm the method's performance against standard ABM-type integrators.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

The supplied code could serve as a base for adding other fractional derivatives such as Riemann-Liouville.
Efficiency gains may allow repeated computation of exponents inside parameter-optimization loops.
The QR approach might reduce sensitivity to round-off in high-dimensional fractional systems.

Load-bearing premise

The quadratic LIL predictor-corrector scheme combined with QR reorthonormalization computes the Lyapunov exponents accurately without introducing large numerical errors or instabilities.

What would settle it

Running the code on a fractional system with independently known exact Lyapunov exponents and finding deviations larger than the method's expected truncation error would falsify the accuracy claim.

Figures

Figures reproduced from arXiv: 2604.09334 by Marius-F. Danca.

**Figure 2.** Figure 2: Phase plot of the numerical solutions of obtained with thhe [PITH_FULL_IMAGE:figures/full_fig_p021_2.png] view at source ↗

read the original abstract

This paper presents an improved Matlab routine, FO_LE, for the numerical computation of Lyapunov exponents of fractional-order systems modeled by Caputo's derivative. It is conceived as an enhanced version of the former FO_Lyapunov and FO_NC_Lyapunov codes for commensurate and non-commensurate orders, respectively. The proposed approach replaces the Gram-Schmidt orthogonalization procedure with QR-based reorthonormalization and uses the new quadratic LIL predictor-corrector scheme for the integration of the extended variational system. Compared with the former implementations, the present routine benefits from the higher order of the fractional integrator LIL and applies to both commensurate and non-commensurate models. Like the previous code, FO_LE retains the full memory structure of the underlying Caputo model. The Matlab code for the LIL solver and for the computation of Lyapunov exponents with FO_LE are provided, while a fast implementation of LIL for commensurate and non-commensurate orders, LIL_nc, is available on MathWorks File Exchange. A benchmark problem with exact solution is used to compare the LIL-based solver with ABM-type methods, whereas the Rabinovich-Fabrikant system illustrates the computation of Lyapunov exponents in different dynamical regimes. The results indicate that the proposed implementation is a compact, robust, and efficient tool for the numerical study of stability and chaos in fractional-order systems.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Desk Editor's Note private letter to a colleague

This is a modest code update for Lyapunov exponents in fractional systems that adds QR reorthonormalization and quadratic LIL but rests on thin validation for the actual exponent values.

read the letter

This paper is mainly an incremental update to the author's earlier Matlab routines for computing Lyapunov exponents in Caputo fractional-order systems. The core changes are swapping Gram-Schmidt for QR-based reorthonormalization and using a quadratic LIL predictor-corrector to integrate the variational equations. The new routine now handles both commensurate and non-commensurate orders while keeping the full-memory structure. The code is released openly, and a separate fast LIL implementation is pointed to on File Exchange. A benchmark problem with an exact solution tests the integrator against ABM methods, and the Rabinovich-Fabrikant system shows LE results across different regimes. These steps make the tool more general and potentially more accurate than the prior versions. The code sharing itself is useful for anyone who needs a ready Matlab implementation in this area. The main limitation is the lack of quantitative checks on the Lyapunov exponents. The comparisons focus on the solver accuracy for the benchmark, but there are no reported error bounds, step-size convergence tests, or direct matches of the computed exponents to known literature values, especially in the non-commensurate case. This leaves the claim of robustness somewhat under-supported for general use. The work is aimed at researchers in fractional chaos who already work in Matlab and want an updated tool. A reader in that narrow group can get practical value from the code and the example. The paper shows straightforward engagement with the methods and does not contain obvious contradictions. It is worth sending for peer review so that specialists can look at the implementation details and suggest stronger tests for the exponent accuracy.

Referee Report

2 major / 2 minor

Summary. The manuscript presents an improved MATLAB routine FO_LE for computing Lyapunov exponents of fractional-order Caputo systems. It upgrades earlier codes (FO_Lyapunov and FO_NC_Lyapunov) by replacing Gram-Schmidt with QR reorthonormalization and adopting a quadratic LIL predictor-corrector integrator for the extended variational system. The implementation supports both commensurate and non-commensurate orders while preserving the full-memory structure, supplies the source code, and demonstrates the solver via a benchmark problem possessing an exact solution (comparing LIL against ABM-type methods) together with an illustrative computation on the Rabinovich-Fabrikant system across different dynamical regimes.

Significance. If the numerical accuracy of the combined LIL-QR procedure is confirmed, the work supplies a compact, higher-order, open-source tool that can facilitate reproducible studies of stability and chaos in fractional-order systems. The explicit provision of both the LIL solver and the FO_LE routine, plus the fast LIL_nc variant on MathWorks File Exchange, constitutes a concrete contribution to reproducibility in the field.

major comments (2)

[Abstract] Abstract and benchmark description: the claim that FO_LE is 'robust' rests on the quadratic LIL predictor-corrector plus QR reorthonormalization producing accurate Lyapunov exponents, yet only the underlying LIL integrator is compared against ABM methods on a benchmark with exact solution; no error norms, step-size convergence rates, or direct LE-value comparisons for the full variational integration are reported.
[Rabinovich-Fabrikant example] Rabinovich-Fabrikant illustration: the manuscript states that the code computes Lyapunov exponents 'in different dynamical regimes' but supplies no tabulated LE values, no comparison against literature results for the same parameter sets, and no discussion of possible instabilities arising from the full-memory implementation in non-commensurate cases.

minor comments (2)

The manuscript should explicitly state the fractional orders and parameter values used in the Rabinovich-Fabrikant demonstration so that readers can reproduce the reported regimes.
Clarify whether the supplied FO_LE code includes built-in options for varying the memory length or step-size adaptation, as these choices directly affect practical use.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the constructive feedback on our manuscript. We address each major comment point by point below, indicating the revisions we intend to make to strengthen the presentation of our results and code.

read point-by-point responses

Referee: [Abstract] Abstract and benchmark description: the claim that FO_LE is 'robust' rests on the quadratic LIL predictor-corrector plus QR reorthonormalization producing accurate Lyapunov exponents, yet only the underlying LIL integrator is compared against ABM methods on a benchmark with exact solution; no error norms, step-size convergence rates, or direct LE-value comparisons for the full variational integration are reported.

Authors: We agree that the robustness claim would be better supported by explicit validation of the full LE computation. The benchmark validates the LIL integrator on the extended variational system (which incorporates the linearized equations), and QR reorthonormalization is a standard, numerically stable procedure whose accuracy properties are well-established in the integer-order LE literature. Nevertheless, we will revise the benchmark section to include error norms, step-size convergence rates for the computed Lyapunov exponents, and direct LE-value comparisons against the exact solution where feasible. We will also adjust the abstract wording to reflect these additions. revision: yes
Referee: [Rabinovich-Fabrikant example] Rabinovich-Fabrikant illustration: the manuscript states that the code computes Lyapunov exponents 'in different dynamical regimes' but supplies no tabulated LE values, no comparison against literature results for the same parameter sets, and no discussion of possible instabilities arising from the full-memory implementation in non-commensurate cases.

Authors: We acknowledge that tabulated values, literature comparisons, and discussion of non-commensurate full-memory effects are currently absent. We will add a table of computed LE values for the reported parameter sets in the Rabinovich-Fabrikant example, include comparisons to published results for commensurate cases where they exist, and insert a concise paragraph addressing the full-memory structure for non-commensurate orders, including any observed or potential numerical considerations. revision: yes

Circularity Check

0 steps flagged

No circularity: numerical code paper with external benchmarks

full rationale

The paper presents an improved numerical implementation (FO_LE) of Lyapunov exponent computation for Caputo fractional-order systems, replacing Gram-Schmidt with QR reorthonormalization and using a quadratic LIL predictor-corrector integrator. It supplies the Matlab code, a benchmark problem with exact solution for comparing LIL against ABM methods, and an illustrative example on the Rabinovich-Fabrikant system. No theoretical derivations, first-principles predictions, or fitted parameters are claimed; the central claim of being a compact, robust tool rests on provided code and numerical illustrations rather than any self-referential reduction. Prior self-citations to FO_Lyapunov/FO_NC_Lyapunov are contextual only and not load-bearing for validity of the new results.

Axiom & Free-Parameter Ledger

0 free parameters · 2 axioms · 0 invented entities

The work relies on established fractional calculus and numerical linear algebra techniques without introducing new physical entities or ad-hoc assumptions beyond the method choice.

axioms (2)

standard math The Caputo fractional derivative is used to model the systems.
This is the standard definition for fractional-order systems in the paper.
domain assumption The variational system can be integrated alongside the original system to compute Lyapunov exponents.
Standard approach in Lyapunov exponent calculation for dynamical systems.

pith-pipeline@v0.9.0 · 5536 in / 1352 out tokens · 38129 ms · 2026-05-10T16:44:16.243212+00:00 · methodology

discussion (0)

Reference graph

Works this paper leans on

27 extracted references

[1]

https://www.mathworks.com/matlabcentral/fileexchange/122377-matlab-code-for-les-of-non-commensurate-fo
[2]

https://www.mathworks.com/matlabcentral/fileexchange/183396-lil
[3]

https://www.mathworks.com/matlabcentral/fileexchange/183478-lil_nc
[4]

https://www.mathworks.com/matlabcentral/fileexchange/my-file-exchange?s_tid=gn_mlc_fx_myfx
[5]

& Gollub, P.G.[1990)] Chaotic Dynamics - An Introduction Chaotic Dynamics- An Introduction, Cambridge University Press, New York

Baker, G.L. & Gollub, P.G.[1990)] Chaotic Dynamics - An Introduction Chaotic Dynamics- An Introduction, Cambridge University Press, New York

1990
[6]

& Strelcyn, J.-M

Benettin, G., Galgani, L., Giorgilli, A. & Strelcyn, J.-M. [1980].`` L yapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems. A method for computing all of them. P art II : Numerical application,'' Meccanica , 15(1), 21--30

1980
[7]

[1967] ``Linear models of dissipation whose Q is almost frequency independent- II ,'' Geophys

Caputo, M. [1967] ``Linear models of dissipation whose Q is almost frequency independent- II ,'' Geophys. J. Roy. Astr. S. , 13(5), 529--539

1967
[8]

& Rugh, H

Christiansen, F . & Rugh, H. H. [1997] ``Computing Lyapunov spectra with continuous Gram - Schmidt orthonormalization,'', Nonlinearity, 10(5), 1063

1997
[9]

& Fe c kan, M

Danca, M.-F. & Fe c kan, M. [2024] ``Memory principle of the MATLAB code for Lyapunov exponents of fractional-order systems'', International Journal of Bifurcation and Chaos, 2450156

2024
[10]

[2026] ``A High Order Method for Caputo Fractional Differential Equations'', submitted

Danca, M.-F. [2026] ``A High Order Method for Caputo Fractional Differential Equations'', submitted

2026
[11]

& Kuznetsov, N

Danca, M.-F. & Kuznetsov, N. [2018] ''Matlab code for Lyapunov exponents of fractional-order systems'', International Journal of Bifurcation and Chaos, 28(5), 1850067

2018
[12]

[2021] ''Matlab code for Lyapunov exponents of fractional-order systems, Part II: The non-commensurate case,'' International Journal of Bifurcation and Chaos, 31(12), 2150187, 2021

Danca, M.-F. [2021] ''Matlab code for Lyapunov exponents of fractional-order systems, Part II: The non-commensurate case,'' International Journal of Bifurcation and Chaos, 31(12), 2150187, 2021

2021
[13]

& Kuznetsov, N

Danca, M.-F., Bourke, P. & Kuznetsov, N. [2019] ``Graphical structure of attraction basins of hidden attractors: the Rabinovich--Fabrikant system'', International Journal of Bifurcation and Chaos, 29(1), 1930001

2019
[14]

& Ford, N

Diethelm, K. & Ford, N. J. [2004] ``Multi-order fractional differential equations and their numerical solution'', Applied Mathematics and Computation, 154, 3, 621-640

2004
[15]

& Ford, N.J

Diethelm, K. & Ford, N.J. [2002] ``Analysis of fractional differential equations,'' J. Math. Anal. Appl. , 265(2),229--248

2002
[16]

& Freed, A.D

Diethelm, K., Ford, N.J. & Freed, A.D. [2002] ``A predictor-corrector approach for the numerical solution of fractional differential equations,'' Nonlinear Dynam. , 29(1), 3--22

2002
[17]

Eckmann, J.-P.& Ruelle D.[1985] ``Ergodic theory of chaos and strange attractors,'' Rev. Mod. Phys. , 57(3), 617--656

1985
[18]

[2018] ``Numerical Solution of Fractional Differential Equations: A Survey and a Software Tutorial'', Mathematics , 6(2), 16

Garrappa, R. [2018] ``Numerical Solution of Fractional Differential Equations: A Survey and a Software Tutorial'', Mathematics , 6(2), 16

2018
[19]

& Mainardi F.[1997] Fractional Calculus

Gorenflo R. & Mainardi F.[1997] Fractional Calculus. In: Carpinteri A., Mainardi F. (eds) Fractals and Fractional Calculus in Continuum Mechanics. International Centre for Mechanical Sciences (Courses and Lectures), vol 378. Springer, Vienna

1997
[20]

& Trujillo, J.J.[2001]`` Differential equations of fractional order

Kilbas, A.A. & Trujillo, J.J.[2001]`` Differential equations of fractional order. Methods, results and problems, I,'' Appl. Anal., 78(1-2): 153-192

2001
[21]

Spanier, J

Oldham, K.B. Spanier, J. [1974] The Fractional Calculus: Theory and Applications of Differentiation and Integration to Arbitrary Order . Academic Press, New York

1974
[22]

[2002] ``Geometric and physical interpretation of fractional integration and fractional differentiation,'' Frac

Podlubny, I. [2002] ``Geometric and physical interpretation of fractional integration and fractional differentiation,'' Frac. Calc. Appl. Anal. , 5(4):367--386

2002
[23]

Podlubny, I.[1999] Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of their Solution and some of their Applications, Academic Press, Dan Diego

1999
[24]

& Meador C

Sarra, S.A. & Meador C. [2011]. ``On the numerical solution of chaotic dynamical systems using extend precision floating point arithmetic and very high order numerical methods,'' Nonlinear Analysis: Modelling and Control , 16(3), 340--352

2011
[25]

& Nagashima, T

Shimada, I. & Nagashima, T. [1979], ``A Numerical Approach to Ergodic Problem of Dissipative Dynamical Systems,'' Prog. Theor. Phys. 61(6) 1605--1616

1979
[26]

& Haeri, M

Tavazoei, M.S. & Haeri, M. [2009] ``A proof for non existence of periodic solutions in time invariant fractional order systems,'' Automatica, 45(8), 1886--1890

2009
[27]

Wang, P., Li, J., & Li, Q [2012] ``Computational uncertainty and the application of a high-performance multiple precision scheme to obtaining the correct reference solution of L orenz equations,'' Numerical Algorithms , 59(1), 147--159

2012