ODEINT - 添加新方程式时的不同结果
ODEINT - Different results when new equations added
希望你一切都好。
这是我的第一个问题，如果有不对的地方，我很抱歉。
我正在研究一些动力系统的数值稳定性和混沌性，更具体地说，是关于由 3 个二阶微分方程组定义的圆形受限三体问题 (CR3BP)。在将这三个二阶微分方程转换为六个一阶微分方程后 as seen here i can finally finally work with them numerically using scipy's ODEINT. Here's an example of an orbit integrated for T = 2^10 with n = 2^18 points (np.linspace(1, 2^10, 2^18)) 这是我的一些代码，要集成的主要函数：
def func(init, t, mu):
#x0, y0, z0, vx0, vy0, vz0 = init
x0, y0, z0, vx0, vy0, vz0 = init[0], init[1], init[2], init[3], init[4], init[5]
dx1dt1 = vx0
dy1dt1 = vy0
dz1dt1 = vz0
dx2dt2 = 2*dy1dt1+d_omega_x(mu, x0, y0, z0)
dy2dt2 = -2*dx1dt1+d_omega_y(mu, x0, y0, z0)
dz2dt2 = d_omega_z(mu, x0, y0, z0)
return np.array([dx1dt1, dy1dt1, dz1dt1, dx2dt2, dy2dt2, dz2dt2])#, dx2dt2, dy2dt2, dz2dt2])
其中 x, y, z, vx, vy, vz = (0.848, 0, 0, 0, 0.0423, 0) 和 mu = 0.01215。
然后是稳定性部分。我正在使用一个名为 Fast Lyapunov Indicator 的混沌检测工具。它基本上由 v'(t)=Df(x)v(t) 定义，其中 Df(x) 是方程组的雅可比矩阵（在本例中为 6x6 矩阵），v(t) 是切向量与 CR3BP 的六个原始方程一起及时演化，然后我取积分 v(t) 的六个分量的范数的 log10 并选择最高值。
人们可能会注意到从 v'(t)=Df(x)v(t) 获得的 6 个“辅助”向量取决于原始的 6 个方程（更具体地说取决于粒子的位置 x、y、z ), 但六个原始方程不依赖于与 v'(t) 定义的切线映射相关的六个新方程和 v(0) 的六个初始条件。
因此，我们不必对 12 个一阶微分方程进行积分。发生的事情是，每次我为 v(0) 设置一个非空初始向量时，对于 CR3BP 的某些初始条件（就像用于生成上述图形的那个）the results obtained are different than those obtained by the integration of only the six original equations 因为系统“崩溃” ”转到 x = y = 0 并且积分器告诉我一些错误而不是“积分成功”，这与应该发生的情况不同。这里，v(0) = (1, 0, 0, 1, 0, 0).
我在两个积分中得到相同结果的唯一情况 is when v(0) = (0, 0, 0, 0, 0, 0)，但我无法获得 Fast Lyapunov 指标的值。
这是新函数的代码片段，其中包含李雅普诺夫指标的六个新方程：
def func2(init, t, mu):
#x0, y0, z0, vx0, vy0, vz0 = init
x0, y0, z0, vx0, vy0, vz0 = init[0], init[1], init[2], init[3], init[4], init[5]
v0, v1, v2, v3, v4, v5 = init[6], init[7], init[8], init[9], init[10], init[11]
#print(init)
dx1dt1 = vx0
dy1dt1 = vy0
dz1dt1 = vz0
dx2dt2 = 2*dy1dt1+d_omega_x(mu, x0, y0, z0)
dy2dt2 = -2*dx1dt1+d_omega_y(mu, x0, y0, z0)
dz2dt2 = d_omega_z(mu, x0, y0, z0)
r1 = r11(mu, x0, y0, z0)
r2 = r22(mu, x0, y0, z0)
jacobiana = [v3,
v4,
v5,
(v0*(mu*(-3*mu - 3*x0 + 3)*(-mu - x0 + 1)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
mu/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(3/2) +
(1 - mu)*(-3*mu - 3*x0)*(-mu - x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2) -
(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(3/2) + 1.0) +
v1*(-3*mu*y0*(-mu - x0 + 1)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
3*y0*(1 - mu)*(-mu - x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v2*(-3*mu*z0*(-mu - x0 + 1)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
3*z0*(1 - mu)*(-mu - x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) + 2*v4),
(v0*(-mu*y0*(-3*mu - 3*x0 + 3)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
y0*(1 - mu)*(-3*mu - 3*x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v1*(3*mu*y0**2/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
mu/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(3/2) +
3*y0**2*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2) -
(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(3/2) + 1.0) +
v2*(3*mu*y0*z0/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) +
3*y0*z0*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) - 2*v3),
(v0*(-mu*z0*(-3*mu - 3*x0 + 3)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
z0*(1 - mu)*(-3*mu - 3*x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v1*(3*mu*y0*z0/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) +
3*y0*z0*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v2*(3*mu*z0**2/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
mu/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(3/2) +
3*z0**2*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2) -
(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(3/2) + 1.0))]
fli = jacobiana
dv1 = fli[0]
dv2 = fli[1]
dv3 = fli[2]
dv4 = fli[3]
dv5 = fli[4]
dv6 = fli[5]
return [dx1dt1, dy1dt1, dz1dt1, dx2dt2, dy2dt2, dz2dt2, dv1, dv2, dv3, dv4, dv5, dv6]
怎么办？这似乎显然是浮点精度的问题，因为每次 运行 代码时我都会得到一些不同的结果。 When I increase the number of points in np.linspace (in this case to 2^20 points), i tend to get correct results，但我不能处理超过一百万个点，而在另一种情况下，我可以用少 4 倍的数据得到正确的结果。我需要对数据应用连续小波变换，因此它变得真的消耗。
再一次，如果问题很长或令人困惑，我很抱歉，如果需要，我可以提供更多相关信息。
总之，非常感谢大家的关注，注意安全
编辑：
正如 Lutz 所指出的，并遵循系统的自然动力学，混沌轨道由 Lyapunov 指标的指数增长值定义——实际上是由范数的 log10 定义的，而不仅仅是范数—— ，结果证明初始向量 v(0) 必须相当小，这样结果才不会溢出。尝试 (1e-10, 0, 0, 0, 0, 0) returns the correct integration. The curve profile for the Lyapunov indicator is also correct, just shifted by a factor log10(1e-10).
再次感谢您！
这可能是由于步长控制也受到快速增长的 v 向量的影响。要么由于刚度而快速减小步长，或者更可能是由于增加步长以匹配主要组件，因此变得不适合原始轨迹的精确集成。这种快速增长是引入 Lyapunov 指数的原因，因为它们以有界的数字捕捉这种增长。
您可以做的是将集成拆分为更小的块，并在每个块的开头规范化 v 向量。人们将不得不试验 v 分量过度控制步长控制需要多长时间。由于耦合是纯乘法的，理论上动态是线性的。因此，如果您将初始 v 缩放为具有标准 1e-100.
，它也会有所帮助
首先检查您使用的错误容限。将它们设置得更窄也有助于稳定计算。通过将最大步长 hmax 设置为外部步长的一半左右，您可能还会取得一些进展。
或者您可以像我在 https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent 中探讨的那样进行 Lyapunov 指数计算。然而，这种方法通过 eigen/singular 向量的 n x n 矩阵和指数乘以时间的 n 乘积增加了维度 n 的系统。
希望你一切都好。
这是我的第一个问题，如果有不对的地方，我很抱歉。
我正在研究一些动力系统的数值稳定性和混沌性，更具体地说，是关于由 3 个二阶微分方程组定义的圆形受限三体问题 (CR3BP)。在将这三个二阶微分方程转换为六个一阶微分方程后 as seen here i can finally finally work with them numerically using scipy's ODEINT. Here's an example of an orbit integrated for T = 2^10 with n = 2^18 points (np.linspace(1, 2^10, 2^18)) 这是我的一些代码，要集成的主要函数：
def func(init, t, mu):
#x0, y0, z0, vx0, vy0, vz0 = init
x0, y0, z0, vx0, vy0, vz0 = init[0], init[1], init[2], init[3], init[4], init[5]
dx1dt1 = vx0
dy1dt1 = vy0
dz1dt1 = vz0
dx2dt2 = 2*dy1dt1+d_omega_x(mu, x0, y0, z0)
dy2dt2 = -2*dx1dt1+d_omega_y(mu, x0, y0, z0)
dz2dt2 = d_omega_z(mu, x0, y0, z0)
return np.array([dx1dt1, dy1dt1, dz1dt1, dx2dt2, dy2dt2, dz2dt2])#, dx2dt2, dy2dt2, dz2dt2])
其中 x, y, z, vx, vy, vz = (0.848, 0, 0, 0, 0.0423, 0) 和 mu = 0.01215。
然后是稳定性部分。我正在使用一个名为 Fast Lyapunov Indicator 的混沌检测工具。它基本上由 v'(t)=Df(x)v(t) 定义，其中 Df(x) 是方程组的雅可比矩阵（在本例中为 6x6 矩阵），v(t) 是切向量与 CR3BP 的六个原始方程一起及时演化，然后我取积分 v(t) 的六个分量的范数的 log10 并选择最高值。
人们可能会注意到从 v'(t)=Df(x)v(t) 获得的 6 个“辅助”向量取决于原始的 6 个方程（更具体地说取决于粒子的位置 x、y、z ), 但六个原始方程不依赖于与 v'(t) 定义的切线映射相关的六个新方程和 v(0) 的六个初始条件。
因此，我们不必对 12 个一阶微分方程进行积分。发生的事情是，每次我为 v(0) 设置一个非空初始向量时，对于 CR3BP 的某些初始条件（就像用于生成上述图形的那个）the results obtained are different than those obtained by the integration of only the six original equations 因为系统“崩溃” ”转到 x = y = 0 并且积分器告诉我一些错误而不是“积分成功”，这与应该发生的情况不同。这里，v(0) = (1, 0, 0, 1, 0, 0).
我在两个积分中得到相同结果的唯一情况 is when v(0) = (0, 0, 0, 0, 0, 0)，但我无法获得 Fast Lyapunov 指标的值。
这是新函数的代码片段，其中包含李雅普诺夫指标的六个新方程：
def func2(init, t, mu):
#x0, y0, z0, vx0, vy0, vz0 = init
x0, y0, z0, vx0, vy0, vz0 = init[0], init[1], init[2], init[3], init[4], init[5]
v0, v1, v2, v3, v4, v5 = init[6], init[7], init[8], init[9], init[10], init[11]
#print(init)
dx1dt1 = vx0
dy1dt1 = vy0
dz1dt1 = vz0
dx2dt2 = 2*dy1dt1+d_omega_x(mu, x0, y0, z0)
dy2dt2 = -2*dx1dt1+d_omega_y(mu, x0, y0, z0)
dz2dt2 = d_omega_z(mu, x0, y0, z0)
r1 = r11(mu, x0, y0, z0)
r2 = r22(mu, x0, y0, z0)
jacobiana = [v3,
v4,
v5,
(v0*(mu*(-3*mu - 3*x0 + 3)*(-mu - x0 + 1)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
mu/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(3/2) +
(1 - mu)*(-3*mu - 3*x0)*(-mu - x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2) -
(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(3/2) + 1.0) +
v1*(-3*mu*y0*(-mu - x0 + 1)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
3*y0*(1 - mu)*(-mu - x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v2*(-3*mu*z0*(-mu - x0 + 1)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
3*z0*(1 - mu)*(-mu - x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) + 2*v4),
(v0*(-mu*y0*(-3*mu - 3*x0 + 3)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
y0*(1 - mu)*(-3*mu - 3*x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v1*(3*mu*y0**2/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
mu/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(3/2) +
3*y0**2*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2) -
(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(3/2) + 1.0) +
v2*(3*mu*y0*z0/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) +
3*y0*z0*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) - 2*v3),
(v0*(-mu*z0*(-3*mu - 3*x0 + 3)/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
z0*(1 - mu)*(-3*mu - 3*x0)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v1*(3*mu*y0*z0/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) +
3*y0*z0*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2)) +
v2*(3*mu*z0**2/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(5/2) -
mu/(y0**2 + z0**2 + (mu + x0 - 1)**2)**(3/2) +
3*z0**2*(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(5/2) -
(1 - mu)/(y0**2 + z0**2 + (mu + x0)**2)**(3/2) + 1.0))]
fli = jacobiana
dv1 = fli[0]
dv2 = fli[1]
dv3 = fli[2]
dv4 = fli[3]
dv5 = fli[4]
dv6 = fli[5]
return [dx1dt1, dy1dt1, dz1dt1, dx2dt2, dy2dt2, dz2dt2, dv1, dv2, dv3, dv4, dv5, dv6]
怎么办？这似乎显然是浮点精度的问题，因为每次 运行 代码时我都会得到一些不同的结果。 When I increase the number of points in np.linspace (in this case to 2^20 points), i tend to get correct results，但我不能处理超过一百万个点，而在另一种情况下，我可以用少 4 倍的数据得到正确的结果。我需要对数据应用连续小波变换，因此它变得真的消耗。
再一次，如果问题很长或令人困惑，我很抱歉，如果需要，我可以提供更多相关信息。
总之，非常感谢大家的关注，注意安全
编辑：
正如 Lutz 所指出的，并遵循系统的自然动力学，混沌轨道由 Lyapunov 指标的指数增长值定义——实际上是由范数的 log10 定义的，而不仅仅是范数—— ，结果证明初始向量 v(0) 必须相当小，这样结果才不会溢出。尝试 (1e-10, 0, 0, 0, 0, 0) returns the correct integration. The curve profile for the Lyapunov indicator is also correct, just shifted by a factor log10(1e-10).
再次感谢您！
这可能是由于步长控制也受到快速增长的 v 向量的影响。要么由于刚度而快速减小步长，或者更可能是由于增加步长以匹配主要组件，因此变得不适合原始轨迹的精确集成。这种快速增长是引入 Lyapunov 指数的原因，因为它们以有界的数字捕捉这种增长。
您可以做的是将集成拆分为更小的块，并在每个块的开头规范化 v 向量。人们将不得不试验 v 分量过度控制步长控制需要多长时间。由于耦合是纯乘法的，理论上动态是线性的。因此，如果您将初始 v 缩放为具有标准 1e-100.
首先检查您使用的错误容限。将它们设置得更窄也有助于稳定计算。通过将最大步长 hmax 设置为外部步长的一半左右，您可能还会取得一些进展。
或者您可以像我在 https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent 中探讨的那样进行 Lyapunov 指数计算。然而，这种方法通过 eigen/singular 向量的 n x n 矩阵和指数乘以时间的 n 乘积增加了维度 n 的系统。