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

[Barrett B]

Moore, Eric (NIH/NIDDK) [F] <eric.moore2 <at> nih.gov> writes:

>
> 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
> 

I thought I did post an example? Here, here's some more code. This simulates
the standard Lorenz system (with parameters rho=28, sigma=10, and beta =
8.0/3). I tested the numerical simulation part, and the graph suggests that
the simulator itself works fine. It is the Lyapunov exponent calculation
that somehow fails, and I cannot figure out why.

==================================

rates = [-0.1, -0.4, -1]
sigma = 10; beta = 8.0/3; rho = 28
t_step = 0.001; t_Fin = 90; N = int(t_Fin/t_step)

#Whether or not to simulate or to calculate max Lyapunov exponent.
Lyap = True; eps_L = 0.001; t_track = 15; i_L = int(t_track/t_step)

#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

#ODE.
def f(X, t):
    return Lorenz(X)

Y0 = np.array([1.0, 1.0, 1.0]) #x, y, z
t = np.arange(0, t_Fin, t_step)

soln = odeint(f, Y0, t)
x = soln[:, 0]; y = soln[:, 1]; z = soln[:, 2]

# Calculate the maximum Lyapunov exponent.
# Source: Clyde-Emmanuel Estorninho Meador. _Numerical Calculation
# of 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], ord=1) +\
        np.linalg.norm(L_f[1], ord=1) + np.linalg.norm(L_f[2], ord=1)
    print L_S/((N - i_L)*3*t_step) #p. 33

============================

After the "if Lyap" statement, each line does the following:
(1; "x_0" line) Toss the early data.
(2; "x_Lap" line) Shift each point in the numerical solution by an
arbitrarily small amount. Separated in case I want to shift x, y, and z by
different amounts.
(3; "delta" line) Integrates one time step via Euler's method.
(4; "L_f" line) Formula from the paper (p. 32).
(5; "L_f_temp" line) (Dummy line--debugging.)
(6; "L_S" line) Formula from the paper (p. 33).
(7; "print" line) Ditto.

I really don't see what's wrong here. I have tried so many times to debug
this, and I still after all these attempts, cannot get a correct Lyapunov
exponent. :(