Порівняння розв'язувачів ОДР: Ейлер, RK4, Верле, адаптивні методи

1. Задача ОДР та поняття похибки

Звичайне диференціальне рівняння (ОДР) виражає похідну за часом вектора стану y як функцію того ж самого стану:

dy/dt = f(t, y) y(t₀) = y₀ Приклади: Маятник: y = [θ, ω], dy/dt = [ω, -(g/L)sin θ] Орбіта: y = [x, y, vx, vy], dy/dt = [vx, vy, -GMx/r³, -GMy/r³] Логістичне: y = N, dy/dt = rN(1 - N/K) Чисельний інтегратор просуває y від t до t+Δt, використовуючи лише обчислення f — він ніколи не бачить розв'язку у замкненій формі.

Важать два типи похибки:

Локальна похибка апроксимації (LTE): похибка, зроблена за один крок, без урахування попередньої похибки. Метод порядку p має LTE ∝ Δt^(p+1).
Глобальна похибка: накопичена похибка після інтегрування від t₀ до T. Для методу p-го порядку глобальна похибка ∝ Δt^p (на один степінь менша за LTE, бо існує T/Δt кроків).

Метод з глобальною похибкою 2-го порядку зменшує свою похибку вдвічі, коли Δt зменшено вдвічі — а метод 4-го порядку зменшує похибку в 16 разів — тобто методи вищих порядків дуже швидко окуповують свою додаткову обчислювальну вартість, коли потрібна точність.

2. Прямий метод Ейлера: простий, але небезпечний

Найпростіший можливий інтегратор: візьми похідну в поточній точці та зроби крок у цьому напрямку.

// Прямий (явний) Ейлер y(t+Δt) = y(t) + Δt \cdot f(t, y(t)) Вартість: 1 обчислення функції на крок Порядок: 1 (глобальна похибка \propto Δt) Приклад — гармонічний осцилятор: y'' = -y \to y = [x, v] x_new = x + Δt \cdot v v_new = v + Δt \cdot (-x) Точний розв'язок: x(t) = A cos(t), енергія = стала Результат Ейлера: енергія зростає необмежено для будь-якого Δt > 0 (орбіта закручується назовні — інтегратор нестійкий для коливань)

Прямий метод Ейлера безумовно збільшує енергію для консервативних коливань — чисельна орбіта набирає енергію на кожному кроці, хоч би яким малим був Δt. За Δt = 0,01 для одиничного гармонічного осцилятора амплітуда зростає на ~0,5% за період орбіти. Після 1000 періодів амплітуда подвоюється. Це робить його непридатним для будь-якої симуляції, що потребує довгострокового збереження енергії.

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

Область стійкості

Для тестового рівняння dy/dt = λy метод Ейлера стійкий, коли |1 + λΔt| ≤ 1. Для суто уявних λ (коливання) це коло ніколи не містить уявної осі — що підтверджує, що метод Ейлера завжди нестійкий для консервативних осциляторів незалежно від розміру кроку часу.

3. Симплектичний Ейлер і чому він важливий

Зміна в один символ порівняно з прямим методом Ейлера дає кардинально кращий інтегратор для гамільтонових (енергозберігальних) систем: використовуй оновлену швидкість при обчисленні нової позиції.

// Симплектичний (напівнеявний) Ейлер v_new = v + Δt \cdot a(x) // спочатку оновлюємо швидкість x_new = x + Δt \cdot v_new // використовуємо НОВУ швидкість для позиції проти прямого Ейлера: v_new = v + Δt \cdot a(x) // те саме x_new = x + Δt \cdot v // використовує СТАРУ швидкість Вартість: 1 обчислення сили на крок (те саме, що й прямий Ейлер) Порядок: 1 (та сама локальна точність) АЛЕ: точно зберігає модифікований гамільтоніан \to відсутній дрейф енергії!

Ключова відмінність геометрична: симплектичний Ейлер зберігає симплектичну структуру фазового простору — площі областей фазового простору зберігаються під дією відображення. Це дискретна версія теореми Ліувілля. Наслідок у тому, що чисельна орбіта симплектичного інтегратора лежить точно на трохи збуреній енергетичній поверхні, коливаючись довкола справжньої енергії, а не дрейфуючи від неї.

Порівняння енергії за 10 000 орбітальних періодів: Прямий Ейлер (Δt=0.01): енергія дрейфує +47% (орбіта втікає) Симплектичний Ейлер (Δt=0.01): енергія коливається \pm0.01% (обмежена) RK4 (Δt=0.01): енергія коливається \pm10⁻⁷ (дуже малий дрейф) Для фізичного рушія гри, що працює на 60 fps (Δt\approx0.016 с), симплектичний Ейлер не дає об'єктам спонтанно набирати енергію — прямий Ейлер спричинив би помітне «підстрибування» в пружинних системах.

Саме тому практично кожен фізичний рушій гри (Box2D, Bullet, інтеграція PhysX в Unity) використовує симплектичний Ейлер, а не прямий. Він коштує стільки ж, але поводиться принципово краще.

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

RK4 обчислює похідну в чотирьох ретельно вибраних точках всередині кожного кроку часу й комбінує їх у зваженому середньому, що скасовує члени похибки до 4-го порядку включно.

// Класичний Рунге-Кутта 4-го порядку (RK4) k₁ = f(t, y) k₂ = f(t + Δt/2, y + Δt/2 \cdot k₁) k₃ = f(t + Δt/2, y + Δt/2 \cdot k₂) k₄ = f(t + Δt, y + Δt \cdot k₃) y(t+Δt) = y(t) + (Δt/6)(k₁ + 2k₂ + 2k₃ + k₄) Вартість: 4 обчислення функції на крок Порядок: 4 (глобальна похибка \propto Δt⁴) Порівняння похибки з прямим Ейлером за того ж Δt: Ейлер (порядок 1): похибка \propto Δt RK4 (порядок 4): похибка \propto Δt⁴ \to ~10⁴\times менша за Δt=0.1

Коефіцієнти (1/6, 1/3, 1/3, 1/6) відповідають правилу Сімпсона — RK4 по суті застосовує правило Сімпсона до інтегральної форми ОДР. Два обчислення в серединній точці, k₂ та k₃, уточнюють одне одного, скасовуючи член похибки, який залишило б одне серединне обчислення (метод середньої точки 2-го порядку).

Практична точність

RK4 з Δt = 0,1 зазвичай точніший, ніж прямий Ейлер з Δt = 0,001 — за тієї самої загальної кількості обчислень функції. Це і є основний аргумент на користь методів вищих порядків: вони купують точність значно дешевше, ніж просте зменшення кроку часу.

Недоліки RK4

Не симплектичний: для довгочасної гамільтонової динаміки RK4 все одно демонструє повільний дрейф енергії (хоча значно повільніший за прямий Ейлер). Leapfrog з тією самою вартістю перевершує RK4 у збереженні енергії впродовж мільйонів кроків.
Чотири обчислення: якщо f дороге (наприклад, гравітація від N тіл — O(N²) на обчислення), 4-кратна вартість має значення.
Фіксований розмір кроку: класичний RK4 не адаптує Δt. Жорсткі задачі вимагають дуже малого Δt усюди, навіть якщо більша частина розв'язку змінюється повільно.

5. Верле та Leapfrog: довгочасна точність

Інтегратор Верле був розроблений у 1967 році для симуляції молекулярної динаміки — контекстів, де потрібно обчислити трильйони кроків часу, зберігаючи повну енергію. Він використовує оборотність у часі та симплектичну структуру ньютонівської механіки.

// Штермер-Верле (початкова форма) x(t+Δt) = 2\cdotx(t) - x(t-Δt) + Δt² \cdot a(t) Потребує ДВОХ попередніх позицій — незручно стартувати та відстежувати швидкості. // Velocity Verlet (еквівалентна, стандартна форма) x(t+Δt) = x(t) + v(t)\cdotΔt + ½\cdota(t)\cdotΔt² a(t+Δt) = f(x(t+Δt)) / m // переобчислюємо силу в новій позиції v(t+Δt) = v(t) + ½ \cdot (a(t) + a(t+Δt)) \cdot Δt Вартість: 1 обчислення сили на крок (a(t+Δt) повторно використовується на наступному кроці) Порядок: 2 для позицій та швидкостей АЛЕ: симплектичний \to похибка енергії обмежена, не дрейфує

Форма Leapfrog

Leapfrog алгебраїчно еквівалентний Velocity Verlet, але зберігає швидкості зі зсувом на півкроку часу («перестрибуючи» одна одну):

// Leapfrog (швидкості з півкроком) v(t + Δt/2) = v(t - Δt/2) + Δt \cdot a(t) x(t + Δt) = x(t) + Δt \cdot v(t + Δt/2) Позиції та швидкості чергуються: x на цілих кроках, v на півкроках. Точно оборотний у часі: запусти назад \to точно повторює орбіту. Симплектичний: об'єм фазового простору зберігається точно. Похибка енергії коливається довкола 0 — обмежена O(Δt²) — назавжди. Без вікового дрейфу навіть після 10⁹ кроків.

Бенчмарк молекулярної динаміки: інтегрування системи з 1000 частинок Леннарда-Джонса протягом 10⁶ кроків: середньоквадратична похибка енергії Leapfrog = 0,003%. RK4 (та сама вартість) = 0,0002%, але зростає — після 10⁸ кроків похибка RK4 перевищує Leapfrog. Для астрономії (орбітальна динаміка N тіл) leapfrog є стандартом з тієї ж причини.

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

Методи Фореста-Рута (4-го порядку) та Йосіди (6-го, 8-го порядку) складають кілька кроків leapfrog з певними додатними та від'ємними підкроками, щоб скасувати похибку вищих порядків, зберігаючи симплектичність. Вони досягають точності RK4/RK8 з довгочасними енергетичними обмеженнями leapfrog — ціною більшої кількості підкроків на повний крок (6–15 обчислень).

6. Адаптивне керування кроком

Інтегратори з фіксованим кроком марнують обчислення на гладких ділянках розв'язку й можуть пропустити швидкі зміни. Адаптивні методи безперервно оцінюють локальну похибку та коригують Δt, щоб утримувати заданий користувачем допуск:

// Вкладена пара Рунге-Кутта (наприклад, Дорманд-Принс RK45) Два розв'язки обчислюються з ТИХ САМИХ k обчислень: y4 — розв'язок 4-го порядку y5 — розв'язок 5-го порядку (1 додаткове обчислення) Оцінка похибки: err = |y5 - y4| (різниця між порядками) Контролер розміру кроку (ПІ-контролер): Якщо err > tol: відхили крок, зменш Δt: Δt_new = Δt \cdot (tol/err)^(1/5) Якщо err < tol: прийми крок, можливо збільш Δt: Δt_new = Δt \cdot (tol/err)^(1/5) \cdot safety

Пара Дорманда-Принса RK45 (основа ode45 у MATLAB та RK45 у SciPy) використовує 6 обчислень на крок з точністю 5-го порядку. Оцінка 4-го порядку дістається «безкоштовно» з тих самих значень k — це забезпечується хитрим вибором коефіцієнтів таблиці Бутчера.

Кеш-Карп та Богацький-Шемпайн

Богацький-Шемпайн (RK23) — дешевша пара 3/2-го порядку, що використовує лише 4 обчислення — ідеальна для застосунків реального часу або коли розв'язок має обмежену гладкість. Кеш-Карп (RK45) був попереднім стандартом MATLAB; Дорманд-Принс трохи точніший за тієї ж вартості й став сучасним стандартом за замовчуванням.

Налаштування допуску: адаптивні розв'язувачі приймають відносний допуск rtol та абсолютний допуск atol. Використовуйте rtol = 1e-6, atol = 1e-9 для точної наукової роботи. rtol = 1e-3 для інтерактивного використання в реальному часі. Дуже жорсткі допуски (< 1e-12) зазвичай не мають сенсу, бо домінує похибка округлення з рухомою комою.

7. Жорсткі ОДР та неявні методи

Жорстка ОДР має компоненти з кардинально різними часовими масштабами. Класичний приклад: мережа хімічних реакцій, де одні речовини реагують за мікросекунди, тоді як інші змінюються годинами. Явний розв'язувач мусить використовувати найменший релевантний Δt (мікросекунди) впродовж усієї симуляції (години) — 10¹⁰ кроків. Неявний розв'язувач може обробити швидку компоненту аналітично й робити великі кроки для повільної динаміки.

// Неявний (зворотний) Ейлер — найпростіший неявний метод y(t+Δt) = y(t) + Δt \cdot f(t+Δt, y(t+Δt)) // використовує МАЙБУТНІЙ стан Це нелінійне рівняння відносно y(t+Δt) — потрібно розв'язувати ітеративно (метод Ньютона на кроці). Кожна ітерація потребує якобіана \partialf/\partialy. Стійкість: A-стійкий — стійкий для ВСІХ λΔt з Re(λ) < 0. Явний RK4 стійкий лише для |λΔt| в обмеженій області. Вартість на крок: значно вища за явний (ітерації Ньютона) Але: Δt може бути в 10³-10⁶\times більшим на жорстких задачах \to чистий ВИГРАШ

Методи BDF та Розенброка

Методи формули диференціювання назад (BDF) Гіра — що використовуються в ode15s MATLAB — є галузевим стандартом для жорсткої хімічної кінетики, симуляції електричних кіл та структурної динаміки з широко змінними жорсткостями пружин. LSODA (використовується в SciPy) автоматично перемикається між нежорсткими (Адамса-Башфорта) та жорсткими (BDF) інтеграторами на основі виявленої жорсткості.

Методи Розенброка лінеаризують неявне рівняння аналітично, замість того щоб використовувати ітерацію Ньютона — дешевші на крок, але обмежені помірними коефіцієнтами жорсткості.

Чи жорстка ваша задача? Обчисліть власні числа якобіана λᵢ. Якщо max|Re(λᵢ)| / min|Re(λᵢ)| перевищує ~10³, вона жорстка. Практичні ознаки: явні розв'язувачі вимагають крихітного Δt, щоб лишатися стійкими, навіть коли розв'язок видається гладким; інтегратор постійно відхиляє та зменшує кроки на тому, що схоже на нудну ділянку траєкторії.

8. Який розв'язувач варто використати?

Сценарій застосування	Рекомендований розв'язувач	Причина
Фізика гри (тверді тіла, пружини)	Симплектичний Ейлер	1 обчислення, без дрейфу енергії, дешевий, ~60 fps
Орбітальна механіка, N тіл	Leapfrog / Velocity Verlet	Симплектичний, 1 обчислення, стійкий упродовж мільйонів кроків
Молекулярна динаміка	Velocity Verlet	Галузевий стандарт, оборотний у часі, ансамбль NVE
Гладка траєкторія, помірна точність	RK4	4-й порядок, простий у реалізації, широко відомий
Висока точність, невідома гладкість	Дорманд-Принс RK45	Адаптивний, автоматичне керування похибкою, `ode45`
Жорстка хімія / електричні кола	BDF (ode15s / LSODA)	A-стійкий, опрацьовує співвідношення масштабів 10⁶:1
Довгочасний гамільтонів (високий порядок)	Йосіда 6-го порядку	Симплектичний + висока точність, без вікового дрейфу
Швидкий прототип / навчання	Прямий Ейлер	5 рядків коду, прямо візуалізує похибки

Підсумок: три ключові властивості

Порядок (точність на крок): вищий порядок → менше потрібних кроків → менше загальної роботи. RK4 (порядок 4) перевершує Ейлера (порядок 1) у 10⁴× за Δt = 0,1.
Симплектичність (довгочасне збереження енергії): Leapfrog та симплектичний Ейлер зберігають модифікований гамільтоніан. Несимплектичні методи (зокрема RK4) дрейфують на збережних величинах упродовж дуже довгих інтегрувань.
A-стійкість (опрацювання жорсткості): неявні методи є A-стійкими — вони не мають обмеження розміру кроку через стійкість. Явні методи вимагають Δt < O(1/|λ_max|) незалежно від того, наскільки повільно насправді розвивається цікава динаміка.

Жоден окремий інтегратор не виграє за всіма трьома. Найкращий вибір залежить від фізики (консервативна проти дисипативної проти жорсткої), потрібної точності та вартості одного обчислення сили/похідної.