龍格庫塔4階方法

在groops源碼數值積分時看到orbitPropagatorRungeKutta4.h，源碼如下：

//4級龍格-庫塔方法，求解位置，速度，加速度，并輸出
inline OrbitArc OrbitPropagatorRungeKutta4::integrateArc(OrbitEpoch startEpoch, Time sampling, UInt posCount, ForcesPtr forces,
                                                         SatelliteModelPtr satellite, EarthRotationPtr earthRotation, EphemeridesPtr ephemerides, Bool timing) const
{
  try
  {
    OrbitArc orbit;
    startEpoch.acceleration = acceleration(startEpoch, forces, satellite, earthRotation, ephemerides);
    orbit.push_back(startEpoch);
    const Double dt = sampling.seconds();

    Single::forEach(posCount-1, [&](UInt /*k*/) //[ ]()語句
    {
      // Evaluate accelerations at 4 positions between current epoch and next
      OrbitEpoch k1 = orbit.back();

      OrbitEpoch k2 = k1;
      k2.time        += seconds2time(dt/2.);
      k2.position    += dt/2. * k1.velocity;
      k2.velocity    += dt/2. * k1.acceleration;
      k2.acceleration = acceleration(k2, forces, satellite, earthRotation, ephemerides);

      OrbitEpoch k3 = k1;
      k3.time        += seconds2time(dt/2.);
      k3.position    += dt/2. * k2.velocity;
      k3.velocity    += dt/2. * k2.acceleration;
      k3.acceleration = acceleration(k3, forces, satellite, earthRotation, ephemerides);

      OrbitEpoch k4 = k1;
      k4.time        += sampling;
      k4.position    += dt * k3.velocity;
      k4.velocity    += dt * k3.acceleration;
      k4.acceleration = acceleration(k4, forces, satellite, earthRotation, ephemerides);

      // Compute final value for this epoch
      OrbitEpoch epoch = k1;
      epoch.time        += sampling;
      epoch.position    += (dt/6.) * (k1.velocity + 2*k2.velocity + 2*k3.velocity + k4.velocity);
      epoch.velocity    += (dt/6.) * (k1.acceleration + 2*k2.acceleration + 2*k3.acceleration + k4.acceleration);
      epoch.acceleration = acceleration(epoch, forces, satellite, earthRotation, ephemerides);
      orbit.push_back(epoch);
    }, timing);

    return orbit;
  }
  catch (std::exception &e)
  {
    GROOPS_RETHROW(e)
  }
}
           

可以看出，源碼中采用經典的RK4龍格-庫塔方法，将函數增量通過4個斜率權重平均計算得到《衛星軌道--模型、方法和應用P113》。

 關于經典龍格庫塔法比較好了解的代碼可看：

https://www.cnblogs.com/JustHaveFun-SAN/archive/2012/02/09/2330655.html

Runge-Kutta公式的思路就是利用區間内一些特殊點的一階導數值的線性組合來替代某點處的n階導數值,這樣就可以僅通過一系列一階導數值來得到某點幂級數展開的預測效果.這和泰勒公式正好是反過來的,泰勒公式是用某點的n階幂級數展開來近似得到小領域内的函數值。

關于龍格-庫塔公式原型，以及s級龍格-庫塔公式中，s次函數計算、步長控制等可參考《衛星軌道--模型、方法和應用P114-P117》，滿足以下關系式：

微分方程求解通常需要預定的輸出點，如果連續兩點間的輸出間隔比步長控制設計的步長大的多，就不會有什麼大問題。但是，如果為了達到預定輸出點，步長經常需要截斷，這時使用龍格-庫塔方法就很沒有效率。因為衛星軌道等應用時輸出稠密，有以下兩種方法：第一時用比較寬的時間步長計算微分方程的解，然後用近似多項式插值得到所要的稠密輸出點。另外一種方法是采用連續方法。

這裡給一個Horn插值的6級龍格-庫塔-Fehlberg方法的系數。