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)

score 1 · Accepted Answer

您正在求解的方程是类型

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）为了更接近中间桥接步骤的意图，抽象代码应该阅读（带有(* .. *)表示注释或不必要的计算）

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

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

1 回答 1

Related

Reference