Aplicación del método Runge-Kutta a las EDO de segundo orden

11

¿Cómo puedo reemplazar el método de Euler por el cuarto orden de Runge-Kutta para determinar el movimiento de caída libre en una magnitud gravitacional no constante (por ejemplo, caída libre desde 10 000 km sobre el suelo)?

Hasta ahora escribí integración simple por el método de Euler:

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

x variable significa la posición actual, v significa velocidad, getMagnitude (x) devuelve la aceleración en la posición x.

Intenté implementar RK4:

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

donde el cuerpo de la función rk4 () es:

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

Pero algo está mal, porque me estoy integrando solo una vez usando RK4 (aceleración). La velocidad de integración con RK4 no tiene sentido porque es lo mismo que v * dt.

¿Podría decirme cómo resolver ecuaciones diferenciales de segundo orden utilizando la integración Runge-Kutta? ¿Debo implementar RK4 calculando los coeficientes k1, l1, k2, l2 ... l4? ¿Cómo puedo hacer eso?

Marcin W.
fuente
Hola @Marcin, edité tu título para reflejar mejor cuál creo que es realmente tu problema. Creo que podemos obtener respuestas más útiles y será más fácil de buscar para otros que vean esta pregunta en el futuro con el nuevo título. Siéntete libre de cambiarlo si no estás de acuerdo.
Doug Lipinski

Respuestas:

17

Parece haber bastante confusión acerca de cómo aplicar los métodos de varios pasos (por ejemplo, Runge-Kutta) a las EDO de segundo orden o superiores o los sistemas de EDO. El proceso es muy simple una vez que lo entiendes, pero quizás no sea obvio sin una buena explicación. El siguiente método es el que encuentro más simple.

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))

Desafortunadamente, C ++ no admite de manera nativa operaciones de vectores como esta, por lo que debe usar una biblioteca de vectores, usar bucles o escribir manualmente las partes separadas. En C ++ puedes usar std::valarraypara lograr el mismo efecto. Aquí hay un ejemplo de trabajo simple con aceleración constante.

#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 Lipinski
fuente
66
" Desafortunadamente, C ++ no admite de manera nativa operaciones vectoriales como esta " Creo que sí, incluso en la biblioteca estándar, pero no necesariamente es fácil de usar con otras bibliotecas de álgebra lineal: en.cppreference.com/w/cpp/numeric/valarray , creo Las bibliotecas comunes de álgebra lineal como Eigen, también deben contar como "soporte".
Kirill
1
@ Kirill, gracias por la sugerencia. Todavía soy relativamente nuevo en C ++ y no he usado valarray antes, ¡también aprendí algo útil! Edición para agregar.
Doug Lipinski
1
Quizás este consejo también sea útil: 1) Use el formato clang para formatear automáticamente su código, es realmente estándar y uniforme. 2) Uso typedef std::valarray<double> Vectorpara tipos de uso común. 3) Use en const int NDIM = 2lugar de #definepor tipo de seguridad y corrección. 4) Desde C ++ 11 puede reemplazar el cuerpo de RHS simplemente con return {X[1], 1}. 5) Es realmente poco común en C ++ (a diferencia de C) declarar primero las variables, luego inicializarlas, preferir declarar variables en el mismo lugar donde las inicializa ( double t = 0., etc.)
Kirill
1
@MarcinW. RHS()calcula el lado derecho de la ecuación diferencial. El vector de estado X es (x, v), entonces dX / dt = (dx / dt, dv / dt) = (v, a). Para su problema (si a = G * M / x ^ 2) RHS debería regresar { X[1], G*M/(X[0]*X[0]) }.
Doug Lipinski
1
@Kirill Lo sé, pero eso solo funciona desde C ++ 11, lo que significa que no funciona con las opciones de compilación predeterminadas en los compiladores más populares. Elegí dejar eso a favor de algo que también funciona con los estándares antiguos y, con suerte, reducir la confusión causada por la incapacidad de compilar el código.
Doug Lipinski