使用 Runge-Kutta 4 计算卫星位置答案

【问题标题】：Satellite position computation using Runge-Kutta 4使用 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。我同意，如果不是恒定的，那么增加可以忽略不计。我只是想看看是否有解决方案，如果有，我可以用它来测试这个未知的差异。

标签： python runge-kutta orbital-mechanics satellite-navigation

【解决方案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。