龙格库塔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方法的系数。