Computing Lyapunov Spectrum of the Lorenz Attraktor with Gram Schmidt and variational equation

General Usage
1

Hello,
I’m new to this forum. I’m studying physics at Cologne university and for my thesis I want to calculate Lyapunovexponents of several systems. I choosed the Lorentz Attraktor at first to program the way to calculate the Lyapunovexponent and to check if my result is ok.
I tried to use a method by using the variational equation by integrating this equation and using Gram-Schmidt orthonormalization.
My code is below. I didn’t explain everything but maybe some can find my fault.
Thank you very much!

    x=X[1]
    y=X[2]
    z=X[3]
    x_=10*(y-x) #16#10
    y_=x*(28-z)-y#45.92,28
    
    z_=x*y-8/3*z #4#8/3
    return([x_,y_,z_])
end
function rungestep_1(X,h)
    #t[i+1]=t[i]+h
        K1=lorentz(X)
       K2=lorentz(X+h*K1/2)
        K3=lorentz(X+h*K2/2)
        K4=lorentz(X+h*K3)
       
    return(X+h/6*(K1+2*K2+2*K3+K4))
 
end
function RungeKutta_2(x0,v_0,N,M,h)
    
    t=zeros(N)
   
     P=zeros(9,N)
    P[:,1]=v_0
    T=zeros(9,N)
  z=M*N
    
     B=zeros(3,z)
      R=zeros(9,z)
    G=zeros(9,N)
    R[:,1]=v_0
     G[:,1]=v_0
    B[:,1]=x0
    s=0
    L=zeros(N)
 for u in 1:N-3  
       
          β=8/3
        σ=10
        ρ=28
         B[:,u+1]=rungestep_1(B[:,u],h)
        for i in 1:M
           
          #  B[:,s+i+1]=rungestep_1(B[:,i+s],h)

        
        
        x=B[1,u]
        y=B[2,u]
        z=B[3,u]
            if i==1
           R[:,i+s]=rungestep_2(G[:,u],h,x,y,z,σ,ρ,β)
            elseif i>1
        R[:,i+s]=rungestep_2(R[:,i+s-1],h,x,y,z,σ,ρ,β)
            end
            
        end
    r=M*u
   g=R[:,r]
       
    l=gramschmidt2(matrix(g))[2]
       P[:,u+1]=vector(gramschmidt2(matrix(g))[1])
        
       
   G[:,u+1]=vector(l)
       L[u+1]=log(norm( [P[2,u],P[5,u],P[8,u] ]))
        s=s+M
       
        
       
end
   return sum(L)/(N*π/100)
end

function gramschmidt2(M)
    v1=M[:,1]
    v2=M[:,2]
    v3=M[:,3]
    w1=v1
    e1=w1/norm(w1)
    w2=v2-proj(e1,v2)
    e2=w2/norm(w2)
    w3=v3-proj(e1,v3)-proj(e2,v3)
    e3=w3/norm(w3)
u=[ [w1[1],w2[1],w3[1]] [w1[2],w2[2],w3[2]] [w1[3],w2[3],w3[3]] ]
    e=[ [e1[1],e2[1],e3[1]] [e1[2],e2[2],e3[2]] [e1[3],e2[3],e3[3]] ]
    return u,e
end
2 Likes
2

I can’t comment on the details of your code, but rather than write your own Runge–Kutta integrator you will be far better off using DifferentialEquations.jl, rather than implementing your own Gram–Schmidt orthogonalization you are much better off using QR factorization via Julia’s qr function, and rather than implementing your own algorithm to compute Lyapunov exponents you will likely be better off using the lyapunovspectrum function in DynamicalSystems.jl.

Even if your educational program requires you to re-implement these algorithms yourself, it would still be a good idea to debug your code by validating it against mature libraries.

3 Likes