ODEINT - 添加新方程时的不同结果

Posted 2023-02-18

【中文标题】ODEINT - 添加新方程时的不同结果

【英文标题】：ODEINT - Different results when new equations added

【发布时间】：2021-08-16 02:30:33

【问题描述】：

希望你一切都好。

这是我的第一个问题，如果有什么不对的地方，我很抱歉。

我正在研究一些动力系统的数值稳定性和混沌性，更具体地说，是关于由 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, 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) 的结果相同，但是我无法获得快速 Lyapunov 指标的值。

以下是新函数的代码 sn-p，其中包含 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).

再次非常感谢您！

【参考方案1】：

这可能是由于步长控制也受到快速增长的v 向量的影响。由于刚度，或者更可能是由于增加步长以匹配主要分量，因此通过快速降低步长，从而变得不适合原始轨迹的精确积分。这种快速增长是引入 Lyapunov 指数的原因，因为它们以非常有界的数字捕捉到这种增长。

您可以做的是将集成拆分为更小的块，并在每个块的开头规范化 v 向量。人们将不得不试验在v 组件过度控制步长控制之前需要多长时间。由于耦合是纯乘法的，理论上动态是线性的。因此，如果您将初始 v 缩放为标准 1e-100，它也会有所帮助。

不过，首先检查您使用的误差容限。将它们设置得更窄也倾向于稳定计算。在将最大步长 hmax 设置为外部步长的一半左右时，您可能还会取得一些进展。

或者您可以像我在https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent 中探讨的那样进行李雅普诺夫指数计算。然而，这种方法通过特征/奇异向量的n x n 矩阵和指数乘以时间的n 乘积来增加维度系统n。

【讨论】：

您好，Lutz，非常感谢您的澄清。事实证明，这是系统动力学的“自然”过程。由于这个问题往往发生在混沌轨道上，而且对于混沌轨道，FLI 往往具有指数分布——我忘了​​更正，FLI 考虑的是范数的 log10，而不仅仅是范数——它完全有意义溢出到倾向于呈​​指数增长的值。我通过为 v(0) 选择较小的值，然后将结果归一化以开始和 y = 0 来纠正问题。非常感谢，再次！