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

[Moore, Eric (NIH/NIDDK) [F]]

> -----Original Message-----
> From: Barrett B [mailto:barrett.n.b at gmail.com]
> Sent: Thursday, July 10, 2014 12:54 AM
> To: scipy-user at scipy.org
> Subject: Re: [SciPy-User] Maximum Lyapunov exponent of my previous
> system
> 
> Barrett B <barrett.n.b <at> gmail.com> writes:
> 
> >
> >
> >
> > I need some help setting up the calculation of the maximum Lyapunov
> exponent of the system I was describing in my previous thread,
> "Heaviside
> function in vector form." Note that I list the actual function being
> integrated, f(X, t), in my first follow-up post on that thread.
> > Here is a brief discussion of the Lyapunov exponent:
> https://en.wikipedia.org/wiki/Lyapunov_exponent
> >
> > A simple Google search may come up with a better description.
> >
> > I tried simply introducing a fourth differential variable and
> starting to
> track it after some arbitrary time t_crit, but because of the nature of
> the
> integrator, that doesn't work. So I need help trying to calculate the
> MLE.
> >
> >
> > _______________________________________________
> > SciPy-User mailing list
> > SciPy-User <at> scipy.org
> > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
> 
> Here is the code that I'm trying:
> 
> ==============
> 
> #Lorenz equations.
> def Lorenz(X):
>     return np.array([sigma*(X[1] - X[0]), #dx/dt
>                      X[0]*(rho - X[2]) - X[1], #dy/dt
>                      X[0]*X[1] - beta*X[2]]) #dz/dt
> 
> # Calculate the maximum Lyapunov exponent.
> # Source: Clyde-Emmanuel Estorninho Meador. _Numerical Calculation
> # Lyapunov Exponents for Three-Dimensional Systems of Ordinary
> # Differential Equations_.
> if Lyap:
>     x_0 = x[i_L: ]; y_0 = y[i_L: ]; z_0 = z[i_L: ] #toss the early data
>     x_Lap = x_0 + eps_L; y_Lap = y_0 + eps_L; z_Lap = z_0 + eps_L
>     delta = np.absolute(Lorenz(np.array([x_Lap, y_Lap, z_Lap])))*t_step
>     L_f = np.log(delta/eps_L)
>     L_f_temp = L_f[0]
>     L_S = np.linalg.norm(L_f[0]) +\
>         np.linalg.norm(L_f[1]) + np.linalg.norm(L_f[2])
>     print L_S/((N - i_L)*3*t_step) #p. 33
> 
> ==============
> 
> With a timestep of 0.002 and a runtime of 90 (i.e., t = np.arange(
> 0, 90, 0.002), the last line above outputs 10.9275287. But according to
> the
> source listed above, I should be getting 0.9057--not even close to my
> result. Help--what am I doing wrong?
> 

Generally, the likely hood of a helpful response will be higher if you post a working example.

>From a few moments looking at Meador's thesis, I don't think you are implementing the algorithm he describes, but I haven't looked in detail and could be mistaken. I assume you've looked at his matlab code in the appendices? It might also be worthwhile to find a copy of his reference 22, where he claims more details are presented.

Eric