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)) 

      .... 

それは中間ステップIだからしかし、私は今のところ示して何を制限したい多くのコードがあります解決に最も関心があります。基本的には、状態ベクトルに慣れていてRK4を使用している人にとっては、位置と速度が中間ステップで更新されるが、加速は更新されないことがわかります。私の質問は、加速をあまりにも更新するために必要な計算に関連しています。それは始まるでしょう：

...しかし、正確に来るものはあまり明確ではありません。アドバイスを歓迎します。

が

 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

あなたが解決されている式はa(x)があなたの​​計算で計算される加速度であるタイプ

ddot x = a(x) 

は次のとおりです。

以下は完全なコードです。実際、一次システムが

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