如何解决3D系统的Lyapunov指数谱
我一直在寻找用于计算Lyapunov指数的Python代码，最后找到了代码LyapunovExponets，但是它很长而且没有向量化，并且没有使用python 3和ODE求解器。我需要帮助将其向量化，并使用solv_ivp代替使用的那个。
我打算将此代码用于也是3D ODE系统的模型，但几乎不可能理解此代码或对其进行修改以适用于一般3D ODE系统。我还希望能够使用其他ODE因为RK方法并不能解决所有问题。代码也很长，所以我希望它更紧凑，当然可以加快速度。这对于研究动力系统的其他人也将是一件好事。当然，有一些Matlab代码可以计算lyapunov指数，但我想使用开源语言编写代码。
# LorenzODELCE.py:
# Estimate the spectrum of Lyapunov Characteristic Exponents
# for the Lorenz ODEs,using the pull-back method.
# Also,estimate the volume-contraction (dissipation) rate and the
# fractal dimenion (latter using the Kaplan-Yorke conjecture).
# Plot out trajectory,for reference.
#
# Comment:
#Notice how much more complicated the code has become,given
#that we're writing out variables in component form.
# This should be rewritten to use vectors,which will be
# much more compact and easier to debug. Equally important,# the code would generalize to any dimension system.
# Import plotting routines
from pylab import *
# The Lorenz 3D ODEs
#Original parameter values: (sigma,R,b) = (10,28,-8/3)
def LorenzXDot(sigma,b,x,y,z):
return sigma * (-x + y)
def LorenzYDot(sigma,z):
return R*x - x*z - y
def LorenzZDot(sigma,z):
return -b*z + x*y
# The tangent space (linearized) flow (aka co-tangent flow)
def LorenzDXDot(sigma,z,dx,dy,dz):
return sigma * (-dx + dy)
def LorenzDYDot(sigma,dz):
return (R-z)*dx - dy - x*dz
def LorenzDZDot(sigma,dz):
return y*dx + x*dy + b*dz
# Volume contraction given by
# Trace(Jacobian(x,z)) = b - sigma - 1
def LorenzODETrJac(sigma,z):
return b - sigma - 1
# As a check,we must have total contraction = Sum of LCEs
# Tr(J) = Sum_i LCEi
# Numerical check: at (sigma,-8/3)
# LCE0 ~ 0.9058
# LCE1 ~ 0.0000
# LCE2 ~ -14.572
# Tr(J) ~ -13.6666
# These use base-2 logs
# The fractal dimension from the LCEs (Kaplan-Yorke conjecture)
# Assume these are ordered: LCE1 >= LCE2 >= LCE3
def FractalDimension3DODE(LCE1,LCE2,LCE3):
# "Close" to zero ... we're estimating here
Err = 0.01
if LCE1 < -Err: # Stable fixed point (-,-,-)
return 0.0
elif abs(LCE1) <= Err:
if LCE2 < -Err: # Limit cycle (0,-)
return 1.0
else: # Torus (0,-)
return 2.0
else: # Chaotic attractor (+,-)
return 2.0 + (LCE1+LCE2) / abs(LCE3)
# 3D fourth-order Runge-Kutta integrator
def RKThreeD(a,c,f,g,h,dt):
k1x = dt * f(a,z)
k1y = dt * g(a,z)
k1z = dt * h(a,z)
k2x = dt * f(a,x + k1x / 2.0,y + k1y / 2.0,z + k1z / 2.0)
k2y = dt * g(a,z + k1z / 2.0)
k2z = dt * h(a,z + k1z / 2.0)
k3x = dt * f(a,x + k2x / 2.0,y + k2y / 2.0,z + k2z / 2.0)
k3y = dt * g(a,z + k2z / 2.0)
k3z = dt * h(a,z + k2z / 2.0)
k4x = dt * f(a,x + k3x,y + k3y,z + k3z)
k4y = dt * g(a,z + k3z)
k4z = dt * h(a,z + k3z)
x += ( k1x + 2.0 * k2x + 2.0 * k3x + k4x ) / 6.0
y += ( k1y + 2.0 * k2y + 2.0 * k3y + k4y ) / 6.0
z += ( k1z + 2.0 * k2z + 2.0 * k3z + k4z ) / 6.0
return x,z
# Tanget space flow (using fourth-order Runge-Kutta integrator)
def TangentFlowRKThreeD(a,df,dg,dh,dz,dt):
k1x = dt * df(a,dz)
k1y = dt * dg(a,dz)
k1z = dt * dh(a,dz)
k2x = dt * df(a,dx+k1x/2.0,dy+k1y/2.0,dz+k1z/2.0)
k2y = dt * dg(a,dz+k1z/2.0)
k2z = dt * dh(a,dz+k1z/2.0)
k3x = dt * df(a,dx+k2x/2.0,dy+k2y/2.0,dz+k2z/2.0)
k3y = dt * dg(a,dz+k2z/2.0)
k3z = dt * dh(a,dz+k2z/2.0)
k4x = dt * df(a,dx+k3x,dy+k3y,dz+k3z)
k4y = dt * dg(a,dz+k3z)
k4z = dt * dh(a,dz+k3z)
dx += ( k1x + 2.0 * k2x + 2.0 * k3x + k4x ) / 6.0
dy += ( k1y + 2.0 * k2y + 2.0 * k3y + k4y ) / 6.0
dz += ( k1z + 2.0 * k2z + 2.0 * k3z + k4z ) / 6.0
return dx,dz
# Simulation parameters
# Integration time step
dt = 0.01
#
# Control parameters for the Lorenz ODEs:
sigma = 10.0
R = 28.0
b = -8.0/3.0
# The number of iterations to throw away
nTransients = 10
# The number of time steps to integrate over
nIterates = 1000
# The main loop that generates the orbit,storing the states
xState = 5.0
yState = 5.0
zState = 5.0
# Iterate for some number of transients,but don't use these states
for n in range(0,nTransients):
xState,yState,zState = RKThreeD(sigma,xState,zState,LorenzXDot,LorenzYDot,LorenzZDot,dt)
# Set up array of iterates and store the current state
x = [xState]
y = [yState]
z = [zState]
for n in range(0,nIterates):
# at each time step calculate new (x,z)(t)
xt,yt,zt = RKThreeD(sigma,x[n],y[n],z[n],dt)
# and append to lists
x.append(xt)
y.append(yt)
z.append(zt)
# Estimate the LCEs
# The number of iterations to throw away
nTransients = 100
# The number of iterations to over which to estimate
# This is really the number of pull-backs
nIterates = 10000
# The number of iterations per pull-back
nItsPerPB = 10
# Initial condition
xState = 5.0
yState = 5.0
zState = 5.0
# Initial tangent vectors
e1x = 1.0
e1y = 0.0
e1z = 0.0
e2x = 0.0
e2y = 1.0
e2z = 0.0
e3x = 0.0
e3y = 0.0
e3z = 1.0
# Iterate away transients and let the tangent vectors align
# with the global stable and unstable manifolds
for n in range(0,nTransients):
for i in range(nItsPerPB):
xState,\
LorenzXDot,dt)
# Evolve tangent vector for maximum LCE (LCE1)
e1x,e1y,e1z = TangentFlowRKThreeD(sigma,\
LorenzDXDot,LorenzDYDot,LorenzDZDot,e1x,e1z,dt)
# Evolve tangent vector for next LCE (LCE2)
e2x,e2y,e2z = TangentFlowRKThreeD(sigma,e2x,e2z,dt)
# Evolve tangent vector for last LCE
e3x,e3y,e3z = TangentFlowRKThreeD(sigma,e3x,e3z,dt)
# Normalize the tangent vector
d = sqrt(e1x*e1x + e1y*e1y + e1z*e1z)
e1x /= d
e1y /= d
e1z /= d
# Pull-back: Remove any e1 component from e2
dote1e2 = e1x * e2x + e1y * e2y + e1z * e2z
e2x -= dote1e2 * e1x
e2y -= dote1e2 * e1y
e2z -= dote1e2 * e1z
# Normalize second tangent vector
d = sqrt(e2x*e2x + e2y*e2y + e2z*e2z)
e2x /= d
e2y /= d
e2z /= d
# Pull-back: Remove any e1 and e2 components from e3
dote1e3 = e1x * e3x + e1y * e3y + e1z * e3z
dote2e3 = e2x * e3x + e2y * e3y + e2z * e3z
e3x -= dote1e3 * e1x + dote2e3 * e2x
e3y -= dote1e3 * e1y + dote2e3 * e2y
e3z -= dote1e3 * e1z + dote2e3 * e2z
# Normalize third tangent vector
d = sqrt(e3x*e3x + e3y*e3y + e3z*e3z)
e3x /= d
e3y /= d
e3z /= d
# Okay,now we're ready to begin the estimation
LCE1 = 0.0
LCE2 = 0.0
LCE3 = 0.0
for n in range(0,nIterates):
for i in range(nItsPerPB):
xState,dt)
# Normalize the tangent vector
d = sqrt(e1x*e1x + e1y*e1y + e1z*e1z)
e1x /= d
e1y /= d
e1z /= d
# Accumulate the first tangent vector's length change factor
LCE1 += log(d)
# Pull-back: Remove any e1 component from e2
dote1e2 = e1x * e2x + e1y * e2y + e1z * e2z
e2x -= dote1e2 * e1x
e2y -= dote1e2 * e1y
e2z -= dote1e2 * e1z
# Normalize second tangent vector
d = sqrt(e2x*e2x + e2y*e2y + e2z*e2z)
e2x /= d
e2y /= d
e2z /= d
# Accumulate the second tangent vector's length change factor
LCE2 += log(d)
# Pull-back: Remove any e1 and e2 components from e3
dote1e3 = e1x * e3x + e1y * e3y + e1z * e3z
dote2e3 = e2x * e3x + e2y * e3y + e2z * e3z
e3x -= dote1e3 * e1x + dote2e3 * e2x
e3y -= dote1e3 * e1y + dote2e3 * e2y
e3z -= dote1e3 * e1z + dote2e3 * e2z
# Normalize third tangent vector
d = sqrt(e3x*e3x + e3y*e3y + e3z*e3z)
e3x /= d
e3y /= d
e3z /= d
# Accumulate the third tangent vector's length change factor
LCE3 += log(d)
# Convert to per-iterate,per-second LCEs and to base-2 logs
IntegrationTime = dt * float(nItsPerPB) * float(nIterates)
LCE1 = LCE1 / IntegrationTime
LCE2 = LCE2 / IntegrationTime
LCE3 = LCE3 / IntegrationTime
# Calculate contraction factor,for comparison.
# For Lorenz ODE,we know this is independent of (x,z).
# Otherwise,we'd have to estimate it along the trajectory,too.
Contraction = LorenzODETrJac(sigma,0.0,0.0)
# Choose a pair of coordinates from (x,z) to show
# Setup the parametric plot:
xlabel('x(t)') # set x-axis label
ylabel('y(t)') # set y-axis label
# Construct plot title
LCEString = '(%g,%g,%g)' % (LCE1,LCE3)
PString = '($\sigma$,b) = (%g,%g)' % (sigma,b)
CString = 'Contraction = %g,Diff = %g' % (Contraction,abs(LCE1+LCE2+LCE3-Contraction))
FString = 'Fractal dimension = %g' % FractalDimension3DODE(LCE1,LCE3)
title('4th order Runge-Kutta Method: Lorenz ODE at ' + PString + ':\n LCEs = ' + LCEString + ',' + CString + '\n ' + FString)
axis('equal')
axis([-20.0,20.0,-20.0,20.0])
# Plot the trajectory in the phase plane
plot(x,'b')
axhline(0.0,color = 'k')
axvline(0.0,color = 'k')
# Use command below to save figure
#savefig('LorenzODELCE',dpi=600)
# Display the plot in a window
show()
我尝试并花了很多时间去做，但是在某个时候迷路了。这是我尝试过的
%matplotlib inline
#%matplotlib notebook
#%pylab
from matplotlib import use
#use("Qt5Agg")
import numpy as np
from numpy import *
from cmath import*
from numpy import linalg as LA
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.signal import find_peaks
from scipy import signal
import scipy as scp
from scipy.integrate import solve_ivp,odeint
import pandas as pd
from scipy import linalg as LA
import nolds
# The Lorenz 3D ODEs
#Original parameter values: (sigma,-8/3)
def Lorenz_Jac(t,X,sigma,b):
LzJac=np.array([[-sigma,0],[R-z,-1,-x],[y,-b]])
return LzJac
def LorenzFun(t,b):
x,z=X
dxdt= sigma * (-x + y)
dydt=R*x - x*z - y
dzdt=b*z + x*y
return[dxdt,dydt,dzdt]
# The tangent space (linearized) flow (aka co-tangent flow)
def LorenzDXDot(t,dxyz,z=X
dx,dz=dxyz
return[sigma * (-dx + dy),(R-z)*dx - dy - x*dz,y*dx + x*dy + b*dz]
# Volume contraction given by
# Trace(Jacobian(x,X):
x,z=X
tr=b - sigma - 1
return tr
# As a check,we must have total contraction = Sum of LCEs
# Tr(J) = Sum_i LCEi
# Numerical check: at (sigma,-8/3)
# LCE0 ~ 0.9058
# LCE1 ~ 0.0000
# LCE2 ~ -14.572
# Tr(J) ~ -13.6666
# These use base-2 logs
# The fractal dimension from the LCEs (Kaplan-Yorke conjecture)
# Assume these are ordered: LCE1 >= LCE2 >= LCE3
def FractalDimension3DODE(LCE1,-)
return 1.0
else: # Torus (0,-)
return 2.0
else: # Chaotic attractor (+,-)
return 2.0 + (LCE1+LCE2) / abs(LCE3)
# 3D fourth-order Runge-Kutta integrator
# Tanget space flow (using fourth-order Runge-Kutta integrator)
# Simulation parameters
# Integration time step
dt = 0.01
#
# Control parameters for the Lorenz ODEs:
sigma = 10.0
R = 28.0
b = -8.0/3.0
# The number of iterations to throw away
nTransients = 10
# The number of time steps to integrate over
nIterates = 1000
# The main loop that generates the orbit,storing the states
X0=np.array([ 5.0,5,5])
# Iterate for some number of transients,but don't use these states
t0=0
tMax=1
t=[t0,tMax]
t_eval_pts=np.linspace(t0,tMax,10*tMax)
sol=solve_ivp(Lorenz,t,X0,method= 'RK45',t_eval= t_eval_pts,args=(sigma,b),jac=Lorenz_Jac,rtol=1e-08,atol=1e-08);
x,z=sol.y
# Set up array of iterates and store the current state
# at each time step calculate new (x,z)(t)
dt=0.1
T=100
t1=[0,T]
solt = solve_ivp(Lorenz,t1,[x[-1],y[-1],z[-1] ],t_eval= np.arange(0,T,dt),atol=1e-08) ;
xt,zt=solt.y
# and append to lists
# Estimate the LCEs
# The number of iterations to throw away
nTransients = 100
# The number of iterations to over which to estimate
# This is really the number of pull-backs
nIterates = 1000
# The number of iterations per pull-back
nItsPerPB = 10
# Initial condition
xState = 5.0
yState = 5.0
zState = 5.0
# Initial tangent vectors
e1 = [1.0,0]
e2 = [0,1,0]
e3=[0,1]
# Iterate away transients and let the tangent vectors align
#with the global stable and unstable manifolds
# Evolve tangent vector for maximum LCE (LCE1)
# Okay,nIterates):
xState,zState =
# Evolve tangent vector for maximum LCE (LCE1)
e1x,e1z =
# Evolve tangent vector for next LCE (LCE2)
e2x,e2z =
# Evolve tangent vector for last LCE
e3x,e3z =
sol1= solve_ivp(Lorenz,LorenzDXDot,atol=1e-08)
e1x,e1z=sol1.y
sol2= solve_ivp(Lorenz,atol=1e-08)
e2x,e2z=sol2.y
sol3= solve_ivp(Lorenz,atol=1e-08)
e3x,e3z=sol3.y
# Normalize the tangent vector
d=LA.norm(e1)
e1x /= d
e1y /= d
e1z /= d
# Accumulate the first tangent vector's length change factor
LCE1 += log(d)
# Pull-back: Remove any e1 component from e2
e2x -= np.dot(e1,e2) * e1x
e2y -= np.dot(e1,e2) * e1y
e2z -= np.dot(e1,e2) * e1z
# Normalize second tangent vector
d=LA.norm(e2)
e2x /= d
e2y /= d
e2z /= d
# Accumulate the second tangent vector's length change factor
LCE2 += log(d)
# Pull-back: Remove any e1 and e2 components from e3
e3x -= np.dot(e1,e3) * e1x + np.dot(e2,e3) * e2x
e3y -= np.dot(e1,e3) * e1y + np.dot(e2,e3) * e2y
e3z -= np.dot(e1,e3) * e1z + np.dot(e2,e3)* e2z
# Normalize third tangent vector
d=LA.norm(e2)
e3x /= d
e3y /= d
e3z /= d
# Accumulate the third tangent vector's length change factor
LCE3 += log(d)
# Convert to per-iterate,per-second LCEs and to base-2 logs
#IntegrationTime = dt * float(nItsPerPB) * float(nIterates)
#LCE1 = LCE1 / IntegrationTime
#LCE2 = LCE2 / IntegrationTime
#LCE3 = LCE3 / IntegrationTime
# Calculate contraction factor,for comparison.
#For Lorenz ODE,z).
#Otherwise,dpi=600)
# Display the plot in a window
show()
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。