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

[Lou Pecora]

Dear Barrett,

I suspect your problem is not so much weak programming skills, but a weak understanding of dynamical systems and ODEs.  Sorry, I do not mean to sound negative. I think you really need to study more about these kinds of systems and how to solve them.  I would guess the RK solver is something you wrote.  It is not how you would do Runge-Kutta.  You could get by quite well with odeint in the scipy library. It's a good routine for solving ODEs.  But you are hampered by not knowing what the Jacobian is.  It provides the other three equations you need from the variational part of the problem.  I tried to explain it, but I guess I didn't get it through.  Please take a step back and put some time in to understand these concepts and methods. Going at it by just hoping someone's code will work is not the answer, but it can cause you lots of problems as you are seeing.  I don't know why you are trying to calculate the Lyapunov exponents, but if this is
 something that is important in your job or school, then a deeper understanding will be the most beneficial. 

I wish you well in your endeavor.

Sincerely yours,

-- Lou Pecora


On Sunday, July 20, 2014 7:22 PM, Barrett B <barrett.n.b at gmail.com> wrote:
 


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!

_______________________________________________
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/20140721/09fbc9cf/attachment.html>