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

[Barrett B]

Lou Pecora <lou_boog2000 <at> yahoo.com> writes:

> 
> 
> Hi, Barrett,
> 
> I added some needed tabs and a calculation of lendx (length of dx) which I
left out. You can see by that why I want you to understand the math.  Even
experienced people like me can make code mistakes.   
> 
> def LypExp_calc(vics, n, T, NT):    lnsum= 0.0      x= vics      for i in
xrange(NT):        temp= NP.RN.random_sample((n,))        templen=
sqrt(dot(temp,temp))        dx= temp/templen       
y=NP.concatenate([x,dx])    	yT= YourODEintegrator(y,T, vecf)    	dxT=
yT[n:2*n] 
> 
> 	lendx= sqrt(NP.dot(dxT,dxT))  #  Line added    	lnsum += log(lendx)    
x= yT[:n]    LypExp= lnsum/(NT*T)    return LypExpdef vecf(y):    n=3    x=
y[:n]    dx= y[n:2*n]    vf= F(x)    Dvf= NP.dot(DF(x),dx)    return
NP.concatenate([vf, Dvf])
> 
>  
> -- Lou Pecora 
> 

Thanks, Lou, but I'm really not sure how to get this up and running. If it's
not obvious by now, my programming skills are very, very weak. :( For
example, this is what I have been using for my integrator:

#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
    
#Time-independent RK4 method.
def RK4(X, timeStep):
    k1 = Lorenz(X)
    k2 = Lorenz(X + timeStep*k1 / 2)
    k3 = Lorenz(X + timeStep*k2 / 2)
    k4 = Lorenz(X + timeStep*k3)       
    return X + timeStep*(k1 + 2*k2 + 2*k3 + k4)/6

And this won't work, because there are six variables to integrate, not
three. Help!