File Exchange

image thumbnail

Calculation Lyapunov Exponents for ODE

version (3.91 KB) by Vasiliy Govorukhin
Realization of alogrithm of chaos detection by determining Lyapunov exponents.


Updated 18 Mar 2004

View License

Lyapunov exponent calcullation for ODE-system. The alogrithm employed in this m-file for determining Lyapunov exponents was proposed in A. Wolf, J. B. Swift, H. L. Swinney, and J. A. Vastano, "Determining Lyapunov Exponents from a Time Series," Physica D, Vol. 16, pp. 285-317, 1985.

For integrating ODE system can be used any MATLAB ODE-suite methods.

This function is a part of MATDS program - toolbox for dynamical system investigation

Input parameters:
n - number of equation
rhs_ext_fcn - handle of function with right hand side of extended ODE-system. This function must include RHS of ODE-system coupled with variational equation (n items of linearized systems, see Example).
fcn_integrator - handle of ODE integrator function, for example: @ode45
tstart - start values of independent value (time t)
stept - step on t-variable for Gram-Schmidt renormalization procedure.
tend - finish value of time
ystart - start point of trajectory of ODE system.
ioutp - step of print to MATLAB main window. ioutp==0 - no print, if ioutp>0 then each ioutp-th point will be print.

Output parameters:
Texp - time values
Lexp - Lyapunov exponents to each time value.

Users have to write their own ODE functions for their specified systems and use handle of this function as rhs_ext_fcn - parameter.

Example. Lorenz system:
dx/dt = sigma*(y - x)
dy/dt = r*x - y - x*z
dz/dt = x*y - b*z

The Jacobian of system:
| -sigma sigma 0 |
J = | r-z -1 -x |
| y x -b |

Then, the variational equation has a form:
F = J*Y
where Y is a square matrix with the same dimension as J.
Corresponding m-file:
function f=lorenz_ext(t,X)
SIGMA = 10; R = 28; BETA = 8/3;
x=X(1); y=X(2); z=X(3);
Y= [X(4), X(7), X(10);
X(5), X(8), X(11);
X(6), X(9), X(12)];
R-z,-1,-x; y, x,-BETA];

% Run Lyapunov exponent calculation:

[T,Res]=lyapunov(3,@lorenz_ext,@ode45,0,0.5,200,[0 1 0],10);

See files: lyapunov.m - algorithm
lorenz_ext - Lorenz system extended rhs
run_lyap - example of calling and result visualization
Govorukhin V.N.
This file is intended for use with MATLAB and was produced for MATDS-program
lyapunov.m is free software. lyapunov.m is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY.

Cite As

Vasiliy Govorukhin (2020). Calculation Lyapunov Exponents for ODE (, MATLAB Central File Exchange. Retrieved .

Comments and Ratings (46)

Dayuan Chen

Hi, I have used this code to calculate Lyapunov exponents for several chaotic systems, however in calculating Hyperchaotic systems, I think I get the wrong Lyapunov exponent since there is only one positive result. Please guide me to work it out. Thank you.

Dayuan Chen

I'm working on Lyapunov exponents. I have used this code in my research work. Now I required to plot Lyapunov exponents over the variation of Bifurcation Parameter instead of time. Kindly guide me how can I do this by using this code.
Thank You.

fxcvxc cxvx


Fan Gao

when i use this program to calculate lyapunov exponent for Lorenz system with 4-dimensional i get the following error
??? Error using ==> odearguments at 122
Solving LORENZ_EXT requires an initial condition vector of
length 20. i need help.


I need to found the maximal lyapunov exponent from a time series, somebody has a program to do that? or cai i used this program for that porpuose



Please, what is the optimal choice for the parameter "stept" (step on t-variable for Gram-Schmidt renormalization procedure)?

I would like to know if it is possible to determine the Lypunov Coefficients with this tool


My ODE systems are: dx/dt =-2x+9/2(abs(x+1)-abs(x-1))+7y


How do you define your ODE function?


I have question,How to calculate Jac matrix in this program if i hava a pieces smooth function in right hand? Ex dx/dt =...+1/2(abs(x+1)-abs(x-1)). Thank for your help

Manel Soria

Thanks !
Note that in Wolf 85 paper, the exponent is defined using log2 and not log as in this code. Just change log by log2 in this line
for k=1:n1 cum(k)=cum(k)+log2(znorm(k)); end;
... if you want to implement the same version as in Wolf 85.
Please correct me if I'm wrong.


what about the fractional order lyapunov exponent? how can I calculate it? Thanks


Hello ,
I am trying to understand Wolf Algorithm to implement for my 13 DOF ODE system of motion equations. It will be great if anyone can help to know how i can calculate Lyapunov exponent for multi DOF system. I found matlab code for henon map, lorenz system etc.



Please, I have a time series signal, named "X". How can I calculate the Lyapunov exponent with this program ? Thank you !



this is great to learn from, might i also suggest users look at:

Pearlyn Lim


Pearlyn Lim

I hav a time series signal, which i saved under the variable "X". How do i use the cides provided to calculate the Lyapunov Exponent? Thank you

Pearlyn Lim

thank you

Behnam Molaee Ardekani

It works good.

For those interested readers who want to know more about the algorhtm of computing Lyapanov Exponent and do not have access to the refered article (i.e. A. Wolf et al 1985) I suggest them to read the following article:
K. Ramasubramanian, M.S. Sriram, A Comparative study of computation of Lyapunov spectra with different algorithms, Physica D, 139 pp 72-86, 2000.

Hint: The programmer has not used all capabilities of matrix manipulations in the MATLAB; instead he/she has used lots of "for" commands in the code.

luuk buur


alaa taha

I am working with the algoritm, with the same example, it works good
can you halpe me to use this program to work whith pendulum system and I greatly apreciate this help

alaa taha

I am working with the algoritm, with the same example, it works good
can you halpe me to use this program to work whith pendulum system and I greatly apreciate this help

Laura Villamizar

It´s excellent tool for the calculation of Lyapunov exponents.

Laura Villamizar

I was working with the algoritm, with the same example, the same conditions; I just write again the code, but it reports an error:
??? Error using ==> +
Function '+' is not defined for values of class 'function_handle'.

Error in ==> C:\MATLAB6p5\lyapunov\lyapunov.m
On line 86 ==> n1=n; n2=n1*(n1+1);

Andrew Lim

thanks for the code, just one question. The result produce by the code does not seems to correspond to the results shown in the book "Chaotic and Fractal Dynamics, Moon,1992, pg318". Why??

Javier Bejarano

José Armando Fernández Gallego

Muy buen algoritmo, ha sido muy útil para encontrar los Exponentes de Lyapunov usando otros sistemas continuos.

rui yue ouyang

Tanmoy Banerjee

Its a very useful for those who are new to the field and trying to cpture the vast field. Thank you.

Li Xian-feng

Good !I hope anybody ,who majoring bifurcation and chaos, join my group,and become friends from now on! I am a Ph. D in China!

mihailescu emilian


Xiang Ji

Excellent. I am puzzling for the calculation of lyapunov exponents for a long time. Thanks for your work.

amar bayou

good afternoon every one
I search a programm regarding lyapunov exponents calculation for a time series
thanks for your help

Hossein Nejat


Yolande Taffoti

Excellent program. Did somebody wrote one program that can calculate the Lyapunov exponent of coupled system? Like two coupled Lorenz attractors with then 6 degree of freedom?

sohel ali

I wrote a program to compute Lyapunov Exponent for Chaotic systems using equations. But I would like to know if anybody has written a program to compute Lyapunov Exponents from Time Series.

Jojo Blanza

Any reference on computing the lyapunov exponents of synchronized chaotic systems.

francesco salvadori

Excellent program!
I'm writing an upgrade for calculation of Lyapunov spectra...
I have a question: why the second Lyapunov exponent do not converge to zero ?

Haigeng Luo

I wrote a matlab programme for caculating the Lyapunov_exponents of chaotic system.Your programme give me a good reference,thank you!But there is a problem:the second lyapunov exponent do not converge to zero,as my result,why?


Good work

MATLAB Release Compatibility
Created with R12
Compatible with any release
Platform Compatibility
Windows macOS Linux