[SciPy-User] Maximum Lyapunov exponent of my previous system

[Lou Pecora]

One mistake in my pseudo code:

The line,

dx= NP.RN.random_sample((n,))/sqrt(n)  # Make random unit vector

should be replaced with,

temp= NP.RN.random_sample((n,))  # Make random unit vector
templen= sqrt(dot(temp,temp))
dx= temp/templen

Sorry.
 
-- Lou Pecora


On Saturday, July 19, 2014 1:13 PM, Lou Pecora <lou_boog2000 at yahoo.com> wrote:
 


(Hi, Andrew)

Barrett,

Here is a sample of python code that you should treat as pseudo code that shows how to code up the method I outlined.  

I have a few important things to say about this whole process.  

(1) I tried to give you a theory view of what you are doing.  The version I showed is the more rigorous form of the dissertation you refer to which uses a rough approximation to the form I used.  You should use the form (with the Jacobian) I showed you.

(2) Calculating the Lyapunov exponent is not as simple as doing something that is well established as, say, an FFT. You really want to understand the math behind it so you can properly test your code and check results.  You don't need to follow all the theorems, but at least have a good sense of what your are doing and why your are using the Jacobian, for example. I strongly recommend against using any code you get as a black box and trusting what comes out.  Know how to write your own. Read the section of the book by Chua I mentioned.  You should also check the book by Schreiber and Kantz on numerical techniques for nonlinear dynamics (I forget the actual title, but a search using author names will bring it up).

(3) Andrew Fraser's code appears to give all Lyapunov exponents.  Andy knows his stuff and it's a good example of how to do that. My code is just to find the largest Lyapunov exponent.

(4) I have not tested my code, but I think it will get you pretty far if you do the above.

The code:

import numpy as NP
import numpy.random as RN

# ---- Calc Lyapanov exponent ------------------------------------
# Input:  
#  vics= initial conditions
#  n= dimension of original ODE, e.g n=3 for Lorenz
#  T= length of time step for each ODE integration.  Test for best size to use to get 
#     dx to grow by about the 1000*initial dx value.  You can try for other
#     growth ratios besides 1000.
#  NT= number of times to integrate the system for T "secs" each time, do this enough to 
#      cover most of the attractor.
#
def LypExp_calc(vics, n, T, NT):
lnsum= 0.0   # Sum of logs of dx ratios
x= vics  # set dynamic variables to initial conditions
for i in xrange(NT):
dx= NP.RN.random_sample((n,))/sqrt(n)  # Make random unit vector
y=NP.concatenate([x,dx])  # Make 2n vector to integrate original ODEs and Variational system (dx ODE)
# Integrate the system for T "secs". vecf is the vector field of
# the original system, e.g. Lorenz ODE and the variational ODE. The combined system is 
# an ODE of 2*n dimensions
yT= YourODEintegrator(y,T, vecf)  # You choose which integrator to put here, this is just
                                  # a place holder
dxT= yT[n:2*n]   # Seperate out the evolved dx part
lendx= sqrt(NP.dot(dxT,dxT))  # Calculate length of the evolved dx
lnsum += log(lendx)  # Add logarithm of dxT length to running sum
x= yT[:n]  # Reset x to new point
LypExp= lnsum/(NT*T)   # Calculate the Lyapanov Exponents
return LypExp
# ---- Combined vector field of original ODEs and Variational ODE ------------------
# You will probably have to write this to conform to what type of function YourODEintegrator
# expects.
# You have to provide the vector field function F and its Jacobian DF
#
def vecf(y):
n=3         # Dimension of your original ODE (set for Lorenz here)
x= y[:n]    # Pull out first n dynamical variables, e.g. (x,y,z) for Lorenz
dx= y[n:2*n]   # Pull out last n dynamical variables, the perturbation
vf= F(x)    # Calculate the vector field
Dvf= NP.dot(DF(x),dx)  # Calculate the variational matrix (the Jacobian of F= DF) and 
                    # multiply time the dx vector
# Return the combined vector field
return NP.concatenate([vf, Dvf])

 
-- Lou Pecora


On Friday, July 18, 2014 6:52 PM, Barrett B <barrett.n.b at gmail.com> wrote:
 


Andrew Fraser <afraser <at> lanl.gov> writes:

> 
> I wrote code to estimate Lyapunov exponents for a book I wrote several 
> years ago.  I haven't looked at it for about a year.  Can you access 
> http://code.google.com/p/hmmds/source/browse/code/applications/synthetic/ ?  
> If so, take a look at the file LyapPlot.py.
> 

It indeed is accessible, thanks. Quick
 questions:

1. What is file f?
2. I can't import lorenz. I see the part near the top about setting up the
pathname--where should I generally search for the Lorenz package?

_______________________________________________
SciPy-User mailing list
SciPy-User at scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20140719/885656e4/attachment.html>