Расчет двумерного орбитального пути с бесконечной степенью детализации (интеграция без Эйлера)

Для двух тел это относительно легко, поскольку уравнения движения описывают конику (эллипс для замкнутой орбиты, гиперболу для «открытой» орбиты). Вы можете использовать уравнение живой природы , чтобы получить параметры орбиты (большую полуось и т. д.) из заданных начальных условий, а остальное следует.

Для эллипса вы также можете выразить положение как функцию времени как

где ваша задача найти$a$,$b$и$\omega$. Для этого вам нужно знать скорость и расстояние в самой дальней точке, но в принципе вы можете решить это для любого места на орбите.

Например, для спутника на эллиптической орбите над массивной планетой (масса$M$), для заданной начальной скорости$v$в радиальном направлении, на расстоянии$r$из центра планеты мы можем узнать следующее (см. http://en.wikipedia.org/wiki/Vis-viva_equation ):

В этом выражении$a$является большой полуосью. Мы можем вычислить его из заданных начальных условий,

Теперь нам нужно$b$; снова, следуя странице Википедии, процитированной выше, и для данных условий (v перпендикулярен радиальному вектору на расстоянии r, поэтому мы знаем угловой момент орбиты$L=mvr$), мы можем написать

С$a$и$b$рассчитано, осталось найти$\omega$что снова следует из начальных условий, так как мы знаем

в исходном положении. И с$a$,$b$и$\omega$известные можно подставить в приведенные выше уравнения.

Между прочим, более практично, когда у вас есть более одного гравитационного объекта, интегрировать "правильно". Я рекомендую использовать интеграцию Рунге-Кутты четвертого порядка, которая достаточно стабильна для такого рода задач. Вы можете увидеть http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf для некоторых идей о том, как это сделать - это немного сложно понять. Некоторое время назад я написал код, который делал это на Visual Basic — он фактически учитывал как гравитацию, так и атмосферное сопротивление. Я воспроизвел это здесь без гарантии...

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

Вызывая RKфункцию несколько раз, вы теперь можете отслеживать орбиту, и вы обнаружите, что она достаточно стабильна. Вы можете легко добавить несколько массивных объектов (просто обновите формулу ускорения), и это все равно будет работать (хотя орбиты станут хаотичными).

Если вы абсолютно уверены в использовании параметрического уравнения, см. https://www.mathopenref.com/coordparamellipse.html .

В противном случае существует множество различных методов интеграции. Выбирайте и выбирайте между ними, исходя из ваших потребностей. Ниже приведен код C++ для симплектической интеграции, интеграции Эйлера и интеграции Рунге-Кутты.

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

... где ...

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