Численное интегрирование подъёма (метод Эйлера)

Урок строит численную модель подъёма ракеты шаг за шагом — простейшим методом Эйлера.

Метод Эйлера — простейший способ численно решать уравнения движения: на каждом маленьком шаге времени считаем ускорение и обновляем скорость и положение по текущим значениям.

Почему численно

Уравнение Циолковского даёт идеальный $\Delta V$, но не показывает, как меняются скорость и высота во времени с учётом гравитации. Аналитически это сложно. Зато численно — легко: разбиваем полёт на маленькие шаги $\Delta t$ и на каждом считаем силы и обновляем состояние. Это и есть метод Эйлера.

Уравнения шага

На каждом шаге: ускорение $a = \dfrac{F}{m} - g$ (тяга минус тяжесть, вертикально). Затем обновляем скорость $v \mathrel{+}= a\,\Delta t$, высоту $h \mathrel{+}= v\,\Delta t$ и массу $m \mathrel{-}= \dot m\,\Delta t$. Повторяем, пока есть топливо.

import math

g = 9.80665
m0 = 5000.0        # стартовая масса, кг
mdot = 20.0        # расход, кг/с
ve = 2800.0        # скорость истечения, м/с
thrust = mdot * ve # тяга, Н
dt = 1.0           # шаг времени, с
burn_time = 150    # время работы, с

t = 0.0
v = 0.0
h = 0.0
m = m0
grav_loss = 0.0
for _ in range(burn_time):
    a = thrust / m - g
    v += a * dt
    h += v * dt
    grav_loss += g * dt
    m -= mdot * dt
    t += dt

print("Тяга:", round(thrust), "Н")
print("После", burn_time, "с:")
print("  масса:", round(m), "кг")
print("  скорость:", round(v, 1), "м/с")
print("  высота:", round(h / 1000, 1), "км")
print("  гравпотери:", round(grav_loss), "м/с")
dV_ideal = ve * math.log(m0 / m)
print("Идеальный ΔV (Циолковский):", round(dV_ideal), "м/с")
print("Проверка v ≈ ΔV - гравпотери:", round(dV_ideal - grav_loss, 1), "м/с")

Вывод:

Тяга: 56000 Н
После 150 с:
  масса: 2000 кг
  скорость: 1086.2 м/с
  высота: 53.2 км
  гравпотери: 1471 м/с
Идеальный ΔV (Циолковский): 2566 м/с
Проверка v ≈ ΔV - гравпотери: 1095.0 м/с

Смотрите, как красиво сходится: реальная скорость (1086 м/с) почти точно равна идеальному $\Delta V$ минус гравитационные потери (2566 − 1471 = 1095). Формула Циолковского и численная модель согласуются — это подтверждает обе.

Профиль полёта по шагам

Распечатаем состояние каждые 30 секунд, чтобы увидеть, как разгоняется ракета (вначале еле-еле — тяговооружённость чуть больше 1, потом всё быстрее).

import math

g = 9.80665
m0 = 5000.0
mdot = 20.0
ve = 2800.0
thrust = mdot * ve
dt = 1.0

t = 0.0; v = 0.0; h = 0.0; m = m0
print("t, с | v, м/с | h, км | масса, кг")
for step in range(150):
    a = thrust / m - g
    v += a * dt
    h += v * dt
    m -= mdot * dt
    t += dt
    if int(t) % 30 == 0:
        print(int(t), "  |", round(v), "  |", round(h/1000, 1), "  |", round(m))

Вывод:

t, с | v, м/с | h, км | масса, кг
30   | 102   | 1.3   | 4400
60   | 250   | 6.4   | 3800
90   | 446   | 16.8  | 3200
120  | 711   | 34.0  | 2600
150  | 1086   | 53.2  | 2000

Ускорение нарастает: за первые 30 с — лишь 102 м/с, за последние — почти 375 м/с. Причина — ракета легчает, тяговооружённость растёт.

Как работает под капотом

Метод Эйлера прост, но грубоват: он предполагает, что ускорение постоянно внутри шага. Чем меньше $\Delta t$, тем точнее результат, но тем больше вычислений. Профессионалы используют более точные методы (Рунге–Кутты), добавляют сопротивление воздуха, переменное $g$ с высотой и угол тангажа (разворот). Но идея та же: разбить полёт на шаги и интегрировать.

Частые ошибки

Брать слишком большой шаг. При большом $\Delta t$ метод Эйлера заметно врёт; уменьшайте шаг для проверки сходимости.
Забывать обновлять массу. Если масса постоянна, ускорение не нарастает и модель неверна.
Не учитывать гравитацию. Без вычитания $g$ получится идеальный $\Delta V$, а не реальная скорость.

Итоги

Метод Эйлера интегрирует движение маленькими шагами по времени.
На каждом шаге: $a=F/m-g$, обновляем $v$, $h$, $m$.
Численная скорость = идеальный $\Delta V$ минус гравитационные потери — модели сходятся.
Меньший шаг точнее; реальные расчёты добавляют сопротивление и Рунге–Кутту.