Berechnung einer zweidimensionalen Orbitalbahn mit unendlicher Granularität (Nicht-Euler-Integration)

Für zwei Körper ist dies relativ einfach, da die Bewegungsgleichungen einen Kegelschnitt beschreiben (eine Ellipse für eine geschlossene Umlaufbahn, eine Hyperbel für eine „offene“ Umlaufbahn). Sie können die vis viva- Gleichung verwenden, um die Parameter der Umlaufbahn (große Halbachse usw.) aus den gegebenen Anfangsbedingungen zu erhalten, und der Rest folgt.

Bei einer Ellipse kann man die Position auch als Funktion der Zeit ausdrücken als

wo Ihre Herausforderung zu finden ist$a$,$b$Und$\omega$. Dazu würden Sie gerne die Geschwindigkeit und Entfernung am weitesten Punkt wissen, aber Sie können dies im Prinzip für jeden Ort entlang der Umlaufbahn lösen.

Zum Beispiel für einen Satelliten in einer elliptischen Umlaufbahn über einem massiven Planeten (Masse$M$), für gegebene Anfangsgeschwindigkeit$v$in radialer Richtung, in einem Abstand$r$vom Mittelpunkt des Planeten können wir Folgendes wissen (siehe http://en.wikipedia.org/wiki/Vis-viva_equation ):

In diesem Ausdruck$a$ist die große Halbachse. Wir können es aus den gegebenen Anfangsbedingungen berechnen,

Jetzt brauchen wir$b$; wieder, nach der oben zitierten Wikipedia-Seite und für die angegebenen Bedingungen (v senkrecht zum radialen Vektor im Abstand r, also kennen wir den Drehimpuls der Umlaufbahn$L=mvr$), wir können schreiben

Mit$a$Und$b$berechnet, müssen wir nur finden$\omega$was wieder aus den Anfangsbedingungen folgt, da wir wissen

an der Ausgangsposition. Und mit$a$,$b$Und$\omega$bekannt, können Sie in die obigen Gleichungen einsetzen.

Übrigens ist es bei mehr als einem Gravitationsobjekt praktischer, die Integration "richtig" zu machen. Ich empfehle die Verwendung einer Runge-Kutta-Integration vierter Ordnung, die für diese Art von Problem recht stabil ist. Unter http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf finden Sie einige Ideen, wie Sie dies tun können - es ist ein bisschen schwer zu folgen. Ich habe tatsächlich vor einiger Zeit einen Code geschrieben, der dies in Visual Basic tat - er berücksichtigte tatsächlich sowohl die Schwerkraft als auch den atmosphärischen Luftwiderstand. Ich reproduziere es hier ohne Gewähr...

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

Wenn Sie die RKFunktion mehrmals aufrufen, können Sie jetzt die Umlaufbahn verfolgen, und Sie werden feststellen, dass sie ziemlich stabil ist. Sie können problemlos mehrere massive Objekte hinzufügen (aktualisieren Sie einfach die Beschleunigungsformel) und es wird immer noch funktionieren (obwohl die Umlaufbahnen chaotisch werden).

Wenn Sie unbedingt eine Parametergleichung verwenden möchten, finden Sie weitere Informationen unter https://www.mathopenref.com/coordparamellipse.html

Ansonsten gibt es viele verschiedene Integrationsmethoden. Wählen Sie je nach Bedarf zwischen ihnen aus. Unten ist C++-Code für symplektische Integration, Euler-Integration und Runge-Kutta-Integration.

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

... Wo ...

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