Zastosowanie metody Runge-Kutta do ODE drugiego rzędu

11

Jak mogę zastąpić metodę Eulera czwartym rzędem Runge-Kutty, aby określić ruch swobodnego spadku o nie stałej wielkości grawitacyjnej (np. Swobodny spadek z 10 000 km nad ziemią)?

Do tej pory napisałem prostą integrację metodą Eulera:

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

x zmienna oznacza aktualną pozycję, v oznacza prędkość, getMagnitude (x) zwraca przyspieszenie na pozycji x.

Próbowałem zaimplementować RK4:

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

gdzie ciało funkcji rk4 () to:

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;
}

Ale coś jest nie tak, ponieważ integruję tylko raz używając RK4 (przyspieszenie). Całkowanie prędkości za pomocą RK4 nie ma sensu, ponieważ jest takie samo jak v * dt.

Czy możesz mi powiedzieć, jak rozwiązać równania różniczkowe drugiego rzędu za pomocą integracji Runge-Kutta? Czy powinienem zaimplementować RK4 obliczając współczynniki k1, l1, k2, l2 ... l4? Jak mogę to zrobić?

Marcin W.
źródło
Cześć @Marcin, zredagowałem twój tytuł, aby lepiej odzwierciedlić to, co według mnie jest twoim problemem. Wydaje mi się, że możemy uzyskać bardziej przydatne odpowiedzi i łatwiej będzie je wyszukać innym, którzy zobaczą to pytanie w przyszłości z nowym tytułem. Jeśli się nie zgadzasz, możesz to zmienić.
Doug Lipiński

Odpowiedzi:

17

Wydaje się, że istnieje sporo zamieszania co do sposobu stosowania metod wieloetapowych (np. Runge-Kutta) do ODE drugiego lub wyższego rzędu lub systemów ODE. Proces jest bardzo prosty, gdy go zrozumiesz, ale być może nie jest oczywisty bez dobrego wyjaśnienia. Poniższa metoda jest najprostsza.

W twoim przypadku równanie różniczkowe, które chciałbyś rozwiązać, to . Pierwszym krokiem jest napisanie ODE drugiego rzędu jako systemu ODE pierwszego rzędu. Odbywa się to jakoF=mx¨

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

Wszystkie równania w tym systemie musi być rozwiązany równocześnie, co znaczy, że należy nie advance , a następnie postęp , powinny być rozszerzone zarówno w tym samym czasie. W językach, które obsługują operacje wektorowe bez pętli, można to łatwo zrobić, tworząc wszystkie niezbędne terminy w wektorach kodu o długości 2. Funkcja obliczająca prawą stronę (szybkość zmian) ODE powinna zwrócić wektor o długości 2 , aby powinny być wektory o długości 2, a zmienna stanu powinien być wektor o długości 2. W MATLAB kod niezbędny do intensyfikacji czasu można zapisać jako:x ( x , v )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

gdzie i zwraca wektor zawierający . Jak widać, poprzez wektoryzację rzeczy nie musisz nawet zmieniać składni kodu RK4 bez względu na to, ile równań znajduje się w twoim systemie ODE.X( ˙ x ( t ) , ˙ v ( t ) )=(x,v)RHS( t, X )(x˙(t),v˙(t))

Niestety C ++ nie obsługuje natywnie takich operacji wektorowych, więc musisz albo użyć biblioteki wektorów, użyć pętli, albo ręcznie napisać osobne części. W C ++ możesz użyć std::valarraytego samego efektu. Oto prosty działający przykład ze stałym przyspieszeniem.

#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;

}
Doug Lipiński
źródło
6
Niestety, C ++ nie natywnie wsparcie operacji wektorowej tak ” Myślę, że to robi, w bibliotece standardowej nawet, ale niekoniecznie łatwe do wykorzystania z innych bibliotek algebry liniowej: en.cppreference.com/w/cpp/numeric/valarray myślę popularne biblioteki algebry liniowej, takie jak Eigen, również powinny być liczone jako „wsparcie”.
Kirill
1
@Kirill, dzięki za wskazówkę. Nadal jestem stosunkowo nowy w C ++ i wcześniej nie korzystałem z Valarray, właśnie nauczyłem się czegoś przydatnego! Edycja do dodania.
Doug Lipiński
1
Być może ta rada będzie również pomocna: 1) Użyj formatu clang do automatycznego formatowania kodu, to naprawdę standardowy i jednolity. 2) Użyj typedef std::valarray<double> Vectordla często używanych typów. 3) Użyj const int NDIM = 2zamiast #definedla bezpieczeństwa typu i poprawności. 4) Od wersji C ++ 11 można po prostu zastąpić korpus RHS return {X[1], 1}. 5) W C ++ (w przeciwieństwie do C) naprawdę rzadko zdarza się, aby najpierw deklarować zmienne, a następnie inicjować je, preferować deklarowanie zmiennych w tym samym miejscu, w którym je inicjujesz ( double t = 0.itp.)
Kirill
1
@MarcinW. RHS()oblicza prawą stronę równania różniczkowego. Wektor stanu X to (x, v), więc dX / dt = (dx / dt, dv / dt) = (v, a). Dla twojego problemu (jeśli a = G * M / x ^ 2) RHS powinien powrócić { X[1], G*M/(X[0]*X[0]) }.
Doug Lipiński
1
@Kirill Wiem, ale to działa tylko od C ++ 11, co oznacza, że ​​nie działa z domyślnymi opcjami kompilatora w najpopularniejszych kompilatorach. Zdecydowałem się to pominąć na rzecz czegoś, co działa również ze starymi standardami i mam nadzieję, że zmniejszy zamieszanie spowodowane niemożnością skompilowania kodu.
Doug Lipiński