Методы Рунге—Кутты: от RK2 до RK4 — схема, порядок, выбор шага

В уроке про метод Эйлера мы посчитали $y^{\prime }=x+y$ с шагом $h=0,2$ и накопили ошибку $0,188$. Инженер говорит: «нужно в сто раз точнее» — и у Эйлера только один ответ: в сто раз больше шагов. Платить сотней вычислений за качество — расточительно, тем более что наклон на каждом шаге Эйлер берёт в одной-единственной точке, в самом начале отрезка. Рунге и Кутта в конце XIX века задались вопросом: а если опросить поле направлений в нескольких точках внутри шага и составить наклоны с весами? Ответ превзошёл ожидания: четыре грамотно расставленных вычисления дают ошибку порядка $h^4$ — шаг крупнее в 16 раз за ту же точность, а выигрыш растёт с каждым порядком.

Лестница порядков: от Эйлера до RK4#

Все методы семейства выглядят одинаково: шаг складывается из взвешенной суммы наклонов $k_i$, вычисленных в специально подобранных точках внутри шага. Разница — в количестве опросов и весах. Эйлер: один наклон в начале, веса $1$. Улучшенный Эйлер (метод Хойна): черновой шаг, затем среднее наклонов в начале и в конце — точность $O(h^2)$. Ещё один второй порядок — метод средней точки: наклон в середине отрезка, оценённый по черновому шагу. Далее RK3, и наконец классика:

Прочитаем схему как маршрут: $k_1$ — наклон в начале шага; $k_2$ — наклон в середине, куда добрались полушагом по $k_1$; $k_3$ — снова середина, но уже с уточнённым подъёмом по $k_2$; $k_4$ — наклон в конце, куда ведёт весь шаг по $k_3$. Итог — средневзвешенное с весами $1,2,2,1$, делённое на 6. Каждый $k_i$ использует предыдущий, а не следующий — перепутать порядок значений $k_2$ и $k_3$ местами — типичная ошибка программирования, которая тихо портит точность до второго порядка.

Веса $1,2,2,1$ не случайны: это правило Симпсона из численного интегрирования — средневзвешенное значений функции в начале, двух серединных и конце. RK4 и есть численное интегрирование равенства $y(x_{k+1})=y(x_k)+\int fdx$ по Симпсону, где подынтегральная функция сама зависит от решения. Вся современная теория одностадийных методов упакована в таблицу Бутчера — компактную запись коэффициентов, из которой читается и порядок, и условия согласования; четырёхстадийная схема выше — самый знаменитый член этого семейства.

Один шаг руками: $y' = y$, $y(0) = 1$, $h = 0{,}5$#

Задача скромная, точное решение $y=e^x$ известно: в точке $x=0,5$ оно равно $e^{0,5}\approx 1,6487$. Прогоним один шаг RK4 на бумажке — и сравним с конкурентами.

Порядок точности: что даёт четвёрка#

Глобальная погрешность RK4 имеет порядок $O(h^4)$: при уменьшении шага вдвое ошибка падает в $2^4=16$ раз. Проверим на нашей задаче: один шаг $h=0,5$ дал ошибку $0,0003$; два шага $h=0,25$ дают около $0,000018$ — падение в 16 раз, закон подтверждается. Сравните с Эйлером, где шаг вдвое мельче покупает лишь вдвое меньшую ошибку: за пятнадцать порядков точности в этом споре спрос четырёхкратный, а выигрыш — шестнадцатикратный за шаг.

Откуда берётся четвёртый порядок? Разложение точного решения по Тейлору до $h^4$ содержит производные решения до четвёртой, а схему можно настроить так, чтобы её собственное разложение совпало с тейлоровским до $h^4$ включительно — условия согласования дают систему на коэффициенты метода. Гений Рунге и Кутты — в подборе таких коэффициентов, при которых совпадение достигается всего четырьмя вычислениями $f$: наивный путь требовал бы дифференцировать $f$ по всем переменным, что на практике катастрофа. Симулятор ниже — то же поле направлений, что и в уроке Эйлера, только теперь вы знаете, что между ломаной и кривой прячется ещё и четырёхточечная схема, которая прилипает к истине на глазах.

Когда формулы нет: RK4 как стандарт#

Уравнение $y^{\prime }=\sin (xy)$ не решается в элементарных функциях — и таких уравнений в прикладных задачах большинство. Именно для них и придумано всё численное хозяйство. Прогоните симулятор: поле направлений искривляется, ломаная Эйлера при грубом шаге заметно «сползает» с истинного русла, а там, где решение круто меняется, мелкий шаг становится вопросом выживания.

Для систем уравнений — а системы ДУ и есть норма в приложениях — схема не меняется: все величины становятся векторами, $f$ возвращает вектор, $k_1,\dots ,k_4$ считаются покомпонентно. Именно в векторном виде RK4 зашит в MATLAB (ode45 — на деле адаптивный Рунге—Кутта—Фельдберг), в SciPy (solve_ivp) и в половине симуляторов физики в играх.

Как выбирать шаг: правило Рунге#

Точную ошибку численник не знает — иначе был бы не численник. Но есть честная оценка. Посчитайте решение дважды: с шагом $h$ и с шагом $h/2$ до одной и той же конечной точки. Расхождение ответов связано с ошибкой мелкого прогона:

Логика та же, что у экстраполяции Ричардсона: ошибка метода порядка $p$ — это $C{h}^p$, и из двух прогонов константа $C$ выпаливается делением. Если $R$ много больше требуемой точности — дробите шаг; если много меньше — можно укрупнить и сэкономить время. Профессиональные солверы делают это автоматически на каждом шаге, и техника управления шагом там устроена ещё хитрее (встроенная оценка ошибки из пары методов разного порядка, как в RK45). Тот же критерий сходимости при дроблении работает и в разностных схемах для уравнений в частных производных — там о разности порядков и устойчивости разностной схемы спрашивают так же регулярно.

• RK4: четыре вычисления $f$, веса $1,2,2,1$, делённые на 6; глобальная ошибка $O(h^4)$
• Шаг вдвое мельче — ошибка в 16 раз меньше; у Эйлера только в 2, у Хойна в 4
• $k_i$ считаются последовательно: каждый использует предыдущий, а не черновой
• Шаг выбирает правило Рунге: два прогона, расхождение, делённое на $2^p-1$, — оценка ошибки

Границы применимости тоже стоит знать. Жёсткие задачи — системы с сильно разными скоростями процессов (химическая кинетика, схемы с очень разными постоянными времени) — требуют шага, диктуемого самым быстрым процессом; явные методы на них разоряются, и берут неявные варианты Рунге—Кутты. Вторая тонкость — долгие расчёты механических систем: классический RK4 не сохраняет энергию и фазовый объём, накопительный дрейф лечится симплектическими методами — отдельной ветвью того же семейства, которой пользуются при интегрировании орбит на миллиарды шагов.