Методы Эйлера и Рунге — Кутты

Где были Солнце и планеты 17 октября 1926 года и с какими скоростями они двигались, сегодня известно очень точно: это хранят эфемериды JPL. Закон движения тоже известен — тяготение по Ньютону. Но формулы, которая выдала бы положения девяти тел через сто лет, нет: закон говорит только, как быстро меняются положения и скорости прямо сейчас. Остаётся шагать — сдвинуть всё на маленький шаг, пересчитать силы, снова сдвинуть. Ниже — сто лет шагов от 1926 года до сегодняшнего дня, и каждый год расчёт сверяется с тем, где планеты были на самом деле.

Метод Рунге — Кутты четвёртого порядка с шагом в сутки ставит Землю 2026 года в 3 400 км от её настоящего места — меньше земного радиуса, Юпитер — в 474 км. А Меркурий уезжает на 48 млн км: сутки для него слишком длинный шаг, он делает оборот за 88 дней. Переключитесь на Эйлера с тем же шагом: Меркурий теряется в первый же год, Земля — во второй, а к 2026-му энергия системы уходит на 13 %. Метод Верле не теряет никого и держит энергию с точностью до десятимиллионных, но Земля у него в 9 млн км от своего места: орбита почти та же, а сама планета сдвинулась по ней на три с половиной градуса. У каждого метода своя ошибка, и об этом статья.

Идея

Состояние системы — набор чисел y: положения и скорости. Закон движения — правило y′ = f(t, y): по состоянию оно говорит, куда и как быстро оно меняется. Метод Эйлера делает шаг h по этой касательной: y ← y + h·f(t, y). Но наклон за время шага меняется, и касательная уводит в сторону. Методы Рунге — Кутты вычисляют наклон ещё в нескольких точках внутри шага и складывают их с весами, подобранными так, чтобы ошибки взаимно гасились. Классический РК4 берёт четыре наклона, и его ошибка за всё время расчёта падает как h⁴ — в шестнадцать раз при шаге вдвое меньше.

Что такое дифференциальное уравнение, поле направлений и почему у пружины и у «хищника — жертвы» нет общих формул, рассказывает глава о дифференциальных уравнениях «Царицы наук»; там же первый шаг Эйлера на числе e. Здесь — алгоритмы: сколько стоит точность, почему точный метод может проиграть неточному на долгом расчёте и когда любая явная схема взрывается.

Шаг Эйлера и цена ошибки

За один шаг Эйлер ошибается на величину порядка h²: касательная расходится с кривой квадратично. Шагов до момента T — T/h, и кажется, что общая ошибка ~ h. Так и есть, но с оговоркой, которую обычно пропускают: ошибки не просто складываются, а ещё и размножаются, потому что каждый следующий шаг стартует уже из неверной точки.

Множитель h — это первый порядок: шаг вдвое меньше, ошибка вдвое меньше. Множитель e^Lt — это расплата за длину расчёта: ошибки, сделанные в начале, раздуваются вместе с расхождением соседних траекторий. Для орбит оценка слишком пессимистична, но направление верное: ошибка Эйлера растёт быстрее, чем время расчёта, а уменьшить её в сто раз стоит в сто раз больше работы. В первой сцене Эйлер с шагом в четверть суток делает вчетверо больше работы и всё равно теряет Меркурий и Землю в первые три года.

У Эйлера на орбитах есть и вторая беда — систематическая. Для колебаний x″ = −ω²x шаг Эйлера переводит (x, v) в (x + hv, v − hω²x), и прямой подстановкой проверяется, что энергия (v² + ω²x²)/2 умножается на 1 + h²ω² — на каждом шаге, без исключений. Метод постоянно подкачивает в систему энергию. С гравитацией выходит то же: в первой сцене Эйлер раскручивает внутренние планеты наружу, и Меркурий уходит с орбиты в первый же год.

Четыре наклона

Идея Рунге: если наклон на шаге меняется, возьмём его не только в начале. Метод средней точки сначала делает полшага по касательной, смотрит наклон там и уже с ним делает весь шаг — ошибка ~ h². РК4 идёт дальше: наклон k₁ в начале, k₂ в середине по k₁, k₃ снова в середине, но по k₂, k₄ в конце по k₃, и шаг

Если f не зависит от y, это ровно формула Симпсона из статьи о численном интегрировании: веса 1, 4, 1 при двух средних наклонах, разделённых пополам. Коэффициенты любого метода Рунге — Кутты удобно записывать таблицей Бутчера: в строке — где взять очередной наклон и с какими весами сложить предыдущие, внизу — веса ответа.

Почему именно так: разложим точное решение и ответ метода в ряд Тейлора по h и потребуем, чтобы совпали члены до h⁴. Получаются восемь уравнений на коэффициенты, и у них много решений; классическое — самое симметричное. Проверка громоздкая, и честнее не повторять её, а измерить порядок. Ниже Земля обходит Солнце ровно год четырьмя методами с разным шагом, а точный ответ даёт уравнение Кеплера, решённое методом Ньютона, — как в статье о бисекции и методе Ньютона.

steps   h, days        Euler     midpoint          RK4     leapfrog
  128     2.854    2.971e+08    1.366e+06        163.5    7.589e+05
  256     1.427     1.82e+08    3.344e+05        9.004    1.899e+05
  512     0.713    1.021e+08    8.258e+04       0.5249    4.747e+04
 1024     0.357    5.418e+07    2.051e+04      0.03162    1.187e+04
order: Euler 0.91, midpoint 2.01, RK4 4.05, leapfrog 2.00

Ошибка РК4 при шаге вдвое меньше падает в 16 раз, средней точки и Верле — в 4, Эйлера — едва вдвое. При 1024 шагах в год РК4 ошибается на 32 метра, а Эйлер — на треть астрономической единицы. Шаг РК4 вчетверо дороже шага Эйлера, но сравнение при равной работе ещё беспощаднее: 256 шагов РК4 стоят столько же, сколько 1024 шага Эйлера, и ошибаются на 9 км против 54 млн.

Шаг по месту: пара Дормана — Принса

Одного шага на всю дорогу мало. Комета Галлея полжизни ползёт за орбитой Нептуна, где за сутки почти ничего не меняется, а у Солнца разворачивается за считанные недели. Постоянный шаг, достаточно мелкий для перигелия, на остальных 75 годах орбиты — пустая трата сил. Решение то же, что у адаптивного интегрирования: оценивать ошибку на каждом шаге и подстраивать шаг под допуск.

Оценить ошибку помогает вложенная пара. Метод Дормана и Принса 1980 года считает семь наклонов и складывает их двумя наборами весов: получаются ответ пятого порядка и ответ четвёртого. Их разница почти целиком — ошибка четвёртого, и это оценка ошибки шага, за которую не нужно ни одного лишнего вычисления. Если она больше допуска — шаг отвергают и повторяют короче; если меньше — принимают, а следующий шаг удлиняют. Ошибка пятого порядка растёт как h⁵, поэтому новый шаг — старый, умноженный на (допуск/ошибка)^1/5, с запасом 0,9. Седьмой наклон берётся в конце шага и служит первым наклоном следующего, так что принятый шаг стоит шесть вычислений.

tolerance 1e-09: 169 steps (+6 rejected), off by 35.7 km; a step of 798 days far out, 1.61 days at the perihelion; 65 steps within 2 au
  RK4 with even steps and as many slopes (262 steps): off by 3.4e+08 km

С допуском 10⁻⁹ метод проходит шестьдесят лет орбиты за 169 шагов: вдали от Солнца шаг больше двух лет, в перигелии 1986 года — полтора дня. 65 шагов из 169 приходятся на ту часть орбиты, что ближе 2 а. е. к Солнцу, хотя комета проводит там около семи месяцев из 76-летнего оборота. Ошибка против точного решения задачи двух тел — 36 км. РК4 с постоянным шагом и тем же числом вычислений — 262 шага по 84 дня — проскакивает перигелий и промахивается на 340 млн км. В первой сцене режим «Дорман — Принс» рисует график шага: он дышит с периодом Меркурия, самого быстрого тела системы, и укорачивается в каждом его перигелии.

Пара Дормана — Принса — метод по умолчанию в scipy.integrate.solve_ivp (под именем RK45) и в ode45 из MATLAB. Это ответ на вопрос «каким методом решать ОДУ, если не знаю, каким»: первым пробуют его.

Верле: не точнее, а честнее с энергией

У движения под действием силы есть особая структура: скорость меняется от положения, положение — от скорости. Метод Верле, он же leap-frog, пользуется этим напрямую. Полшага меняется только скорость — силами в старом положении. Затем целый шаг меняется только положение — уже новой скоростью. Затем ещё полшага скорости — силами в новом положении. Силу из конца шага можно взять в начало следующего, поэтому вычисление сил одно на шаг, как у Эйлера, а порядок второй.

По таблице с Землёй Верле того же порядка, что метод средней точки. Его сила в другом. Запустите сцену ниже: четыре планеты-гиганта с Солнцем на сто тысяч лет вперёд, и у трёх методов одинаковая работа — РК4 делает шаг в 40 дней, Эйлер и Верле — в 10.

Эйлер теряет десятую часть энергии за 10 лет — планеты уходят с орбит. У РК4 ошибка энергии растёт прямо пропорционально времени: за сто лет 8·10⁻⁷, за тысячу 9·10⁻⁶, за сто тысяч 9·10⁻⁴. У Верле она с первых лет колеблется около нескольких миллионных и никуда не уходит. Через 825 лет РК4 перестаёт выигрывать, а к концу проигрывает больше чем в сто раз. Поэтому орбиты на миллионы и миллиарды лет считают методами семейства Верле: в 1991 году Джек Уиздом и Мэтью Холман построили на этой идее схему для многих тел, и её потомками считают эволюцию планетных систем на миллиарды лет вперёд.

Для многих тел то же верно для объёма в пространстве всех положений и скоростей, и даже больше: Верле сохраняет так называемую симплектическую структуру — то, что сохраняет и сама механика. Отсюда и поведение энергии. Можно доказать, что Верле точно решает чуть-чуть другую механическую задачу, с энергией, которая отличается от настоящей на величину порядка h², — и эту «почти энергию» сохраняет очень долго; доказательство длинное, его можно найти в обзоре Хайрера, Любиха и Ваннера «Геометрическое численное интегрирование на примере метода Штёрмера — Верле» (2003). Отсюда колебания вместо дрейфа. Но и плата видна в первой сцене: фаза. Орбита у Верле остаётся правильной формы, а планета на ней спешит или отстаёт — Земля за сто лет сместилась на три с половиной градуса.

Жёсткие уравнения: когда шаг держит не точность

Возьмём задачу попроще планет: кирпичная стена 40 см в самую холодную неделю зимы 2024/25 в Алматы. Внутри 20 °C, снаружи — настоящая температура по часам, от +2,3 до −21,6 °C. Разрежем стену на слои по сантиметру: температура каждого слоя тянется к среднему соседей, и получается система из 39 обычных дифференциальных уравнений. Решения у неё плавные — температура в толще стены меняется за часы и сутки, — так что шаг в час выглядит скромным. Попробуйте.

Явный Эйлер с шагом 100 секунд честно проходит неделю. Со 115 секундами расчёт разносит через сутки: температуры слоёв начинают прыгать то вверх, то вниз через один и растут, пока не уйдут в бесконечность. РК4 терпит шаг чуть длиннее, до 158 секунд, а дальше взрывается так же. А неявный метод Эйлера шагает по часу и ошибается в тепловом потоке на 0,04 %. Это и называется жёсткостью: решение медленное, а шаг явного метода ограничен чем-то другим.

Чем — видно на простейшем уравнении y′ = λy с λ < 0, решение которого e^λt гаснет. Шаг Эйлера умножает y на 1 + hλ, и если hλ < −2, этот множитель меньше −1: вместо затухания — колебания с растущим размахом. В стене такие «λ» есть у каждого способа изогнуть профиль температуры, и самый быстрый — изгиб через слой, зигзаг. Он гаснет меньше чем за минуту, в решении его давно нет, но шаг явного метода обязан с ним считаться.

Для нашей стены Δx = 1 см, α = 4,4·10⁻⁷ м²/с, и предел — 113 секунд. Мельче сетка — жёстче система: вдвое тоньше слои — вчетверо короче допустимый шаг. У РК4 множитель — многочлен 1 + z + z²/2 + z³/6 + z⁴/24 от z = hλ, и он не больше единицы по модулю при −2,785 ≤ z ≤ 0: предел в 1,39 раза дальше, но тоже есть. Высокий порядок от жёсткости не спасает.

Неявный шаг стоит дороже: новые температуры входят в уравнение с обеих сторон, и на каждом шаге приходится решать систему — здесь трёхдиагональную, двумя проходами по слоям, вперёд и назад. Это частный случай метода Гаусса, и для такой матрицы он линейный по числу слоёв. Неделя — 168 таких шагов вместо 6048 явных.

brick: alpha = 4.419e-07 m2/s; explicit Euler is stable only for h < dx^2/(2 alpha) = 113 s
reference, backward Euler with h = 10 s: 49.839 W/m2, 8.373 kWh per m2 a week
explicit h =    100 s:  6048 steps, 49.840 W/m2, off by 0.00%
explicit h =    115 s: blew up after 849 steps (27.1 hours of the week)
backward h =   3600 s:   168 steps, 49.820 W/m2, off by 0.04%

Через квадратный метр такой стены за морозную неделю уходит 8,4 киловатт-часа: модель простая, без сопротивления воздуха у поверхностей, но порядок верный. Жёсткие системы встречаются везде, где рядом живут быстрые и медленные процессы: химические реакции, у которых скорости различаются в миллионы раз, электрические схемы с крошечными паразитными ёмкостями, тонкие сетки в расчётах тепла и прочности. Термин «жёсткие уравнения» ввели Чарлз Кёртисс и Джозеф Хиршфельдер в 1952 году, а условие, которое делает метод пригодным для любой жёсткости, — A-устойчивость — сформулировал Гермунд Дальквист в 1963-м.

Код: гиганты с 1926 года

Первая сцена считает в браузере все тела сразу. Программа ниже делает то, на что хватит десяти секунд песочницы: Солнце вместе с внутренними планетами считает одним телом в их общем центре масс и ведёт четыре планеты-гиганта сто лет — тремя методами при одинаковом числе вычислений сил. В конце она сверяет положения с JPL. Данные те же: /data/algorithms/solar-1926.json.

method     step, d  forces      Jupiter       Saturn       Uranus      Neptune   energy drift
Euler           10    3700    2.133e+09    3.152e+09    4.795e+08    8.291e+07       -3.9e-01
leapfrog        10    3700    3.123e+06    3.436e+05    2.241e+04         7495       -5.0e-06
RK4             40    4000    3.075e+04         1774        293.9          255        9.0e-07
RK4             10   14800    1.304e+04         2478        293.1        255.4        9.3e-10

Последняя строка — РК4 с шагом вчетверо мельче: энергия сохранена в тысячу раз лучше, а Юпитер всё равно в 13 тысячах километров от настоящего. Это уже не ошибка метода, а ошибка модели: мы сложили внутренние планеты в одну точку с Солнцем. В первой сцене, где все планеты считаются отдельно, этот предел для Юпитера — около 470 км, для Земли — около 5500, для Меркурия — около 8 тысяч: дальше мешает то, чего нет в законе Ньютона, — поправки теории относительности, астероиды, а для кометы Галлея ещё и реактивные струи газа у Солнца. Сколько ни уменьшай шаг, ниже ошибки модели не опустишься, и разумный шаг — тот, при котором ошибка метода заметно меньше её.

Поэтому Земля в первой сцене при шаге в сутки и оказалась ближе этого предела, в 3 400 км: это совпадение, а не точность. Обе ошибки направлены вдоль орбиты. Модельная Земля к 2026 году отстаёт от настоящей на 5 500 км, а РК4 с шагом в сутки уводит её вперёд примерно на 2 тысячи, и ошибки отчасти гасят друг друга. С шагом в полсуток ошибка метода почти исчезает, и расхождение вырастает до 5 400 км, а с шагом в двое суток Земля оказывается в 58 тысячах км от места. Поэтому шаг проверяют не совпадением с наблюдениями, а сходимостью: при меньшем шаге ответ не должен заметно меняться.

Сколько это стоит

Работа метода — число вычислений правой части f; у гравитации N тел каждое стоит O(N²). Шагов T/h, а h выбирают из двух ограничений: точности и устойчивости.

Чтобы ошибка уменьшилась в 10 000 раз, Эйлеру нужно в 10 000 раз больше шагов, РК4 — в 10. Поэтому метод первого порядка годится там, где точность не нужна, а нужна устойчивость и простота, как в играх, — и там выбирают полунеявного Эйлера, который сначала обновляет скорость, а положение двигает уже новой скоростью: он тоже сохраняет площадь, как Верле. Для точного расчёта на умеренное время — адаптивный РК высокого порядка. Для долгого — Верле и его родственники. Для жёсткой системы — неявный метод, и тогда главная стоимость шага — решение системы уравнений.

Где это работает

• Научные библиотеки. scipy.integrate.solve_ivp по умолчанию берёт RK45 — пару Дормана — Принса, для жёстких задач предлагает неявные Radau, BDF и LSODA. В MATLAB та же пара — ode45.
• Молекулярная динамика. Белки и жидкости считают на миллиарды шагов, и энергия не должна уплывать: GROMACS по умолчанию шагает методом leap-frog, LAMMPS — скоростным Верле.
• Игры и анимация. Физические движки обновляют сначала скорость, потом положение по новой скорости — полунеявный Эйлер. Так советует, например, автор движка Box2D Эрин Катто.
• Небесная механика. Долгие расчёты планетных систем ведут симплектическими методами; эфемериды для навигации космических аппаратов — методами высокого порядка с переменным шагом.
• Химия, электроника, тепло. Кинетика реакций, симуляторы схем, расчёты теплопроводности на мелкой сетке — жёсткие системы, и там работают неявные методы.

История

Шаг по касательной описал Леонард Эйлер в «Интегральном исчислении» 1768 года. Сто с лишним лет спустя Карл Рунге в 1895 году предложил брать наклон и внутри шага, Карл Хойн в 1900-м построил методы второго и третьего порядка, а Мартин Кутта в 1901 году — тот самый метод четвёртого порядка с весами 1, 2, 2, 1. Таблицы коэффициентов и общую теорию порядков в 1960-х развил Джон Бутчер, его именем таблицы и называют. Пара с оценкой ошибки, ставшая стандартом, — работа Джона Дормана и Питера Принса 1980 года.

У метода Верле история длиннее его имени. Исаак Ньютон в «Началах» 1687 года, доказывая второй закон Кеплера, рассуждал ровно этими шагами: тело летит по прямой, получает толчок к Солнцу, снова летит по прямой, — Ричард Фейнман разбирал этот чертёж в лекциях 1964 года. В 1907 году Карл Штёрмер считал так траектории заряженных частиц в магнитном поле Земли, изучая полярные сияния, а в 1967-м Лу Верле — движение атомов аргона в первых компьютерных опытах с жидкостью. Именно в молекулярной динамике метод и получил своё нынешнее имя.

Задачи

В первой задаче нужно довести тело по орбите до заданного момента — настоящие состояния Земли, кометы и МКС, и одна орбита сильно вытянута. Во второй — та же стена, но с шагом, на котором явный метод взрывается: придётся сделать неявный шаг и решить трёхдиагональную систему.