📐 Ейлер проти Рунге-Кутти: який метод інтегрування обрати для вашої симуляції?

Задача інтегрування

Більшість фізичних симуляцій описуються звичайними диференціальними рівняннями (ЗДР). Кут маятника θ задовольняє рівняння d²θ/dt² = −(g/L)·sin(θ). Гравітаційна N-тільна система задовольняє другий закон Ньютона для кожного тіла. Рідина описується рівняннями Нав'є-Стокса. Щоб симулювати ці системи, нам потрібно інтегрувати ЗДР у часі.

Основна складність: комп'ютери можуть представляти лише дискретні стани. Ми обираємо крок часу Δt і запитуємо: «маючи стан у момент t, яким буде стан у момент t + Δt?» Відповідь завжди є наближенням — і якість цього наближення залежить від методу інтегрування. Невдалий вибір призводить до дрейфу енергії, нестабільності або відвертого чисельного вибуху.

Метод Ейлера: простий, але небезпечний

Метод Леонарда Ейлера, що сягає 1768 року, є найпростішим підходом. Для ЗДР першого порядку dy/dt = f(y, t) він робить крок уперед так:

y(t + Δt) ≈ y(t) + Δt · f(y(t), t)

Це розклад Тейлора першого порядку — він використовує лише похідну у поточній точці для екстраполяції. Локальна похибка відсікання (похибка на один крок) пропорційна Δt², а глобальна похибка, накопичена за фіксований проміжок часу, становить O(Δt) — звідси й назва «метод першого порядку».

Практична проблема — дрейф енергії. Для консервативних систем, таких як маятники чи орбіти, метод Ейлера систематично додає енергію до системи. Маятник, симульований методом Ейлера, буде повільно розкручуватися назовні, нескінченно нарощуючи амплітуду. Зменшення Δt сповільнює дрейф — але ніколи його не зупиняє. Це фундаментальна властивість методу, а не помилка, яку можна виправити меншим кроком часу.

Рунге-Кутта 4-го порядку: робоча конячка

Метод Рунге-Кутти четвертого порядку (RK4) обчислює похідну в чотирьох точках усередині кожного кроку часу й комбінує їх зі зваженим середнім:

k1 = f(y,       t)
k2 = f(y + k1·Δt/2,  t + Δt/2)
k3 = f(y + k2·Δt/2,  t + Δt/2)
k4 = f(y + k3·Δt,    t + Δt)

y(t+Δt) = y(t) + (Δt/6)·(k1 + 2k2 + 2k3 + k4)

Локальна похибка відсікання становить O(Δt⁵), а глобальна похибка — O(Δt⁴). Подвоєння кроку часу помножує похибку на 16. Зменшення його вдвічі ділить похибку на 16. Ця драматична чутливість до розміру кроку є причиною того, чому RK4 настільки потужний — помірні кроки часу забезпечують високу точність.

RK4 потребує чотирьох обчислень функції на крок (проти одного для Ейлера), тож коштує приблизно у 4 рази дорожче за крок. Але оскільки він досягає точності при більших кроках часу, загальна вартість досягнення заданого рівня точності значно нижча. Для більшості гладких фізичних задач RK4 є типовим вибором за замовчуванням.

Проблема дрейфу енергії: симплектичні інтегратори

RK4 значно кращий за Ейлера, але він все ще має фундаментальну ваду для консервативних механічних систем: він не є симплектичним. Симплектичний інтегратор зберігає геометричну структуру гамільтонової механіки — зокрема, зберігає об'єм фазового простору (теорема Ліувілля). Несимплектичні методи порушують цю властивість і породжують вікові (secular) дрейфи енергії.

Найпростіший симплектичний інтегратор — це метод «жаб'ячих стрибків» (leapfrog, або метод Штермера-Верле). Для системи з позицією x і швидкістю v:

// Оновлення швидкості на половину кроку
v_half = v(t) + (Δt/2) · a(x(t))

// Оновлення позиції на повний крок
x(t+Δt) = x(t) + Δt · v_half

// Обчислення прискорення в новій позиції
a_new = f(x(t+Δt))

// Завершення кроку швидкості
v(t+Δt) = v_half + (Δt/2) · a_new

Leapfrog має точність лише другого порядку — гірше за RK4 на крок. Але оскільки він симплектичний, він зберігає енергію в середньому на довгих часових масштабах. Енергія коливається навколо істинного значення, а не дрейфує монотонно. Для довготривалих симуляцій консервативних систем — планетних орбіт, молекулярної динаміки, прискорювачів частинок — симплектичні інтегратори є необхідними.

Симплектичні методи вищого порядку

Форест і Рут (1990) та Йошида (1990) показали, як побудувати симплектичні інтегратори вищого порядку, комбінуючи кроки leapfrog з ретельно підібраними коефіцієнтами. Симплектичний інтегратор четвертого порядку потребує приблизно 4–7 підкроків leapfrog на повний крок, досягаючи точності рівня RK4 при збереженні симплектичної структури. Ці методи є стандартом у довготривалих кодах молекулярної динаміки та N-тільної астрономії.

Матриця компромісів:

• Ейлер: Швидкий, першого порядку, з дрейфом енергії. Використовуйте лише для швидких прототипів чи неконсервативних систем.
• RK4: Чотири обчислення на крок, четвертого порядку, з незначним дрейфом енергії. Найкращий вибір для точних коротко- та середньострокових симуляцій гладких систем.
• Leapfrog/Verlet: Два обчислення на крок, другого порядку, зберігає енергію в середньому. Найкращий вибір для довготривалих симуляцій консервативних систем.
• Симплектичний 4-го порядку: 5–7 обчислень на крок, четвертого порядку, симплектичний. Поєднує найкраще з обох світів для довгих, вимогливих симуляцій.

Адаптивний розмір кроку

Методи з фіксованим кроком часу мають труднощі, коли розв'язок швидко змінюється в одних областях і повільно — в інших. Адаптивні методи оцінюють локальну похибку, порівнюючи розв'язки різних порядків (наприклад, RK4 проти RK5), і динамічно коригують Δt. Метод Дормана-Принса (використовується у функції ode45 MATLAB та solve_ivp SciPy) досягає цього лише шістьма обчисленнями функції. Розмір кроку збільшується, коли похибка мала (для збільшення швидкості), і зменшується, коли похибка велика (для збереження точності).