Классический RK4: вывод коэффициентов
RK4 — самый известный метод семейства. Он берёт четыре наклона: один в начале, два в середине и один в конце, а собирает их в формулу с загадочными на первый взгляд весами 1, 2, 2, 1. Разберёмся, откуда они и почему дают именно четвёртый порядок.
Классический RK4 — явный метод Рунге-Кутты четвёртого порядка с четырьмя наклонами и весами (1, 2, 2, 1)/6, для которого локальная ошибка одного шага есть O(h⁵), а накопленная за отрезок — O(h⁴).
Запишем формулу целиком, а потом расшифруем каждый кусочек:
k₁ = f(tₙ, yₙ) наклон в НАЧАЛЕ k₂ = f(tₙ + h/2, yₙ + (h/2)·k₁) наклон в СЕРЕДИНЕ, по k₁ k₃ = f(tₙ + h/2, yₙ + (h/2)·k₂) наклон в СЕРЕДИНЕ, по k₂ k₄ = f(tₙ + h, yₙ + h·k₃) наклон в КОНЦЕ, по k₃ yₙ₊₁ = yₙ + (h/6)·(k₁ + 2·k₂ + 2·k₃ + k₄)
Четыре наклона по очереди
Прочитаем формулу как маршрут разведчика. k₁ — это наклон прямо здесь и сейчас, в левом конце шага; то же, что считает Эйлер. k₂ — мы делаем пробный полушаг по k₁ до середины отрезка и меряем наклон там. k₃ — снова идём в середину, но на этот раз отталкиваемся уже от улучшенного наклона k₂; это уточнённая оценка наклона в центре. k₄ — полным шагом по k₃ добираемся до правого конца и меряем наклон там. Итог — взвешенное среднее этих четырёх замеров. Середина учитывается вдвое сильнее краёв, потому что на неё приходятся две независимые оценки (k₂ и k₃), и именно поведение в середине лучше всего характеризует весь отрезок.
Откуда веса 1, 2, 2, 1
Связь с правилом Симпсона
Самый наглядный способ понять веса — рассмотреть случай, когда f зависит только от t: y' = g(t). Тогда решение — это просто интеграл от g, и наш шаг превращается в численное интегрирование на отрезке [tₙ, tₙ+h]. Подставим: k₁ = g(tₙ), k₂ = k₃ = g(tₙ+h/2), k₄ = g(tₙ+h). Формула RK4 даёт:
yₙ₊₁ - yₙ = (h/6)·(g(tₙ) + 4·g(tₙ+h/2) + g(tₙ+h))
А это — в точности правило Симпсона для интеграла! (Два средних наклона слились в один с весом 2 + 2 = 4.) Правило Симпсона приближает функцию параболой и интегрирует её точно для многочленов до третьей степени, давая ошибку O(h⁵) на шаге. Веса 1, 2, 2, 1 — это и есть веса Симпсона 1, 4, 1, «расщеплённые» на два серединных замера. Так что RK4 — это обобщение Симпсона на случай, когда подынтегральная функция зависит ещё и от самого решения y.
Условия порядка
Когда f зависит и от y, всё сложнее: чтобы воспровести ряд Тейлора решения вплоть до члена h⁴, коэффициенты a, b, c должны удовлетворять системе из восьми нелинейных алгебраических уравнений — условий порядка для четвёртого порядка. Кутта в 1901 году выписал их и нашёл удобное решение: c = (0, 1/2, 1/2, 1), b = (1/6, 2/6, 2/6, 1/6), а матрица a — с единственными ненулевыми поддиагональными элементами 1/2, 1/2, 1. Это решение не единственное (есть и другие RK4), но оно самое простое и симметричное, поэтому стало «классическим». Веса не подобраны на глаз — они однозначно вытекают из требования совпасть с рядом Тейлора до h⁴ включительно.
Локальная и глобальная ошибки
Раз RK4 точно воспроизводит члены ряда вплоть до h⁴, первая отброшенная величина — порядка h⁵. Это локальная ошибка — ошибка, которую вносит один-единственный шаг: она есть O(h⁵). За весь отрезок мы делаем число шагов, пропорциональное 1/h, и локальные ошибки накапливаются; в итоге глобальная ошибка теряет один порядок и составляет O(h⁴). Практический смысл этого огромен: уменьшив шаг вдвое, мы уменьшаем итоговую ошибку в 2⁴ = 16 раз. Для сравнения, у Эйлера тот же приём даёт всего двукратное улучшение. Поэтому RK4 при умеренном h часто выдаёт точность, которой Эйлеру не достичь и при шаге в тысячу раз меньшем.
Как работает под капотом
Проследим один шаг руками на задаче y' = y, y₀ = 1, h = 0.5, чтобы числа из формулы стали осязаемыми. Поскольку f(t, y) = y, каждый наклон равен текущему значению y в соответствующей точке:
k₁ = y₀ = 1
k₂ = y₀ + (h/2)·k₁ = 1 + 0.25·1 = 1.25
k₃ = y₀ + (h/2)·k₂ = 1 + 0.25·1.25 = 1.3125
k₄ = y₀ + h·k₃ = 1 + 0.5·1.3125 = 1.65625
сумма с весами = k₁ + 2k₂ + 2k₃ + k₄
= 1 + 2.5 + 2.625 + 1.65625 = 7.78125
y₁ = 1 + (0.5/6)·7.78125 = 1 + 0.6484375 = 1.6484375Точное значение e⁰·⁵ ≈ 1.6487213. Ошибка одного шага — около 3·10⁻⁴ при огромном шаге h = 0.5. Эйлер на этом же шаге дал бы y₁ = 1.5 с ошибкой почти 0.15 — в пятьсот раз хуже. Проверим арифметику кодом, который печатает все четыре наклона:
import math
f = lambda t, y: y # y' = y
t, y, h = 0.0, 1.0, 0.5
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_next = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
print(f"k1 = {k1:.6f}")
print(f"k2 = {k2:.6f}")
print(f"k3 = {k3:.6f}")
print(f"k4 = {k4:.6f}")
print(f"y1 = {y_next:.7f}")
print(f"точно = {math.exp(0.5):.7f}")
print(f"ошибка = {abs(math.exp(0.5) - y_next):.2e}")
Вывод:
k1 = 1.000000 k2 = 1.250000 k3 = 1.312500 k4 = 1.656250 y1 = 1.6484375 точно = 1.6487213 ошибка = 2.84e-04
Частые ошибки
- Перепутать аргументы k₃ и k₄. Третий наклон считается в середине (t+h/2) по сдвигу (h/2)·k₂, а четвёртый — в конце (t+h) по полному сдвигу h·k₃. Лёгкая описка типа h/2 вместо h в k₄ молча роняет порядок до второго.
- Поставить веса 1, 1, 1, 1 или 1, 2, 2, 1 без деления на 6. Множитель h/6 обязателен: сумма весов внутри скобок равна 6 (1+2+2+1), и деление на 6 возвращает среднее. Без него шаг вырастет в шесть раз.
- Использовать k₂ там, где нужен k₃. k₃ строится поверх k₂ (уточняет середину). Если оба наклона середины взять одинаковыми (k₃ = k₂), порядок упадёт.
- Ждать локальной ошибки O(h⁴). O(h⁴) — это глобальная ошибка. Локальная (за один шаг) на порядок мельче, O(h⁵); путаница искажает оценки точности.
Итоги
- RK4 берёт четыре наклона: в начале (k₁), дважды в середине (k₂, k₃) и в конце (k₄).
- Итог собирается с весами (1, 2, 2, 1)/6 — середина учитывается вдвое сильнее краёв.
- Веса напрямую связаны с правилом Симпсона: при f, зависящей только от t, RK4 совпадает с формулой Симпсона.
- Веса не подобраны вручную, а следуют из условий порядка — требования совпасть с рядом Тейлора до h⁴.
- Локальная ошибка шага O(h⁵), глобальная — O(h⁴): деление шага пополам улучшает точность в 16 раз.