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

Question

我的问题与 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)

Answer 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) 为了更接近中间桥接步骤的意图，抽象代码应该是 ( with (* .. *) 表示评论或不必要的计算)

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