Применение метода Рунге-Кутты к ОДУ второго порядка

11

Как я могу заменить метод Эйлера на 4-й порядок Рунге-Кутты, чтобы определить движение свободного падения с не постоянной гравитационной величиной (например, свободное падение с высоты 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.

Не могли бы вы рассказать, как решать дифференциальные уравнения второго порядка с помощью интегрирования Рунге-Кутты? Должен ли я реализовать RK4 путем вычисления коэффициентов k1, l1, k2, l2 ... l4? Как я могу это сделать?

Марчин В.
источник
Привет, @Marcin, я отредактировал твой заголовок, чтобы лучше отразить то, что я думаю о твоей проблеме. Я думаю, что мы можем получить больше полезных ответов, и это будет более доступным для поиска для тех, кто увидит этот вопрос в будущем с новым названием. Не стесняйтесь изменить его обратно, если вы не согласны.
Дуг Липински

Ответы:

17

Похоже, существует некоторая путаница в отношении того, как применять многошаговые (например, Рунге-Кутта) методы к ODE 2-го или более высокого порядка или системам ODE. Процесс очень прост, если вы понимаете это, но, возможно, не очевиден без хорошего объяснения. Следующий метод, который я считаю самым простым.

F=mx¨

[x˙v˙]=[vF/m]

vxk1k4(x,v)

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 )(x˙(t),v˙(t))

К сожалению, 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;

}
Дуг Липински
источник
6
« К сожалению, C ++ изначально не поддерживает векторные операции, подобные этой », я думаю, что это возможно, даже в стандартной библиотеке, но не обязательно легко использовать с другими библиотеками линейной алгебры: en.cppreference.com/w/cpp/numeric/valarray Я думаю, Общие библиотеки линейной алгебры, такие как Eigen, также должны учитываться как «поддержка».
Кирилл
1
@Kirill, спасибо за совет. Я все еще относительно новичок в C ++, и раньше я не использовал valarray, я тоже научился чему-то полезному! Редактирование добавить.
Дуг Липински
1
Возможно, этот совет также будет полезен: 1) Используйте формат clang для автоматического форматирования кода, он действительно стандартный и единый. 2) Используйте typedef std::valarray<double> Vectorдля часто используемых типов. 3) Используйте const int NDIM = 2вместо #defineдля безопасности типа и правильности. 4) Начиная с C ++ 11, вы можете просто заменить тело RHS return {X[1], 1}. 5) В C ++ очень редко (в отличие от C) сначала объявляют переменные, а затем инициализируют их, предпочитают объявлять переменные в том же месте, где вы их инициализируете ( double t = 0.и т. Д.)
Кирилл
1
@MarcinW. RHS()вычисляет правую часть дифференциального уравнения. Вектор состояния X имеет вид (x, v), поэтому dX / dt = (dx / dt, dv / dt) = (v, a). Для вашей проблемы (если a = G * M / x ^ 2) RHS должен вернуться { X[1], G*M/(X[0]*X[0]) }.
Даг Липински
1
@Kirill Я знаю, но это работает только с C ++ 11, что означает, что он не работает с опциями компилятора по умолчанию на самых популярных компиляторах. Я решил оставить это в пользу чего-то, что также работает со старыми стандартами и, надеюсь, уменьшит путаницу, вызванную невозможностью компилировать код.
Даг Липински