Практично кожна фізична симуляція зводиться до інтегрування системи звичайних диференціальних рівнянь. Вибір чисельного методу визначає точність, стійкість та обчислювальну вартість — і хибний вибір для жорсткої системи може означати результати, що відрізняються на порядки величини або відверто розбігаються. Ця стаття оглядає ключові методи: від підручникового методу Ейлера через класичний 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^порядку.