问题描述
我想使用标准方法数值计算洛伦兹系统的李雅普诺夫谱这在论文，第81页中有所描述.
I want to numerically compute the Lyapunov Spectrum of the Lorenz System by using the standard method which is described in this Paper, p.81.
一个基本上需要整合洛伦兹系统和切向向量(我为此使用了 Runge-Kutta 方法).切向向量的演化方程由洛伦兹系统的雅可比矩阵给出.每次迭代后，需要对向量应用 Gram-Schmidt 方案并存储其长度.三个李雅普诺夫指数然后由存储长度的平均值给出.
One basically need integrate the Lorenz system and the tangential vectors (i used the Runge-Kutta method for this). The evolution equation of the tangential vectors are given by the Jacobi matrix of the Lorenz system. After each iterations one needs to apply the Gram-Schmidt scheme on the vectors and store its lengths. The three Lyapunov exponents are then given by the averages of the stored lengths.
我在python(使用3.7.4版)中实现了上面解释的方案，但我没有得到正确的结果.
I implemented the above explained scheme in python (used version 3.7.4), but I don't get the correct results.
我认为错误在于 der 向量的 Rk4 方法，但我找不到任何错误......轨迹 x、y、z 的 RK4 方法正常工作(由图指示)和实现的 Gram-施密特方案也正确实施.
I thing the bug lies in the Rk4-Method for der vectors, but i cannot find any error...The RK4-method for the trajectories x,y,z works correctly (indicated by the plot) and the implemented Gram-Schmidt scheme is also correctly implemented.
我希望有人能看一下我的短代码，也许会发现我的错误
I hope that someone could look through my short code and maybe find my error
更新代码
from numpy import array, arange, zeros, dot, log
import matplotlib.pyplot as plt
from numpy.linalg import norm
# 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 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
def g(d, r):
dx = d[0]
dy = d[1]
dz = d[2]
M = jacobian(r)
dfx = dot(M, dx)
dfy = dot(M, dy)
dfz = dot(M, dz)
return array([dfx, dfy, dfz], float)
# Initial conditions
d = array([[1,0,0], [0,1,0], [0,0,1]], float)
r = array([19.0, 20.0, 50.0], float)
sigma, rho, beta = 10, 45.92, 4.0
T = 10**5 # time steps
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):
# RK4 - Method
k1 = dt * f(r)
k11 = dt * g(d, r)
k2 = dt * f(r + 0.5 * k1)
k22 = dt * g(d + 0.5 * k11, r + 0.5 * k1)
k3 = dt * f(r + 0.5 * k2)
k33 = dt * g(d + 0.5 * k22, r + 0.5 * k2)
k4 = dt * f(r + k3)
k44 = dt * g(d + k33, r + k3)
r += (k1 + 2 * k2 + 2 * k3 + k4) / 6
d += (k11 + 2 * k22 + 2 * k33 + k44) / 6
# Gram-Schmidt-Scheme
orth_1 = d[0]
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)
k11 = dt * g(d, r)
k2 = dt * f(r + 0.5 * k1)
k22 = dt * g(d + 0.5 * k11, r + 0.5 * k1)
k3 = dt * f(r + 0.5 * k2)
k33 = dt * g(d + 0.5 * k22, r + 0.5 * k2)
k4 = dt * f(r + k3)
k44 = dt * g(d + k33, r + k3)
r += (k1 + 2 * k2 + 2 * k3 + k4) / 6
d += (k11 + 2 * k22 + 2 * k33 + k44) / 6
orth_1 = d[0] # Gram-Schmidt-Scheme
l1 += log(norm(orth_1))
d[0] = orth_1 / norm(orth_1)
orth_2 = d[1] - dot(d[1], d[0]) * d[0]
l2 += log(norm(orth_2))
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])
l3 += log(norm(orth_3))
d[2] = orth_3 / norm(orth_3)
# Correct Solution (2.16, 0.0, -32.4)
lya1 = l1 / (dt * T)
lya2 = l2 / (dt * T) - lya1
lya3 = l3 / (dt * T) - lya1 - lya2
lya1, lya2, lya3
# my solution T = 10^5 : (1.3540301507934012, -0.0021967491623752448, -16.351653561383387)
以上代码根据Lutz的建议更新.结果看起来好多了，但仍然不是 100% 准确.
The above code is updated according to Lutz suggestions.The results look much better but they are still not 100% accurate.
正确的解决方案 (2.16, 0.0, -32.4)
Correct Solution (2.16, 0.0, -32.4)
我的解决方案 (1.3540301507934012, -0.0021967491623752448, -16.351653561383387)
My solution (1.3540301507934012, -0.0021967491623752448, -16.351653561383387)
正确的解决方案来自 Wolf's Paper, p.289.在第 290-291 页，他描述了他的方法，看起来与我在本文开头提到的论文(论文，第 81 页)中的完全相同.
The correct solutions are from Wolf's Paper, p.289. On page 290-291 he describes his method, which looks exactly the same as in the paper that i mentioned in the beginning of this post (Paper, p.81).
所以我的代码中肯定有另一个错误......
So there must be another error in my code...
推荐答案
您需要将点系统和雅可比系统作为(前向)耦合系统求解.在完全完成的原始源代码中，所有内容都在组合系统的一个 RK4 调用中更新.
You need to solve the system of point and Jacobian as the (forward) coupled system that it is. In the original source exactly that is done, everything is updated in one RK4 call for the combined system.
例如在第二阶段，您可以混合操作以获得组合的第二阶段
So for instance in the second stage, you would mix the operations to have a combined second stage
k2 = dt * f(r + 0.5 * k1)
M = jacobian(r + 0.5 * k1)
k22 = dt * g(d + 0.5 * k11, r + 0.5 * k1)
您也可以在 g 函数内部委托 M 的计算，因为这是唯一需要它的地方，并且您增加了变量范围内的局部性.
You could also delegate the computation of M inside the g function, as this is the only place where it is needed, and you increase locality in the scope of variables.
注意我把d的update从k1改成了k11，这应该是数值结果错误的主要来源.
Note that I changed the update of d from k1 to k11, which should be the main source of the error in the numerical result.
对最新代码版本 (2/28/2021) 的补充说明:
正如评论中所说，代码看起来会按照算法的数学规则执行.有两个误读会阻止代码返回接近引用的结果:
As said in the comments, the code looks like it will do what the mathematics of the algorithm prescribes. There are two misreadings that prevent the code from returning a result close to the reference:
- 论文中的参数是
sigma=16. - 论文使用的不是自然对数，而是二进制的，即幅度演化给出为
2^(L_it).因此，您必须将计算出的指数除以log(2).
- The parameter in the paper is
sigma=16. - The paper uses not the natural logarithm, but the binary one, that is, the magnitude evolution is given as
2^(L_it). So you have to divide the computed exponents bylog(2).
使用https://scicomp.stackexchange.com/questions/中的方法36013/numerical-computation-of-lyapunov-exponent 我得到指数
[2.1531855610566595, -0.00847304754613621, -32.441308372177566]
与参考 (2.16, 0.0, -32.4) 足够接近.
这篇关于已知 ODE 的 Lyapunov 谱 - Python 3的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！