Практично кожна фізична симуляція зводиться до інтегрування системи звичайних диференціальних рівнянь. Вибір чисельного методу визначає точність, стійкість та обчислювальну вартість — і хибний вибір для жорсткої системи (системи, де одні величини змінюються набагато швидше за інші, через що доводиться брати дуже малий крок, щоб не втратити стійкість) може означати результати, що відрізняються на порядки величини або відверто розбігаються. Ця стаття оглядає ключові методи: від підручникового методу Ейлера через класичний RK4, до адаптивного методу Дорманда–Принса, багатокрокового методу Адамса–Башфорта та симплектичних інтеграторів, що зберігають геометрію.
1. Задача Коші (задача з початковими умовами)
Ми прагнемо розв'язати:
де y ∈ ℝⁿ та f: ℝ × ℝⁿ → ℝⁿ. Усі класичні фізичні системи (другий закон Ньютона, рівняння Лоренца, хімічна кінетика) відповідають цій формі після зведення рівнянь вищих порядків до систем першого порядку.
Чисельний інтегратор апроксимує y(t) у послідовності дискретних часових точок tₖ = t₀ + k·h (фіксований крок) або в адаптивно вибраних tₖ. Основний виклик — балансування трьох конкурентних вимог: точності (локальна похибка апроксимації на крок), стійкості (похибка не зростає необмежено) та вартості (кількість обчислень функції на крок).
2. Метод Ейлера (порядок 1)
Локальна похибка апроксимації (LTE): O(h²). Глобальна похибка: O(h). Найпростіший метод, який майже ніколи не використовується на практиці — він вимагає надзвичайно малого h для прийнятної точності. Явний метод Ейлера також є A-нестійким: для тестового рівняння dy/dt = λy з Re(λ) < 0 метод вимагає h < 2/|λ|, щоб уникнути зростаючих коливань, що дуже обмежувально для жорстких задач.
Симплектичний метод Ейлера — просте покращення
Для гамільтонових систем (q — координата, p — імпульс)
симплектичний метод Ейлера спершу оновлює p, а потім використовує нове p
для оновлення q:
pₙ₊₁ = pₙ + h·f(qₙ)
qₙ₊₁ = qₙ + h·pₙ₊₁/m
Попри лише перший порядок точності, він
зберігає симплектичну структуру (площу у фазовому
просторі), тож похибка енергії лишається обмеженою впродовж експоненційно
довгих часів — на відміну від методів Рунге-Кутта, що повільно накопичують дрейф енергії.
3. Класичний метод Рунге–Кутта (RK4, порядок 4)
«Золотий стандарт» — робоча конячка фізичної симуляції:
LTE: O(h⁵). Глобальна похибка: O(h⁴). Чотири обчислення функції на крок. Константа похибки дуже мала: для гармонічного осцилятора похибка на одиницю часу становить приблизно 1/180 · h⁴ · max|y⁽⁵⁾|. Для більшості нежорстких задач з гладкими розв'язками RK4 з h ≈ 0,01–0,1 є дуже точним. Область стійкості у комплексній площині h·λ — це дископодібна область, що покриває більшу частину лівої півплощини аж до |hλ| ≈ 2,8.
// RK4 мовою JavaScript
function rk4(f, t, y, h) {
const k1 = f(t, y);
const k2 = f(t + h/2, y.map((v, i) => v + h/2 * k1[i]));
const k3 = f(t + h/2, y.map((v, i) => v + h/2 * k2[i]));
const k4 = f(t + h, y.map((v, i) => v + h * k3[i]));
return y.map((v, i) => v + h/6 * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]));
}4. Адаптивне керування кроком: Дорманд–Принс RK45
RK4 з фіксованим кроком марнує обчислення функції в гладких областях і робить надто грубі кроки у швидкозмінних. Адаптивні методи керують розміром кроку, щоб утримувати оцінену локальну похибку в межах допуску (rtol, atol).
Метод Дорманда–Принса (ode45 у Matlab,
RK45 у SciPy) використовує пару вкладених формул RK зі спільними
тими самими 6 стадіями (+ одне повторне використання FSAL «first same as last»): одну 5-го порядку
й одну 4-го порядку. Їхня різниця оцінює локальну похибку без
додаткових витрат.
Прийом FSAL повторно використовує k₇ з кроку n як k₁ кроку n+1, знижуючи вартість з 6 до 5 обчислень на прийнятий крок. Більшість практичних симуляцій використовують rtol = 1e-6, atol = 1e-9 як значення за замовчуванням для наукової точності.
| Метод | Порядок | Обчислень/крок | Адаптивний? | Найкраще застосування |
|---|---|---|---|---|
| Ейлер | 1 | 1 | Ні | Лише для навчання |
| Гойн (RK2) | 2 | 2 | Ні | Дуже грубі апроксимації |
| RK4 (класичний) | 4 | 4 | Ні | Нежорсткі, відомі гладкі ОДР |
| Дорманд–Принс RK45 | 4(5) | 6 (5 з FSAL) | Так | Універсальний за замовчуванням для нежорстких |
| Кеш–Карп RK45 | 4(5) | 6 | Так | Гладкі нежорсткі з жорстким допуском |
| DOP853 | 8(5,3) | 13 | Так | Високоточні нежорсткі |
| Верле / leapfrog | 2 | 1 | Ні | Гамільтонові, N-тіл, МД |
| Штермер–Верле | 2 | 1 | Ні | Симплектичні, збереження енергії |
5. Багатокрокові методи: Адамса–Башфорта та Адамса–Мултона
Однокрокові методи (Рунге–Кутта) відкидають усі раніше обчислені значення функції після кожного кроку. Багатокрокові методи повторно використовують останні k обчислень функції, досягаючи високого порядку лише з одним чи двома новими обчисленнями на крок.
Явна p-крокова формула Адамса–Башфорта (AB-p):
AB-4 (4-кроковий) має порядок 4 лише з 1 новим обчисленням на
крок (проти 4 для RK4). Неявний Адамса–Мултона (AM-p)
досягає порядку p+1 з p кроками. Поєднання явного предиктора +
неявного коректора — це схема PECE (Predict–Evaluate–Correct–Evaluate),
що використовується в ode113 Matlab.
Обмеження: багатокрокові методи потребують стартової процедури (використання RK4 для перших k кроків) і не можуть легко змінювати розмір кроку (потребують перезапуску або інтерполяції).
6. Жорсткі системи та неявні методи
Система є жорсткою, коли найбільше власне число якобіана ∂f/∂y за модулем значно перевищує обернену величину часового масштабу, що нас цікавить: коефіцієнт жорсткості S = max|Re(λᵢ)| · T ≫ 1.
Приклади: хімічна кінетика зі швидкими/повільними реакціями (S може сягати 10¹⁵), RC-коло з широко рознесеними сталими часу, рівняння дифузії після просторової дискретизації (власні числа ∝ −1/h²).
A-стійкість та L-стійкість
Інтегратор є A-стійким, якщо для всіх λ з Re(λ)≤0
виконується |R(hλ)| ≤ 1 для всіх h>0 — коефіцієнт підсилення ніколи не перевищує 1.
Явні методи RK не є A-стійкими (їхні області стійкості
в лівій півплощині обмежені). Неявні методи можуть бути
A-стійкими.
L-стійкість (сильніша): додатково R(∞)=0, тож жорсткі
компоненти згасають миттєво. Неявний метод Ейлера (порядок 1), метод трапецій
(порядок 2, A-стійкий, але не L-стійкий), методи SDIRK та методи
Розенброка забезпечують A/L-стійке інтегрування.
Для жорстких систем неявна схема неявного методу Ейлера
вимагає розв'язання нелінійної системи на кожному кроці (ітерація Ньютона), але
дозволяє довільно великий h без нестійкості. Методи
BDF-k (формула диференціювання назад, методи Гіра),
що використовуються в ode15s MATLAB, досягають порядку до 6,
залишаючись жорстко стійкими.
7. Симплектичні інтегратори для гамільтонових систем
Для консервативних механічних систем з гамільтоніаном H = T(p) + V(q) методи Рунге–Кутта вносять повільний віковий дрейф енергії навіть за високого порядку. Симплектичний інтегратор точно зберігає симплектичну 2-форму ω = Σ dqᵢ ∧ dpᵢ, що на практиці означає:
- Похибка енергії обмежена (коливальна), а не дрейфуюча
- Модифікований гамільтоніан H̃ = H + O(hᵖ) зберігається точно
- Ідеальні для гравітації N тіл, молекулярної динаміки та довгострокової орбітальної механіки
Інтегратор Верле/leapfrog є симплектичним 2-го порядку:
Симплектичні методи вищих порядків (4-го порядку Йосіди, 4-го порядку Рута–Фореста) застосовують відображення Верле з хитро масштабованими підкроками. Ціна: відсутність адаптивного кроку — крок мусить лишатися фіксованим, щоб симплектичність зберігалася.
8. Інтерактив: порівняння методів на згасаючому осциляторі
Згасаючий гармонічний осцилятор ÿ + 2γẏ + ω₀²y = 0 має точний розв'язок y(t) = e^{−γt}·cos(ωd·t), де ωd = √(ω₀²−γ²). Налаштуйте h (розмір кроку) та γ (згасання), щоб побачити, як кожен метод відслідковує справжній розв'язок. Велике h показує, як метод Ейлера розбігається, тоді як RK4 лишається точним. Абсолютна похибка енергії з часом показана на нижній панелі.
Точний
Чисельний
Похибка енергії
(масштабована)
Верле з γ=0 зберігає точну енергію; усі явні методи RK
накопичують дрейф енергії, пропорційний h^порядку.