Cálculo de una trayectoria orbital bidimensional con granularidad infinita (integración no Euler)

Para dos cuerpos esto es relativamente fácil ya que las ecuaciones de movimiento describen una cónica (una elipse para una órbita cerrada, una hipérbola para una órbita "abierta"). Puede usar la ecuación de vis viva para obtener los parámetros de la órbita (semieje mayor, etc.) a partir de las condiciones iniciales dadas, y el resto sigue.

Para una elipse, también puede expresar la posición en función del tiempo como

donde tu reto es encontrar$a$,$b$y$\omega$. Para esto, le gustaría saber la velocidad y la distancia en el punto más lejano, pero en principio puede resolver esto para cualquier lugar a lo largo de la órbita.

Por ejemplo, para un satélite en órbita elíptica sobre un planeta masivo (masa$M$), para una velocidad inicial dada$v$en la dirección radial, a una distancia$r$desde el centro del planeta, podemos saber lo siguiente (ver http://en.wikipedia.org/wiki/Vis-viva_equation ):

En esta expresión,$a$es el semieje mayor. Podemos calcularlo a partir de las condiciones iniciales dadas,

ahora necesitamos$b$; nuevamente, siguiendo la página de wikipedia citada anteriormente, y para las condiciones dadas (v perpendicular al vector radial a la distancia r, por lo que conocemos el momento angular de la órbita$L=mvr$), podemos escribir

Con$a$y$b$calculado, sólo tenemos que encontrar$\omega$que se sigue de nuevo de las condiciones iniciales, ya que sabemos

en la posición inicial. Y con$a$,$b$y$\omega$conocido se puede sustituir en las ecuaciones anteriores.

Por cierto, es más práctico, cuando tienes más de un objeto gravitacional, hacer la integración "correctamente". Recomiendo usar una integración de Runge-Kutta de cuarto orden, que es bastante estable para este tipo de problema. Puede ver http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf para obtener algunas ideas sobre cómo hacer esto; es un poco difícil de seguir. De hecho, escribí un fragmento de código que hizo esto en Visual Basic hace un tiempo; en realidad, estaba teniendo en cuenta tanto la gravedad como la resistencia atmosférica. Lo reproduzco aquí sin garantía...

Option Explicit

Const gm = 9.8 * 6300000# * 6300000# ' approximate value of GM on earth surface
                                     ' the # signifies a double precision "integer"
Const C_d = 0.5 * 0.03 * 0.03 * 0.2  ' drag coefficient for particular object
Const R_e = 6300000#                 ' approximate radius of earth
Const mass = 25                      ' mass of object
Const k1 = 12000#                    ' constants used to fit density of atmosphere
Const k2 = 22000#                    ' again, approximate!

Function rho(h)   
  ' return density as a function of height
  If h < 0 Then rho = 1.2 Else rho = 1.225 * Exp(-(h / k1 + (h / k2) ^ (3 / 2)))
End Function

Sub accel(x, y, vx, vy, ByRef ax, ByRef ay)
' compute acceleration of object at location x,y with velocity vx, vy
' return values in ax, ay
  Dim r As Double, v2 As Double, v As Double, r3 As Double, h As Double
  Dim density
  r = Sqr(x * x + y * y) ' sqrt() in most languages... this is VBA
  v2 = vx * vx + vy * vy
  v = Sqr(v2)
  r3 = 1# / r ^ 3
  density = rho(r - R_e)
  ' don't need the second term if you don't care about drag!
  ax = -gm * x * r3 - vx ^ 3 * C_d * density / (v * mass)
  ay = -gm * y * r3 - vy ^ 3 * C_d * density / (v * mass)
End Sub

Function rk(ByRef x, ByRef y, ByRef vx, ByRef vy, dt)
' implement one Runge-Kutta fourth order stop
' this requires calculating the acceleration at verious locations
' for a single "step" in the algorithm
  Dim ax As Double, ay As Double
  Dim kx1, kx2, kx3, kx4, ky1, ky2, ky3, ky4, lx1, lx2, lx3, lx4, ly1, ly2, ly3, ly4, dt2
  ' get acceleration at initial point
  accel x, y, vx, vy, ax, ay
  ' half time step:
  dt2 = dt / 2
  kx1 = dt2 * ax
  ky1 = dt2 * ay
  lx1 = dt2 * vx
  ly1 = dt2 * vy
  ' get acceleration at new location
  accel x + lx1, y + ly1, vx + kx1, vy + ky1, ax, ay
  kx2 = dt2 * ax
  ky2 = dt2 * ay
  lx2 = dt2 * (vx + kx1)
  ly2 = dt2 * (vy + ky1)
  ' get acceleration at half way point:
  accel x + lx2, y + ly2, vx + kx2, vy + ky2, ax, ay
  kx3 = dt * ax
  ky3 = dt * ay
  lx3 = dt * (vx + kx2)
  ly3 = dt * (vy + ky2)
  ' compute acceleration for third combination of position and velocity:
  accel x + lx3, y + ly3, vx + kx3, vy + ky3, ax, ay
  kx4 = dt2 * ax
  ky4 = dt2 * ay
  lx4 = dt2 * (vx + kx3)
  ly4 = dt2 * (vy + ky3)
  ' compute best approximation to new x, y, vx, vy
  x = x + (lx1 + 2# * lx2 + lx3 + lx4) / 3#
  y = y + (ly1 + 2# * ly2 + ly3 + ly4) / 3#
  vx = vx + (kx1 + 2# * kx2 + kx3 + kx4) / 3#
  vy = vy + (ky1 + 2# * ky2 + ky3 + ky4) / 3#

End Function

Al llamar a la RKfunción varias veces, ahora puede rastrear la órbita y encontrará que es bastante estable. Puede agregar fácilmente múltiples objetos masivos (solo actualice la fórmula de aceleración) y seguirá funcionando (aunque las órbitas se volverán caóticas).

Si está absolutamente decidido a usar una ecuación paramétrica, consulte https://www.mathopenref.com/coordparamellipse.html

De lo contrario, hay muchos métodos de integración diferentes. Escoja y elija entre ellos, en función de sus necesidades. A continuación se muestra el código C++ para la integración simpléctica, la integración de Euler y la integración de Runge-Kutta.

// https://en.wikipedia.org/wiki/Symplectic_integrator
// Also known as Verlet integration
void proceed_symplectic2(custom_math::vector_3 &pos, custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
    static const double c[2] = { 0, 1 };
    static const double d[2] = { 0.5, 0.5 };

//  pos += vel * c[0] * dt; // first element c[0] is always 0
    vel += acceleration(pos, vel, ang_vel) * d[0] * dt;

    pos += vel * c[1] * dt;
    vel += acceleration(pos, vel, ang_vel) * d[1] * dt;
}

// https://en.wikipedia.org/wiki/Symplectic_integrator
void proceed_symplectic4(custom_math::vector_3 &pos, custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
    static double const cr2 = pow(2.0, 1.0 / 3.0);

    static const double c[4] =
    {
        1.0 / (2.0*(2.0 - cr2)),
        (1.0 - cr2) / (2.0*(2.0 - cr2)),
        (1.0 - cr2) / (2.0*(2.0 - cr2)),
        1.0 / (2.0*(2.0 - cr2))
    };

    static const double d[4] =
    {
        1.0 / (2.0 - cr2),
        -cr2 / (2.0 - cr2),
        1.0 / (2.0 - cr2),
        0.0
    };

    pos += vel * c[0] * dt;
    vel += acceleration(pos, vel, ang_vel) * d[0] * dt;

    pos += vel * c[1] * dt;
    vel += acceleration(pos, vel, ang_vel) * d[1] * dt;

    pos += vel * c[2] * dt;
    vel += acceleration(pos, vel, ang_vel) * d[2] * dt;

    pos += vel * c[3] * dt;
//  vel += acceleration(pos, vel, ang_vel) * d[3] * dt; // last element d[3] is always 0
}

void proceed_Euler(custom_math::vector_3 &pos, custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
    vel += acceleration(pos, vel, ang_vel) * dt;
    pos += vel * dt;
}

inline void proceed_RK2(custom_math::vector_3 &pos, custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
    custom_math::vector_3 k1_velocity = vel;
    custom_math::vector_3 k1_acceleration = acceleration(pos, k1_velocity, ang_vel);
    custom_math::vector_3 k2_velocity = vel + k1_acceleration * dt*0.5;
    custom_math::vector_3 k2_acceleration = acceleration(pos + k1_velocity * dt*0.5, k2_velocity, ang_vel);

    vel += k2_acceleration * dt;
    pos += k2_velocity * dt;
}

void proceed_RK4(custom_math::vector_3 &pos, custom_math::vector_3 &vel, const custom_math::vector_3 &ang_vel)
{
    static const double one_sixth = 1.0 / 6.0;

    custom_math::vector_3 k1_velocity = vel;
    custom_math::vector_3 k1_acceleration = acceleration(pos, k1_velocity, ang_vel);
    custom_math::vector_3 k2_velocity = vel + k1_acceleration * dt*0.5;
    custom_math::vector_3 k2_acceleration = acceleration(pos + k1_velocity * dt*0.5, k2_velocity, ang_vel);
    custom_math::vector_3 k3_velocity = vel + k2_acceleration * dt*0.5;
    custom_math::vector_3 k3_acceleration = acceleration(pos + k2_velocity * dt*0.5, k3_velocity, ang_vel);
    custom_math::vector_3 k4_velocity = vel + k3_acceleration * dt;
    custom_math::vector_3 k4_acceleration = acceleration(pos + k3_velocity * dt, k4_velocity, ang_vel);

    vel += (k1_acceleration + (k2_acceleration + k3_acceleration)*2.0 + k4_acceleration)*one_sixth*dt;
    pos += (k1_velocity + (k2_velocity + k3_velocity)*2.0 + k4_velocity)*one_sixth*dt;
}

... dónde ...

custom_math::vector_3 acceleration(const custom_math::vector_3 &pos, const custom_math::vector_3 &vel, custom_math::vector_3 &ang_vel)
{
    custom_math::vector_3 grav_dir = sun_pos - pos;

    float distance = grav_dir.length();

    grav_dir.normalize();
    custom_math::vector_3 accel = grav_dir * (G*sun_mass / pow(distance, 2.0));

    return accel;
}