我想通过使用此Lorenz System中描述的标准方法来数值计算Paper, p.81的Lyapunov谱。
基本上需要集成Lorenz系统和切向矢量(为此我使用了Runge-Kutta方法)。切向向量的演化方程由Lorenz系统的Jacobi矩阵给出。每次迭代后,需要对向量应用Gram-Schmidt方案并存储其长度。然后,通过存储长度的平均值给出三个Lyapunov指数。
我在python(使用版本3.7.4)中实现了上述方案,但没有得到正确的结果。
我的错误在于der向量的Rk4-方法,但我找不到任何错误...轨迹x,y,z的RK4-方法正确运行(由图表示),并且已实现Gram-施密特方案也得到了正确实施。
我希望有人可以浏览我的短代码并找到我的错误
from numpy import array, arange, zeros, dot, log
import matplotlib.pyplot as plt
from numpy.linalg import norm
# Lorenz-System: x' = sigma*(y-x)
# y' = x*(rho - z) - y
# z' = x*y - beta*z
# Evolution equation of tracjectories and tangential vectors
def f(r):
x = r[0]
y = r[1]
z = r[2]
fx = sigma * (y - x)
fy = x * (rho - z) - y
fz = x * y - beta * z
return array([fx,fy,fz], float)
def g(d, r):
dx = d[0]
dy = d[1]
dz = d[2]
dfx = dot(M, dx)
dfy = dot(M, dy)
dfz = dot(M, dz)
return array([dfx, dfy, dfz], float)
# Jacobi Matrix of Lorenz System
def jacobian(r):
M = zeros([3,3])
M[0,:] = [- sigma, sigma, 0]
M[1,:] = [rho - r[2], -1, - r[0] ]
M[2,:] = [r[1], r[0], -beta]
return M
# Initial conditions
d = array([[1,0,0], [0,1,0], [0,0,1]], float)
r = array([10.0, 1.0, 0.0], float)
sigma, rho, beta = 10, 45.92, 4.0
T = 10**5 # time steps and increment
dt = 0.01 # time increment
Teq = 10**4 # Transient time
l1, l2, l3 = 0, 0, 0 # Lengths
xpoints, ypoints, zpoints = [], [], []
# Transient
for t in range(Teq):
k1 = dt * f(r) #RK4 for trajectories x,y,z
k2 = dt * f(r + 0.5 * k1)
k3 = dt * f(r + 0.5 * k2)
k4 = dt * f(r + k3)
r += (k1 + 2 * k2 + 2 * k3 + k4) / 6
M = jacobian(r)
k11 = dt * g(d, r) #RK4 for tangent vectors
k22 = dt * g(d + 0.5 * k1, r)
k33 = dt * g(d + 0.5 * k2, r)
k44 = dt * g(d + k3, r)
d += (k11 + 2 * k22 + 2 * k33 + k44) / 6
xpoints.append(r[0]), ypoints.append(r[1]), zpoints.append(r[2])
orth_1 = d[0] # Gram-Schmidt-Scheme
d[0] = orth_1 / norm(orth_1)
orth_2 = d[1] - dot(d[1], d[0]) * d[0]
d[1] = orth_2 / norm(orth_2)
orth_3 = d[2] - (dot(d[2], d[1]) * d[1]) - (dot(d[2], d[0]) * d[0])
d[2] = orth_3 / norm(orth_3)
for t in range(T):
k1 = dt * f(r)
k2 = dt * f(r + 0.5 * k1)
k3 = dt * f(r + 0.5 * k2)
k4 = dt * f(r + k3)
r += (k1 + 2 * k2 + 2 * k3 + k4) / 6
k11 = dt * g(d, r)
k22 = dt * g(d + 0.5 * k1, r)
k33 = dt * g(d + 0.5 * k2, r)
k44 = dt * g(d + k3, r)
d += (k11 + 2 * k22 + 2 * k33 + k44) / 6
orth_1 = d[0]
d[0] = orth_1 / norm(orth_1)
l1 += log(norm(orth_1)) # add lengths up
orth_2 = d[1] - dot(d[1], d[0]) * d[0]
d[1] = orth_2 / norm(orth_2)
l2 += log(norm(orth_2))
orth_3 = d[2] - (dot(d[2], d[1]) * d[1]) - (dot(d[2], d[0]) * d[0])
d[2] = orth_3 / norm(orth_3)
l3 += log(norm(orth_3))
plt.plot(xpoints, ypoints)
plt.show()
# Correct Solution (2.16, 0.0, -32.4)
lya1 = l1 / (dt * T)
lya2 = l2 / (dt * T)
lya3 = l3 / (dt * T)
lya1, lya2, lya3
#My result (-3.286149130588001e-15, 12.182263788663551, -9.503287701985853)
您需要将点和雅可比矩阵解为(正)耦合系统。在原始源代码中,所做的一切恰恰是在组合系统的一次RK4
调用中更新的。
例如,在第二阶段,您将混合操作以合并第二阶段
k2 = dt * f(r + 0.5 * k1)
M = jacobian(r + 0.5 * k1)
k22 = dt * g(d + 0.5 * k11, r + 0.5 * k1)
您还可以在M
函数中委派g
的计算,因为这是唯一需要它的地方,并且可以增加变量范围内的局部性。
注意,我将d
的更新从k1
更改为k11
,这应该是数值结果中错误的主要来源。