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

[Barrett B]

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?