python - ODEINT - 添加新方程时的不同结果
问题描述
希望你们一切都好。
这是我的第一个问题,所以如果有什么不对的地方,我很抱歉。
我正在研究一些动力系统的数值稳定性和混沌性,更具体地说,是关于由 3 个二阶微分方程组定义的圆形受限三体问题 (CR3BP)。在将这三个二阶微分方程转换为六个一阶微分方程之后,我终于可以使用 scipy 的 ODEINT 对它们进行数值处理。这是一个针对 T = 2^10 与 n = 2^18 点 (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 的某些初始条件(就像用于生成上述图形的那个),获得的结果与通过仅对六个原始方程进行积分,因为系统“崩溃”到 x = y = 0 并且积分器告诉我一些错误而不是“积分成功”,这与应该发生的情况不同。这里,v(0) = (1, 0, 0, 1, 0, 0)。
我对两个积分的结果相同的唯一情况是当 v(0) = (0, 0, 0, 0, 0, 0)时,但我无法获得快速 Lyapunov 指标的值。
以下是包含 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]
该怎么办?这显然是浮点精度的问题,因为每次运行代码时都会得到一些不同的结果。当我增加 np.linspace 中的点数(在本例中为 2^20 点)时,我往往会得到正确的结果,但在另一种情况下,我不能处理超过一百万个点可以用少 4 倍的数据得到正确的结果。我需要对数据应用连续小波变换,因此它变得非常消耗。
再一次,如果问题相当长或令人困惑,我很抱歉,如果需要,我可以提供更多关于它的信息。
无论如何,非常感谢您的关注,注意安全。
编辑:
正如 Lutz 所指出的并遵循系统的自然动力学,混沌轨道由 Lyapunov 指标的指数增加值定义——这实际上由范数的 log10 定义,而不仅仅是范数——结果是初始向量 v(0) 必须相当小,这样结果才不会溢出。尝试 (1e-10, 0, 0, 0, 0, 0)返回正确的积分。Lyapunov 指标的曲线轮廓也是正确的,只是移动了一个因子 log10(1e-10)。
再次非常感谢您!
解决方案
这可能是由于步长控制也受到快速增长的v
向量的影响。由于刚度,或者更可能是由于增加步长以匹配主要分量,因此通过快速降低步长,从而变得不适合原始轨迹的精确积分。这种快速增长是引入李雅普诺夫指数的原因,因为它们以非常有界的数字捕捉到这种增长。
您可以做的是将集成拆分为更小的块,并v
在每个块的开始处对向量进行归一化。人们将不得不试验在v
组件过度控制步长控制之前需要多长时间。由于耦合是纯乘法的,理论上动态是线性的。因此,如果您将初始值缩放为v
norm也会有所帮助1e-100
。
但是,首先检查您使用的误差容限。将它们设置得更窄也倾向于稳定计算。您可能还会在将最大步长设置hmax
为外部步长的一半左右方面取得一些进展。
或者您可以像我在https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent中探索的那样进行 Lyapunov 指数计算。然而,这种方法通过特征/奇异向量n
的矩阵和指数乘以时间的乘积来增加维度系统。n x n
n
推荐阅读
- mongodb - MongoDB 提示()失败 - 不确定是否是因为索引仍在索引
- excel - 偶尔会遇到运行时错误 1004
- python - dict_keys 对象不支持索引
- postgresql - Postgres 11 Standby 永远赶不上
- c - 指向数组的初始化指针
- hyperledger-fabric - 一次生成大量资产时出现超时错误
- css - 上下文菜单颜色
- html - HTML5 视频在 macOS Safari 上播放,而不是 iOS Safari
- r - 使用 ggplot2 和 facet_wrap 绘制估计值而不重新拟合模型
- android - 用 kotlin 学习 Android MVVM 架构组件