Рунге-Кутты RK4: рабочая лошадка

Урок про метод Рунге-Кутты 4-го порядка (RK4) — золотой стандарт численного решения ОДУ.

Метод RK4 делает шаг, усредняя четыре пробных наклона внутри интервала, и достигает 4-го порядка точности: y_{n+1} = y_n + h/6·(k₁ + 2k₂ + 2k₃ + k₄).

Идея: усреднить наклоны

Слабость Эйлера в том, что он берёт наклон лишь в начале шага и слепо доверяет ему весь шаг. RK4 поступает умнее: он пробует наклон в четырёх точках интервала и усредняет их с продуманными весами. k₁ — наклон в начале (как у Эйлера). k₂ — наклон в середине, если шагнуть туда по k₁. k₃ — снова в середине, но уже по k₂ (уточнение). k₄ — наклон в конце по k₃. Итоговый шаг — взвешенное среднее (k₁ + 2k₂ + 2k₃ + k₄)/6, где середина (она важнее) учтена с удвоенным весом.

def rk4(f, t0, y0, h, n):
    t, y = t0, y0
    точки = [(t, y)]
    for _ in range(n):
        k1 = f(t, y)
        k2 = f(t + h/2, y + h/2 * k1)
        k3 = f(t + h/2, y + h/2 * k2)
        k4 = f(t + h, y + h * k3)
        y = y + h/6 * (k1 + 2*k2 + 2*k3 + k4)   # взвешенное среднее
        t = t + h
        точки.append((t, y))
    return точки

import math
# Та же задача y' = y, y(0)=1, тот же шаг 0.25 — сравним с Эйлером
for t, y in rk4(lambda t, y: y, 0.0, 1.0, 0.25, 4):
    print(f"t={t:.2f}  RK4={y:.8f}  точное={math.exp(t):.8f}  "
          f"ошибка={abs(y - math.exp(t)):.2e}")

Вывод:

t=0.00  RK4=1.00000000  точное=1.00000000  ошибка=0.00e+00
t=0.25  RK4=1.28401693  точное=1.28402542  ошибка=8.49e-06
t=0.50  RK4=1.64869947  точное=1.64872127  ошибка=2.18e-05
t=0.75  RK4=2.11695803  точное=2.11700002  ошибка=4.20e-05
t=1.00  RK4=2.71820994  точное=2.71828183  ошибка=7.19e-05

Сравните с Эйлером из прошлого урока: при том же шаге 0.25 Эйлер ошибался к t=1 на 0.28, а RK4 — на 0.00007. Это в четыре тысячи раз точнее при той же длине шага и всего вчетверо большем числе вычислений функции. Вот почему RK4 — стандарт.

Четвёртый порядок в цифрах

RK4 имеет четвёртый порядок: уменьшение шага вдвое снижает ошибку в 2⁴ = 16 раз. Это родственно тому, как Симпсон относится к трапециям в интегрировании — и не случайно: веса 1, 2, 2, 1 в RK4 перекликаются с симпсоновскими. Четвёртый порядок — удачный компромисс: выше точность на шаг при умеренной сложности (4 вызова f). Методы ещё выше порядка существуют, но для большинства задач RK4 оптимален.

  шаг RK4 = усреднение 4 пробных наклонов:
   k1 ─ наклон в начале
   k2 ─ наклон в середине (шагнув по k1)   } середина с весом 2
   k3 ─ наклон в середине (шагнув по k2)   }
   k4 ─ наклон в конце (шагнув по k3)
   y_new = y + h·(k1 + 2k2 + 2k3 + k4)/6

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

Веса 1/6, 2/6, 2/6, 1/6 и точки выбраны так, чтобы разложение шага RK4 совпало с рядом Тейлора решения вплоть до члена h⁴ — отсюда четвёртый порядок. На практике используют не фиксированный шаг, а адаптивный RK с контролем ошибки: метод Рунге-Кутты-Фельберга (RKF45) или Дорманда-Принса (за ним стоит scipy.integrate.solve_ivp по умолчанию) считают шаг двумя порядками сразу, по их разнице оценивают локальную ошибку и автоматически подбирают h — крупный на гладких участках, мелкий на резких. Важная оговорка: RK4 — явный метод и, как Эйлер, неустойчив на жёстких уравнениях; для них берут неявные методы (BDF, Radau).

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

Считать k₂, k₃, k₄ независимо от предыдущих. Каждый следующий наклон вычисляется в точке, смещённой по предыдущему наклону — порядок вычислений важен.
Применять RK4 к жёсткой задаче. Он явный и взорвётся при крупном шаге; нужны неявные методы.
Брать фиксированный шаг для задачи с резкими участками. Адаптивный RK (RKF45/Дорманд-Принс) тратит шаги экономно и контролирует ошибку.

Итоги

RK4 усредняет 4 пробных наклона за шаг; четвёртый порядок (шаг ÷2 → ошибка ÷16).
При том же шаге RK4 на порядки точнее Эйлера ценой лишь 4 вычислений f вместо одного.
На практике — адаптивный вариант (RKF45, Дорманд-Принс) с автоподбором шага.
RK4 явный ⇒ неустойчив на жёстких уравнениях; там нужны неявные методы.