龙格库塔4阶方法

最新推荐文章于 2025-10-29 10:14:11 发布

原创最新推荐文章于 2025-10-29 10:14:11 发布 · 3k 阅读

10 ·

CC 4.0 BY-SA版权

文章标签：

#算法

这段源码展示了4级龙格-库塔(Runge-Kutta 4)方法在卫星轨道数值积分中的应用。通过计算4个中间点的加速度并加权平均，得到最终的位置和速度更新。博客提到了对于密集输出点的处理策略，包括插值方法，并提及了更高阶的Runge-Kutta-Fehlberg方法。此外，讨论了微分方程求解中步长控制和效率问题。

在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方法的系数。