Чисельне інтегрування для фізики — Ейлер, Верле та RK4
Фізичні симуляції розвивають свій стан у часі, апроксимуючи неперервні
рівняння дискретними часовими кроками. Вибір інтегратора визначає
стабільність, збереження енергії та точність. Цей посібник порівнює
чотири інтегратори на системі пружина-маса — від найпростішого прямого
методу Ейлера до еталонного RK4.
1Задача: інтегрування ЗДР
Фізичне тіло має стан: положення x та швидкість v.
Другий закон Ньютона дає прискорення: a = F/m. Інтегрування
означає просування стану вперед на малий часовий крок dt:
// Стан: положення x, швидкість v // Прискорення: a = F(x, v) / m //
Пружина-маса: F = -k * x - b * v (пружина + згасання) const k = 50; //
жорсткість пружини const b = 0.5; // коефіцієнт згасання const m = 1; //
маса function acceleration(x, v) { return (-k * x - b * v) / m; } //
Початкові умови let x = 1.0; // розтягнуто на 1 одиницю від рівноваги
let v = 0.0; // у спокої const dt = 0.016; // 16мс / кадр ≈ 60 к/с
2Прямий (явний) метод Ейлера — простий,
але нестабільний
Використовує поточну похідну для екстраполяції стану. Швидкий і
простий, але додає енергію системі — за великих
dt пружини зрештою «вибухають».
function stepEuler(x, v, dt) { const a = acceleration(x, v); const
xNew = x + v * dt; // положення використовує СТАРУ швидкість const vNew = v + a *
dt; // швидкість використовує СТАРЕ прискорення return [xNew, vNew]; } // Енергія
з часом дрейфує вгору — амплітуда зростає // СТАБІЛЬНО лише коли: dt
< 2 / sqrt(k/m) ≈ 0.28с для k=50, m=1
Прямий метод Ейлера має точність першого порядку (похибка ∝ dt). Для
пружини з k=50 і m=1 критичний часовий крок становить ~0.28с. Ігровий
цикл на 60 к/с (dt≈0.016с) безпечно вкладається в цю межу, але жорсткі
пружини за низької частоти кадрів вийдуть з-під контролю.
3Симплектичний Ейлер — зберігає енергію
Одна невелика зміна: спочатку оновлюємо швидкість за поточною силою, а
потім оновлюємо положення, використовуючи нову швидкість. Це
робить інтегратор симплектичним — він зберігає модифікований
гамільтоніан, і осцилятори не набирають енергії.
function stepSymplecticEuler(x, v, dt) { const a = acceleration(x, v);
const vNew = v + a * dt; // спочатку швидкість (НОВА) const xNew = x + vNew
* dt; // положення використовує НОВУ швидкість return [xNew, vNew]; } // Енергія
обмежена — пружина стабільно коливається завжди // Саме це більшість
ігрових фізичних рушіїв (Bullet, Box2D) використовують усередині // Така сама вартість, як у
прямого Ейлера — просто два рядки переставлено місцями!
4Інтегрування Верле — на основі положення
Швидкісний Верле зберігає попереднє положення замість явної швидкості,
даючи точність другого порядку та добре збереження енергії. Популярний
у симуляціях тканини, молекулярній динаміці (GROMACS) та динаміці на
основі положень.
// Швидкісний Верле (варіант «чехарда») function stepVerlet(x, v, dt) {
const a = acceleration(x, v); const xNew = x + v * dt + 0.5 * a * dt *
dt; const aNew = acceleration(xNew, v); // прискорення в новому
положенні const vNew = v + 0.5 * (a + aNew) * dt; // середнє
прискорення return [xNew, vNew]; } // Класичний Верле (без явної
швидкості — використовує x_prev): let xPrev = x - v * dt; // ініціалізуємо x_prev
function stepClassicVerlet(x, xPrev, dt) { const a = acceleration(x,
0); // лише сили, незалежні від швидкості const xNew = 2 * x - xPrev + a
* dt * dt; const vApprox = (xNew - xPrev) / (2 * dt); // центральна
різниця return [xNew, x, vApprox]; }
Класичний Верле не може напряму використовувати сили, що залежать від
швидкості (як-от згасання) — для систем зі згасанням застосовуйте
швидкісний Верле або варіант «чехарда».
5Рунге-Кутта 4 — еталонний стандарт
RK4 обчислює чотири оцінки нахилу за крок і бере зважене середнє.
Точність четвертого порядку означає, що зменшення dt удвічі
зменшує похибку у 16 разів. Ціною є чотири обчислення сили за крок.
function stepRK4(x, v, dt) { // Кожне k = [dx, dv] — пара похідних
const k1 = deriv(x, v ); const k2 = deriv(x + k1[0]*dt/2, v +
k1[1]*dt/2); const k3 = deriv(x + k2[0]*dt/2, v + k2[1]*dt/2); const
k4 = deriv(x + k3[0]*dt, v + k3[1]*dt ); const xNew = x + (dt/6) *
(k1[0] + 2*k2[0] + 2*k3[0] + k4[0]); const vNew = v + (dt/6) * (k1[1]
+ 2*k2[1] + 2*k3[1] + k4[1]); return [xNew, vNew]; } // Функція похідної
[ẋ = v, v̇ = a(x,v)] function deriv(x, v) { return [v,
acceleration(x, v)]; } // RK4 у 4 рази дорожчий за Ейлера, але
дозволяє значно більший dt // зберігаючи точність — ідеальний для
орбітальної механіки та жорстких систем
6Вибір інтегратора та фіксований
часовий крок
Інтегратор
Порядок
Енергія
Обчислень/крок
Сфера застосування
Прямий Ейлер
1-й
Набирає
1
Лише для навчання
Симплектичний Ейлер
1-й
Зберігає
1
Ігри, частинки
Швидкісний Верле
2-й
Зберігає
2
Тканина, тіла
RK4
4-й
~Зберігає
4
Орбіти, точність
Завжди використовуйте фіксований часовий крок з
акумулятором, щоб відокремити рендеринг від фізики:
const FIXED_DT = 1 / 120; // фізика на 120 Гц let accumulator = 0; let
prevTime = performance.now() / 1000; function loop(nowMs) { const now
= nowMs / 1000; const frameTime = Math.min(now - prevTime, 0.1); //
обмежуємо до 100мс (втрата фокусу вкладки) prevTime = now; accumulator += frameTime;
while (accumulator >= FIXED_DT) { [x, v] = stepSymplecticEuler(x, v,
FIXED_DT); accumulator -= FIXED_DT; } // Інтерполяція для плавного
рендерингу (alpha = залишкова частка) const alpha = accumulator /
FIXED_DT; const xRender = x; // або lerp(xPrev, x, alpha) для плавності
між кадрами mesh.position.x = xRender; requestAnimationFrame(loop); }
Шаблон з акумулятором — це галузевий стандарт (Gaffer On Games:
«Fix Your Timestep!»). Без нього фізика працює по-різному на 30 к/с та
144 к/с, і результати симуляції змінюються залежно від частоти
оновлення дисплея.