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

[Lou Pecora]

Do your calculations agree with published Lyapunov exponents?

 
-- Lou Pecora


On Monday, July 21, 2014 2:56 PM, Barrett B <barrett.n.b at gmail.com> wrote:
 


Alright, guys, I FINALLY got it to work. Thanks to all for the help. The
current code is (ignore the annoying line returns associated with comments):


import numpy as np
#import matplotlib.pyplot as plt
#import matplotlib.mlab as mlab
from scipy.integrate import *

#System constants
rates = [-0.1, -0.4, -1]
sigma = 10; beta = 8.0/3; rho = 28

#Simulation constants
timeStep = 1e-3;
tStartTrack = 20; #Time when Lyapunov exponent calculation begins.
nStartTrack = int(tStartTrack/timeStep) #init step for Lyap expo
tFinal = 4000; #Final timestep.
nFinal = int(tFinal/timeStep) #a.k.a N
Lyap = True #Whether to calculate the largest Lyapunov exponent.

#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. Lorenz is time independent, but odeint requires the time to be input.
def f(X, t):
    return Lorenz(X)
    
#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

#Script starts here.
Y0 = np.array([1.0, 1.0, 1.0]) #x, y, z
if Lyap:
    t = np.arange(0, tStartTrack, timeStep) #Will use RK4 the rest of the way.
else:
    t = np.arange(0, tFinal, timeStep)

soln = odeint(f, Y0, t)

if Lyap:
    L_s = 0.0; saves = 0; timesUsed = 0
    offsetArray = np.array([1e-6, 1e-6, 1e-6])
    tInterval = 10.0 #Time length of each interval. Reset soln_1 after each one.
    numberOfIntervals = int((tFinal - tStartTrack)/tInterval)

    pointsPerInterval = int(tInterval/timeStep) #Data points per interval.
    soln_0 = soln[nStartTrack-1, :] #The original solution
    soln_1 = soln_0 + offsetArray #The perturbed solution
    #Second norm ("None") of the distance, before some integration:
    d_0 = np.linalg.norm(offsetArray, ord=None)

    print 'Beginning Lyapunov calculation'
    
    for currentInterval in range (numberOfIntervals):
        if currentInterval>0:
            #Rescale soln_1 such that its dist from soln_0 = that of eps_L's.
            newPerturb = (soln_1 - soln_0) *\
                np.linalg.norm(offsetArray, ord=None)/d_1
            soln_1 = soln_0 + newPerturb
        
        for j in range(pointsPerInterval):
            #Update status.
            statusPct = 100.0*(currentInterval + 1)/numberOfIntervals

            #Use RK4 to integrate this ONE timestep.
            soln_1 = RK4(soln_1, timeStep)
            soln_0 = RK4(soln_0, timeStep)
            
        #Second norm ("None") of the distance, after some integration:
        d_1 = np.linalg.norm(soln_0 - soln_1, ord=None)

        #Keep a running total of the natural log of the ratio d_1/d_0.
        if d_1/d_0 > 1e-20: #just in case.
            L_s += np.log(d_1/d_0)
            timesUsed += 1
        else:
            saves += 1
    
    #Divide by the number of data points and the timestep.
    #This is the largest Lyapunov exponent.
    print L_s/tInterval/numberOfIntervals

_______________________________________________
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/e9d18491/attachment.html>