如何解决ODEINT - 添加新方程时的不同结果
希望你一切都好。
这是我的第一个问题，如果有什么不对的地方，我很抱歉。
我正在研究一些动力系统的数值稳定性和混沌性，更具体地说，是关于由 3 个二阶微分方程组定义的圆形受限三体问题 (CR3BP)。在将这三个二阶微分方程转换为六个一阶微分方程 as seen here 后，我终于可以使用 scipy 的 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,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,z0)
dy2dt2 = -2*dx1dt1+d_omega_y(mu,z0)
dz2dt2 = d_omega_z(mu,z0)
return np.array([dx1dt1,dy1dt1,dz1dt1,dx2dt2,dy2dt2,dz2dt2])#,dz2dt2])
其中 x,y,z,vx,vy,vz = (0.848,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,1,0).
我对两个积分 is when v(0) = (0,0) 的结果相同的唯一情况，但是我无法获得快速李雅普诺夫指标的值。
def func2(init,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,z0)
r1 = r11(mu,z0)
r2 = r22(mu,z0)
jacobiana = [v3,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,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，但我无法处理超过 100 万个点，而对于另一种情况，我可以用少 4 倍的数据获得正确的结果。我需要对数据应用连续小波变换，因此它变得真正消耗。
再次，如果问题很长或令人困惑，我很抱歉，如果需要，我可以提供更多信息。
无论如何，非常感谢您的关注，注意安全。
编辑：
正如 Lutz 所指出的，并遵循系统的自然动力学，其中混沌轨道由 Lyapunov 指标的指数增加值定义 - 实际上由范数的 log10 定义，不仅是范数——，它结果证明初始向量 v(0) 必须相当小，这样结果就不会溢出。尝试 (1e-10,0) returns the correct integration. The curve profile for the Lyapunov indicator is also correct,just shifted by a factor log10(1e-10).
再次感谢您！
解决方法
这可能是由于步长控制也受到快速增长的 v 向量的影响。通过由于刚度而迅速减小步长，或者更有可能是由于增加步长以匹配主要分量，因此变得不适合原始轨迹的精确积分。这种快速增长是引入李雅普诺夫指数的原因，因为它们以很好的有界数捕捉这种增长。
您可以做的是将集成拆分为更小的块，并在每个块的开头对 v 向量进行归一化。人们必须试验 v 组件过度支配步长控制需要多长时间。由于耦合是纯粹的乘法，因此理论上动态是线性的。因此，如果您将初始 v 缩放为规范 1e-100，它也会有所帮助。
首先检查您使用的容错。将它们设置得更窄也有助于稳定计算。将最大步长 hmax 设置为外部步长的一半左右，您也可能会取得一些进展。
或者您可以像我在 https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent 中探讨的那样进行李雅普诺夫指数计算。然而，这种方法通过本征/奇异向量的 n 矩阵和指数乘以时间的 n x n 乘积来增加维度 n 的系统。
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
