希望你们都好
这是我的第一个问题,如果有什么不对劲,我很抱歉
我正在研究一些动力系统的数值稳定性和混沌性,更具体地说,是关于由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指示符的值
以下是新函数的代码片段,其中包括六个新的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,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
缩放为norm1e-100
,这也会有所帮助但是,首先检查您使用的误差公差。将它们设置得更窄也有助于稳定计算。如果将最大步长
hmax
设置为外部步长的一半左右,也可能会取得一些进展或者你可以像我在https://scicomp.stackexchange.com/questions/36013/numerical-computation-of-lyapunov-exponent中探讨的那样进行李雅普诺夫指数计算。然而,这种方法通过本征/奇异向量的
n x n
矩阵和指数乘以时间的n
乘积来增加维数n
系统相关问题 更多 >
编程相关推荐