使用 Runge-Kutta 4 计算卫星位置
Posted
技术标签:
【中文标题】使用 Runge-Kutta 4 计算卫星位置【英文标题】：Satellite position computation using Runge-Kutta 4 【发布时间】：2016-12-26 10:00:57 【问题描述】：我的问题与 Runge-Kutta 4 (RK4) 方法和轨道卫星状态向量所需的正确迭代步骤有关。 以下代码（在 Python 中）根据此链接 (http://www.navipedia.net/index.php/GLONASS_Satellite_Coordinates_Computation) 的描述描述了运动：
if total_step_number != 0:
for i in range(1, total_step_number+1):
#Calculate k1
k1[0] = (-cs.GM_GLONASS * XYZ[0] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[0] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ[0]) + (2 * cs.OMEGAE_DOT * XYZDot[1])
k1[1] = (-cs.GM_GLONASS * XYZ[1] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[1] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ[1]) - (2 * cs.OMEGAE_DOT * XYZDot[0])
k1[2] = (-cs.GM_GLONASS * XYZ[2] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[2] * (3 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[2]
#Intermediate step to bridge k1 to k2
XYZ2[0] = XYZ[0] + (XYZDot[0] * h / 2) + (k1[0] * h**2 / 8)
XYZDot2[0] = XYZDot[0] + (k1[0] * h / 2)
XYZ2[1] = XYZ[1] + (XYZDot[1] * h / 2) + (k1[1] * h**2 / 8)
XYZDot2[1] = XYZDot[1] + (k1[1] * h / 2)
XYZ2[2] = XYZ[2] + (XYZDot[2] * h / 2) + (k1[2] * h**2 / 8)
XYZDot2[2] = XYZDot[2] + (k1[2] * h / 2)
radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))
....
还有更多代码，但我想限制我现在展示的内容，因为这是我最感兴趣解决的中间步骤。基本上，对于那些熟悉状态向量并使用 RK4 的人来说，您可以看到在中间步骤中更新了位置和速度，但没有更新加速度。我的问题与更新加速度所需的计算有关。它会开始：
XYZDDot[0] = ...
XYZDDot[1] = ...
XYZDDot[2] = ...
...但是后面的具体内容不是很清楚。欢迎任何建议。
以下是完整代码：
for j in h_step_values:
h = j
if h > 0:
one_way_iteration_steps = one_way_iteration_steps -1
elif h < 0:
one_way_iteration_steps = one_way_iteration_steps +1
XYZ = initial_XYZ
#if total_step_number != 0:
for i in range(0, one_way_iteration_steps):
#Calculate k1
k1[0] = (-cs.GM_GLONASS * XYZ[0] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[0] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ[0]) + (2 * cs.OMEGAE_DOT * XYZDot[1])
k1[1] = (-cs.GM_GLONASS * XYZ[1] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[1] * (1 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ[1]) - (2 * cs.OMEGAE_DOT * XYZDot[0])
k1[2] = (-cs.GM_GLONASS * XYZ[2] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ[2] * (3 - (5*(XYZ[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[2]
#Intermediate step to bridge k1 to k2
XYZ2[0] = XYZ[0] + (XYZDot[0] * h / 2) + (k1[0] * h**2 / 8)
XYZDot2[0] = XYZDot[0] + (k1[0] * h / 2)
XYZDDot2[0] = XYZDDot[0] + (k1[0] * h / 2)
XYZ2[1] = XYZ[1] + (XYZDot[1] * h / 2) + (k1[1] * h**2 / 8)
XYZDot2[1] = XYZDot[1] + (k1[1] * h / 2)
XYZ2[2] = XYZ[2] + (XYZDot[2] * h / 2) + (k1[2] * h**2 / 8)
XYZDot2[2] = XYZDot[2] + (k1[2] * h / 2)
radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))
#Calculate k2
k2[0] = (-cs.GM_GLONASS * XYZ2[0] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[0] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ2[0]) + (2 * cs.OMEGAE_DOT * XYZDot2[1])
k2[1] = (-cs.GM_GLONASS * XYZ2[1] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[1] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ2[1]) - (2 * cs.OMEGAE_DOT * XYZDot2[0])
k2[2] = (-cs.GM_GLONASS * XYZ2[2] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[2] * (3 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[2]
#Intermediate step to bridge k2 to k3
XYZ2[0] = XYZ[0] + (XYZDot[0] * h / 2) + (k2[0] * h**2 / 8)
XYZDot2[0] = XYZDot[0] + (k2[0] * h / 2)
XYZ2[1] = XYZ[1] + (XYZDot[1] * h / 2) + (k2[1] * h**2 / 8)
XYZDot2[1] = XYZDot[1] + (k2[1] * h / 2)
XYZ2[2] = XYZ[2] + (XYZDot[2] * h / 2) + (k2[2] * h**2 / 8)
XYZDot2[2] = XYZDot[2] + (k2[2] * h / 2)
radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))
#Calculate k3
k3[0] = (-cs.GM_GLONASS * XYZ2[0] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[0] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ2[0]) + (2 * cs.OMEGAE_DOT * XYZDot2[1])
k3[1] = (-cs.GM_GLONASS * XYZ2[1] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[1] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ2[1]) - (2 * cs.OMEGAE_DOT * XYZDot2[0])
k3[2] = (-cs.GM_GLONASS * XYZ2[2] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[2] * (3 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[2]
#Intermediate step to bridge k3 to k4
XYZ2[0] = XYZ[0] + (XYZDot[0] * h) + (k3[0] * h**2 / 2)
XYZDot2[0] = XYZDot[0] + (k3[0] * h)
XYZ2[1] = XYZ[1] + (XYZDot[1] * h) + (k3[1] * h**2 / 2)
XYZDot2[1] = XYZDot[1] + (k3[1] * h)
XYZ2[2] = XYZ[2] + (XYZDot[2] * h) + (k3[2] * h**2 / 2)
XYZDot2[2] = XYZDot[2] + (k3[2] * h)
radius = np.sqrt((XYZ2[0]**2)+(XYZ2[1]**2)+(XYZ2[2]**2))
#Calculate k4
k4[0] = (-cs.GM_GLONASS * XYZ2[0] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[0] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[0] + (cs.OMEGAE_DOT**2 * XYZ2[0]) + (2 * cs.OMEGAE_DOT * XYZDot2[1])
k4[1] = (-cs.GM_GLONASS * XYZ2[1] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[1] * (1 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[1] + (cs.OMEGAE_DOT**2 * XYZ2[1]) - (2 * cs.OMEGAE_DOT * XYZDot2[0])
k4[2] = (-cs.GM_GLONASS * XYZ2[2] / radius**3) \
+ ((3/2) * cs.C_20 * cs.GM_GLONASS * cs.SEMI_MAJOR_AXIS_GLONASS**2 * XYZ2[2] * (3 - (5*(XYZ2[2]**2) / (radius**2))) / radius**5) \
+ XYZDDot[2]
for p in range(3):
XYZ[p] = XYZ[p] + XYZDot[p] * h + h**2 * ((k1[p] + 2*k2[p] + 2*k3[p] + k4[p]) / 12)
XYZDot[p] = XYZDot[p] + (h * (k1[p] + 2*k2[p] + 2*k3[p] + k4[p]) / 6)
radius = np.sqrt((XYZ[0])**2 + (XYZ[0])**2 + (XYZ[0])**2)
【问题讨论】：
为什么不（当前步速 - 前一步速度）/time_interval？但是如果加速度没有出现在你的方程中，你为什么要计算它？ @Ev.Kounis，加速度确实出现在方程中（它们被称为 XYZDDot）。与加速度、速度和位置的唯一区别是加速度采用初始值，但此后从未更新。而“#Intermediate step...”部分处理速度和位置的更新。有人告诉我，更新加速几乎没有什么区别，而且会花费更多的处理时间（无论如何可能可以忽略不计）。 但是，目前我在相邻的数据块之间经历了“分裂”：一个向前传播 15 分钟，另一个向后传播 15 分钟。两端的点应该（非常接近）匹配。但是有一个偏移量，我不知道为什么迭代会这样。 加速度不是恒定的吗？毕竟这是空间 @Ev.Kounis。我同意，如果不是恒定的，那么增加可以忽略不计。我只是想看看是否有解决方案，如果有，我可以用它来测试这个未知的差异。 【参考方案1】：你要解的方程是这样的
ddot x = a(x)
其中a(x) 是在k1 计算中计算的加速度。事实上，一阶系统是
dot v = a(x)
dot x = v
RK4 的实现因此开始于
k1 = a(x)
l1 = v
k2 = a(x+l1*h/2) = a(x+v*h/2)
l2 = v+k1*h/2
等等。 l1,l2,... 的使用似乎隐含在代码中，将这些线性组合直接插入它们出现的位置。
简而言之，你没有错过加速度计算，它是代码片段的主要部分。
更新：(8/22)为了更接近中间桥接步骤的意图，抽象代码应为（(* .. *) 表示 cmets 或不必要的计算）
k1 = a(x) (* l1 = v *)
x2 = x + v*h/2 (* v2 = v + k1*h/2 *)
k2 = a(x2) (* l2 = v2 *)
x3 (* = x + l2*h/2 *)
= x + v*h/2 + k1*h^2/4 (* v3 = v + k2*h/2 *)
k3 = a(x3) (* l3 = v3 *)
x4 (* = x + l3*h *)
= x + v*h + k2*h^2/2 (* v4 = v + k3*h *)
k4 = a(x4) (* l4 = v4 *)
delta_v = ( k1+2*(k2+k3)+k4 ) * h/6
delta_x (* = ( l1+2*(l2+l3)+l4 ) * h/6 *)
= v*h + (k1+k2+k3) * h^2/6
【讨论】：
所以从你的回答来看，迭代应该是，例如：XYZDDot[0] = XYZDDot[0] + (k1[0] * h / 2) ? 不，XYZDDOT 是一个测量和传输的外部常数。至少在积分区间内，月球和太阳的引力是恒定的。发生明显变化的是相对于地球的位置，并且为每个中间步骤重新计算所产生的力。
所以实际上加速不需要更新，并且保持原样。如果是这种情况，那么不清楚还有哪些其他因素导致前向积分相差 15 分钟，而后向积分（在下一个星历块上）相差 15 分钟。这两个值应该几乎相同，但它们不是。
添加了完整的代码，所以你可以看看这是否符合你对RK4流程的解释。如您所见，加速度没有任何更新（XYZDDot 的 XYZDDot2）。
向前飞行 15 分钟和向后飞行 15 分钟应该会得到不同的结果。或者是结合起来，以便第二次整合应该回到初始点？ -- 不管怎样，虽然表面上类似，但实现的方法却不是RK4。因此不清楚它的顺序是 4 还是 2。以上是关于使用 Runge-Kutta 4 计算卫星位置的主要内容，如果未能解决你的问题，请参考以下文章