Шпаргалка з чисельних методів

🎯 Навіщо потрібні чисельні методи

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

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

🌀 Інтегратори ОДР

Стан симуляції (позиція, швидкість) розвивається за диференціальним рівнянням першого порядку $\dot{\mathbf{y}} = f(t, \mathbf{y})$. Інтегратор просуває $\mathbf{y}$ від $t$ до $t+h$, використовуючи обчислення $f$.

Явний (прямий) метод Ейлера

Найпростіший інтегратор: одне обчислення функції на крок, тривіальний у реалізації, але систематично накопичує енергію в осцилюючих системах (маятник, змодельований методом Ейлера, повільно розкручується назовні). Локальна похибка $O(h^2)$, глобальна — $O(h)$.

function eulerStep(state, force, h) {
  // state = { x, v }, force(x, v) повертає прискорення
  const a = force(state.x, state.v);
  state.x += h * state.v;
  state.v += h * a;
  return state;
}

Напівнеявний (симплектичний) метод Ейлера

Поміняйте порядок — спершу оновіть швидкість, а потім використайте нову швидкість для оновлення позиції. Та сама вартість, що й у прямого Ейлера, все ще лише перший порядок точності, але метод симплектичний: енергія коливається навколо істинного значення замість того, щоб дрейфувати геть. Це стандартний інтегратор майже у кожному ігровому фізичному рушії реального часу (Box2D, Cannon-es, Rapier).

function symplecticEulerStep(state, force, h) {
  const a = force(state.x, state.v);
  state.v += h * a;        // спершу швидкість
  state.x += h * state.v;  // потім позиція, з новою v
  return state;
}

Швидкісний метод Верле (Velocity Verlet)

Робочий кінь для систем частинок, молекулярної динаміки та розв'язувачів тканини/канатів. Точність другого порядку, симплектичний (обмежена похибка енергії) та оборотний у часі. Коштує два обчислення прискорення на крок, коли прискорення залежить від швидкості; дешевше — коли ні.

function verletStep(state, accel, h) {
  const a0 = accel(state.x);
  state.x += h * state.v + 0.5 * h * h * a0;
  const a1 = accel(state.x);
  state.v += 0.5 * h * (a0 + a1);
  return state;
}

Рунге-Кутта 4-го порядку (RK4)

Чотири обчислення функції на крок дають точність четвертого порядку — глобальна похибка $O(h^4)$ проти $O(h)$ у Ейлера. Чудово підходить для гладкої, неколізійної динаміки (орбітальна механіка, траєкторії камери), де потрібна висока точність при більшому кроці. RK4 не симплектичний: за дуже довгих інтегрувань (мільйони орбіт) енергія все одно може дрейфувати, лише набагато повільніше, ніж у методі Ейлера.

function rk4Step(y, t, h, f) {
  const k1 = f(t, y);
  const k2 = f(t + h / 2, add(y, scale(k1, h / 2)));
  const k3 = f(t + h / 2, add(y, scale(k2, h / 2)));
  const k4 = f(t + h, add(y, scale(k3, h)));
  return add(y, scale(
    add(add(k1, scale(k2, 2)), add(scale(k3, 2), k4)),
    h / 6
  ));
}

⚖ Стійкість і крок інтегрування

Метод є стійким для заданого кроку $h$, якщо малі похибки не зростають необмежено протягом багатьох кроків. Для тестового рівняння $\dot{y} = \lambda y$ (де $\lambda$ комплексне, $\text{Re}(\lambda) \le 0$) кожен явний метод має обмежену область абсолютної стійкості у комплексній площині $h\lambda$.

Для демпфованої системи пружина-маса з жорсткістю $k$ і масою $m$ це перетворюється на практичне правило:

Перевищте це — і жорстка пружина (з великим k), змодельована явним інтегратором, почне коливатися з амплітудою, що зростає, і врешті "вибухне" — класична помилка "тканина, що трясеться і вибухає". Три способи виправлення, за пріоритетом:

• Фіксований малий крок часу — відв'яжіть частоту оновлення фізики від частоти кадрів рендерингу (наприклад, фізика на 60 Гц через акумулятор, навіть при 30 fps рендерингу).
• Підкроки (substepping) — виконуйте N менших кроків фізики на один кадр рендерингу замість одного великого.
• Неявне інтегрування — зворотний Ейлер або position-based dynamics (PBD) безумовно стійкі для лінійних жорстких доданків ціною розв'язання лінійної (або ітеративної) системи на кожному кроці.

// Цикл з фіксованим кроком і акумулятором (відв'язує фізику від частоти рендерингу)
let accumulator = 0;
const FIXED_DT = 1 / 60;

function frame(now) {
  accumulator += (now - lastTime) / 1000;
  lastTime = now;
  while (accumulator >= FIXED_DT) {
    physicsStep(FIXED_DT);
    accumulator -= FIXED_DT;
  }
  render(accumulator / FIXED_DT); // альфа інтерполяції
  requestAnimationFrame(frame);
}

🔍 Пошук коренів

Використовується для запитів часу зіткнення (time-of-impact), зворотної кінематики та розв'язання неявних рівнянь зв'язків $g(\mathbf{x}) = 0$.

Метод бісекції

Вимагає, щоб $f(a)$ і $f(b)$ мали протилежні знаки. Завжди збігається, ділячи навпіл інтервал на кожній ітерації — надійно, але повільно ($O(\log_2(\epsilon^{-1}))$ ітерацій для точності $\epsilon$).

Метод Ньютона-Рафсона

Збігається набагато швидше за бісекцію, коли взагалі збігається, але може розходитись або циклити, якщо початкове наближення далеко від кореня або $f'(x_n) \approx 0$. Активно використовується у неперервному виявленні зіткнень для уточнення точного часу удару.

function newtonRaphson(f, fPrime, x0, tol = 1e-6, maxIter = 50) {
  let x = x0;
  for (let i = 0; i < maxIter; i++) {
    const fx = f(x);
    if (Math.abs(fx) < tol) return x;
    const dfx = fPrime(x);
    if (Math.abs(dfx) < 1e-12) break; // похідна занадто плоска, вихід
    x -= fx / dfx;
  }
  return x; // найкраще наближення — виклик повинен перевірити результат
}

Бісекція проти Ньютона — коли що використовувати

🔢 Лінійні системи

Розв'язання $A\mathbf{x} = \mathbf{b}$ трапляється по всій симуляції: розв'язувачі зв'язків (контакти твердих тіл, канати), неявне інтегрування та якобіани зворотної кінематики.

Ітерація Гаусса-Зейделя

Фізичні рушії (Box2D, Bullet, Cannon-es) розв'язують контактні та шарнірні зв'язки жменькою ітерацій Гаусса-Зейделя (або проєктивного Гаусса-Зейделя для нерівностей, наприклад, неперетину) на кадр замість точного розв'язку — кілька ітерацій наближеної відповіді щокадру кращі, ніж один точний розв'язок, що не встигає в дедлайн.

// N проходів Гаусса-Зейделя по розрідженій системі зв'язків
function solveConstraints(constraints, iterations = 10) {
  for (let iter = 0; iter < iterations; iter++) {
    for (const c of constraints) {
      c.solve(); // оновлює швидкості на місці, використовує останні значення
    }
  }
}

LU-розклад (щільний, прямий)

Для невеликих щільних систем (наприклад, розв'язку тензора інерції 3×3 або 4×4) прямий LU-розклад дешевий і точний: факторизуйте $A = LU$ один раз, а потім розв'язуйте дві трикутні системи для кожної правої частини через пряму/зворотну підстановку — $O(n^3)$ на факторизацію, $O(n^2)$ на кожен розв'язок.

📈 Інтерполяція

Лінійна інтерполяція (lerp)

Кубічний сплайн Ерміта / Катмулла-Рома

Сплайн Катмулла-Рома проходить точно через кожну контрольну точку (на відміну від Безьє) і є стандартним вибором для траєкторій камери та процедурних доріг/трас у симуляціях на цьому сайті.

Сферична лінійна інтерполяція (slerp)

На відміну від покомпонентної lerp-інтерполяції кватерніонів, slerp інтерполює вздовж дуги великого кола, тому швидкість обертання лишається сталою — це критично важливо для плавного змішування орієнтації камери та персонажа.

∫ Чисельне диференціювання та інтегрування

Центральна різниця (похідна)

Використовується для чисельного наближення якобіана або нормального вектора (наприклад, градієнта поля знакованої відстані), коли аналітична похідна недоступна. Центральна різниця скасовує доданок похибки першого порядку, який несе наївна пряма різниця $\frac{f(x+h)-f(x)}{h}$.

Метод трапецій (інтегрування)

Метод Сімпсона (інтегрування)

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

🌊 Основи ШПФ

Швидке перетворення Фур'є (ШПФ) перетворює сигнал з часового/ просторового домену в частотний за $O(n\log n)$ замість $O(n^2)$ наївного дискретного перетворення Фур'є (ДПФ) — основа симуляції океанських хвиль (Tessendorf/FFT ocean), аудіореактивної візуалізації та розв'язувачів рідин на GPU.

Алгоритм Кулі-Тьюкі рекурсивно розбиває суму на доданки з парними та непарними індексами, зменшуючи розмір задачі вдвічі на кожному рівні рекурсії — тому $N$ зазвичай має бути степенем двійки для класичної реалізації radix-2, яка використовується у більшості WebGL-шейдерів океану.

На практиці океанські та водні симуляції на цьому сайті виконують ШПФ на compute-шейдері GPU (або через прохід ping-pong render target у WebGL), щоб перетворити спектр висотного поля на зміщення — сітка 256×256 коштує приблизно 16 проходів radix-2 "метеликів" (butterflies) на вимір.

📊 Таблиця порядків похибки

Зменшення кроку $h$ удвічі для методу порядку $p$ знижує похибку приблизно у $2^p$ разів. Використовуйте це, щоб емпірично перевірити збіжність методу (екстраполяція Річардсона): зменшіть $h$ удвічі, перезапустіть і переконайтеся, що похибка падає в очікуваному співвідношенні.

⚠ Типові пастки

• Змінний крок часу з методом Ейлера — запуск прямого Ейлера з довільним deltaTime від браузера змушує дрейф енергії залежати від частоти кадрів, тож одна й та сама симуляція поводиться по-різному на різних машинах.
• Катастрофічне скорочення (cancellation) — віднімання двох майже рівних чисел з рухомою комою (наприклад, b*b - 4*a*c у квадратній формулі, коли дискримінант ≈ 0) посилює похибку округлення; використовуйте чисельно стійку форму квадратної формули, коли точність важлива.
• Ігнорування жорсткості — одна дуже жорстка пружина або зв'язок в іншій м'якій системі змушує всю симуляцію використовувати крихітний крок (або призводить до нестійкості); ізолюйте жорсткі зв'язки і розв'язуйте їх неявно.
• Неправильний інтегратор для завдання — використання RK4 для твердих тіл, що стикаються, витрачає цикли даремно (він не симплектичний, а зіткнення все одно порушують припущення про гладкість); використання Ейлера для точної орбітальної механіки накопичує помітний дрейф.
• Аліасинг без низькочастотного фільтра — дискретизація висотного поля на основі ШПФ на занадто грубій сітці або відсутність віконної функції для сигналу вносить високочастотні артефакти.
• Порівняння чисел з рухомою комою на рівність — завжди використовуйте допуск епсилон (Math.abs(a - b) < 1e-6) замість a === b для будь-якої величини, отриманої ітеративними чисельними методами.