计算科学 - 将 Runge-Kutta 方法应用于二阶 ODE - 吾爱随笔录

将 Runge-Kutta 方法应用于二阶 ODE

计算科学数值分析 C++ 龙格库塔积分方程

2021-12-16 03:33:26

如何用 Runge-Kutta 4 阶替换 Euler 方法来确定非恒定重力大小的自由落体运动（例如，从地面以上 10 000 公里的自由落体）？

到目前为止，我通过欧拉方法编写了简单的积分：

while()
{
    v += getMagnitude(x) * dt;
    x += v * dt;
    time += dt;
}

x 变量表示当前位置，v 表示速度，getMagnitude(x) 返回 x 位置的加速度。

我尝试实现 RK4：

while()
{
    v += rk4(x, dt) * dt; // rk4() instead of getMagintude()
    x += v * dt;
    time += dt;
}

其中 rk4() 函数体是：

inline double rk4(double tx, double tdt)
{
   double k1 = getMagnitude(tx);
   double k2 = getMagnitude(tx + 0.5 * tdt * k1);
   double k3 = getMagnitude(tx + 0.5 * tdt * k2);
   double k4 = getMagnitude(tx + tdt * k3);

   return (k1 + 2*k2 + 2*k3 + k4)/6.0;
}

但是出了点问题，因为我只使用 RK4（加速）集成了一次。使用 RK4 积分速度没有意义，因为它与 v * dt 相同。

你能告诉我如何使用龙格-库塔积分求解二阶微分方程吗？我应该通过计算 k1、l1、k2、l2 ... l4 系数来实现 RK4 吗？我怎样才能做到这一点？

1个回答

关于如何将多步（例如 Runge-Kutta）方法应用于二阶或更高阶 ODE 或 ODE 系统，似乎存在相当多的困惑。一旦你理解了这个过程就非常简单，但如果没有很好的解释，它可能并不明显。下面的方法是我觉得最简单的一种。

在您的情况下，您要解决的微分方程是 $F = m\ddot{x}$ . 第一步是将这个二阶 ODE 编写为一阶 ODE 系统。这是这样做的

[\begin{matrix} \dot{x} \\ \dot{v} \end{matrix}] = [\begin{matrix} v \\ F / m \end{matrix}]

$\left[\begin{array}{c} \dot{x} \\ \dot{v} \end{array}\right] = \left[\begin{array}{c} v \\ F/m \end{array}\right]$

这个系统中的所有方程必须同时求解，也就是说你不应该前进 $v$ 然后前进 $x$ ，它们应该同时被推进。在支持不带循环的向量运算的语言中，这很容易通过在长度为 2 的代码向量中添加所有必要项来完成。计算 ODE 右侧（变化率）的函数应该返回长度为 2 的向量, k1tok4应该是长度为 2 的向量，以及你的状态变量 $(x,v)$ 应该是一个长度为 2 的向量。在 MATLAB 中，时间步长的必要代码可以写成：

while (t<TMAX)
    k1 = RHS( t, X );
    k2 = RHS( t + dt / 2, X + dt / 2 * k1 );
    k3 = RHS( t + dt / 2, X + dt / 2 * k2 );
    k4 = RHS( t + dt, X + dt * k3 );
    X = X + dt / 6 * ( k1 + 2 * k2 + 2 * k3 + k4 );
    t = t + dt;
end

在哪里X $=(x,v)$ 并RHS( t, X )返回一个包含 $(\dot{x}(t),\dot{v}(t))$ . 如您所见，通过向量化事物，您甚至不需要更改 RK4 代码的语法，无论您的 ODE 系统中有多少方程。

~~不幸的是，C++ 本身并不支持这样的向量操作，因此您需要使用向量库、使用循环或手动写出单独的部分。~~在 C++ 中你可以使用std::valarray来达到同样的效果。这是一个恒定加速度的简单工作示例。

#include <valarray>
#include <iostream>

const size_t NDIM = 2;

typedef std::valarray<double> Vector;

Vector RHS( const double t, const Vector X )
{
  // Right hand side of the ODE to solve, in this case:
  // d/dt(x) = v;
  // d/dt(v) = 1;
  Vector output(NDIM);
  output[0] = X[1];
  output[1] = 1;
  return output;
}

int main()
{

  //initialize values

  // State variable X is [position, velocity]
  double init[] = { 0., 0. };
  Vector X( init, NDIM );

  double t = 0.;
  double tMax=5.;
  double dt = 0.1;

  //time loop
  int nSteps = round( ( tMax - t ) / dt );
  for (int stepNumber = 1; stepNumber<=nSteps; ++stepNumber)
  {

    Vector k1 = RHS( t, X );
    Vector k2 = RHS( t + dt / 2.0,  X + dt / 2.0 * k1 );
    Vector k3 = RHS( t + dt / 2.0, X + dt / 2.0 * k2 );
    Vector k4 = RHS( t + dt, X + dt * k3 );

    X += dt / 6.0 * ( k1 + 2.0 * k2 + 2.0 * k3 + k4 );
    t += dt;
  }
  std::cout<<"Final time: "<<t<<std::endl;
  std::cout<<"Final position: "<<X[0]<<std::endl;
  std::cout<<"Final velocity: "<<X[1]<<std::endl;

}

其它你可能感兴趣的问题

上一篇有限体积代码的数据结构：数组与类下一篇使用 petsc4py 与 c/c++/fortran 的效率