🔬 Числова стійкість: чому симуляції «вибухають» і як це виправити

Плаваюча кома: це не справжні числа

64-бітні числа з плаваючою комою (подвійна точність IEEE 754) представляють числа приблизно з 15–17 значущими десятковими цифрами точності. Вони розташовані нерівномірно: проміжок між представними числами поблизу 1,0 становить приблизно 2,2 × 10⁻¹⁶ (машинний епсилон), тоді як поблизу 10⁶ проміжок становить приблизно 2,2 × 10⁻¹⁰. Саме тому порядок операцій має значення для чисел навіть у «точній» арифметиці.

Класичний збій — це катастрофічне скорочення: віднімання двох майже рівних чисел призводить до втрати значущих цифр. Розгляньмо обчислення (1 + x) - 1 для x = 1e-16. У точній арифметиці це дорівнює x. У плаваючій комі 1 + 1e-16 округлюється точно до 1,0 (нижче машинного епсилону), тому результат дорівнює 0,0 — похибка 100%. Це проявляється в багатьох фізичних формулах:

// Погано: катастрофічне скорочення для малих x
const bad = (1 - Math.cos(x)) / (x * x);

// Добре: використовуємо тригонометричну тотожність
// (1 - cos x) = 2 sin²(x/2)
const good = 2 * Math.sin(x/2)**2 / (x * x);

Інші небезпеки плаваючої коми: підсумовування великої кількості чисел у наївному циклі накопичує похибки округлення O(n · ε). Сумування Кагана виправляє це, підтримуючи поточний компенсаційний доданок, зменшуючи накопичену похибку до O(ε) незалежно від n. Для фізичних симуляцій, що підсумовують сили від тисяч частинок, сумування Кагана може стати вирішальним фактором між стабільною енергією та її дрейфом.

Числа обумовленості: погано обумовлені задачі

Число обумовленості задачі вимірює, наскільки сильно змінюється результат відносно малих змін вхідних даних. Число обумовленості 100 означає, що відносна похибка входу 10⁻¹⁵ може призвести до відносної похибки виходу 10⁻¹³ — все ще прийнятно. Число обумовленості 10¹⁵ означає, що сама точність плаваючої коми робить результат безглуздим.

Розв'язання лінійної системи Ax = b має число обумовленості κ(A) = ||A|| · ||A⁻¹||. Майже вироджені матриці (де два рядки майже паралельні) мають величезні числа обумовленості. Саме тому пряме обертання матриці небажане для числових обчислень — натомість слід використовувати LU-розкладання, QR-факторизацію або ітеративні методи, що враховують структуру задачі.

Для симуляцій погана обумовленість часто виникає через екстремальні співвідношення параметрів. Симуляція рідини, що поєднує швидкі акустичні хвилі (що поширюються зі швидкістю 340 м/с) з повільним потоком (0,1 м/с), має число Маха 0,0003 і співвідношення часових масштабів акустика/потік 3400. Явна покрокова інтеграція за часом мусить розрізняти швидкі акустичні хвилі, навіть коли нас цікавить лише повільний потік — надзвичайно марнотратно. Число обумовленості пов'язаної системи пропорційне цьому співвідношенню.

Жорсткість: коли явні методи не спрацьовують

Жорстке звичайне диференціальне рівняння містить процеси з дуже різними часовими масштабами. Класичний приклад — хімічна кінетика: швидка реакція з часовим масштабом 10⁻⁹ с і повільна реакція з часовим масштабом 10³ с пов'язані між собою. Явний інтегратор мусить робити кроки за часом < 10⁻⁹ с, щоб залишатися стійким — мільярди кроків для моделювання однієї мілісекунди, навіть якщо швидкий процес практично миттєво досягає квазірівноваги і не додає корисної динаміки на повільному часовому масштабі.

Область стійкості явного методу — це обмежена область на комплексній площині. Для методу Ейлера стійкість вимагає |1 + Δt·λ| ≤ 1 для всіх власних значень λ якобіана. Для жорсткої системи з λ = −10⁶ стійкість вимагає Δt ≤ 2/10⁶ = 2 мкс. Явний RK4 має більшу область стійкості, але те саме фундаментальне обмеження.

Неявні методи включають майбутній стан в обчислення похідної, що робить умови стійкості значно менш обмежувальними. Метод оберненого Ейлера y(t+Δt) = y(t) + Δt · f(y(t+Δt)) вимагає розв'язання нелінійної системи на кожному кроці (методом Ньютона), але є A-стійким: стійким за будь-якого Δt для задач із власними значеннями від'ємної дійсної частини. Це дозволяє робити величезні кроки за часом для жорстких задач, часто роблячи неявні методи в 10–1000 разів швидшими, попри вищу вартість кожного кроку.

Діагностика нестійкої симуляції

Коли симуляція «вибухає», першопричина майже завжди одна з трьох:

• Крок за часом занадто великий: перевірте Δt · max(|власні значення|). Якщо це перевищує область стійкості вашого інтегратора, зменшіть Δt або перейдіть на неявний метод.
• Катастрофа плаваючої коми: логуйте проміжні значення й шукайте раптові стрибки чи поширення NaN. Поширені причини: ділення на близьку до нуля величину (відстані зіткнення частинок), квадратні корені від від'ємних чисел (шум плаваючої коми, через який величина, що мала бути нулем, стає трохи від'ємною), або логарифми чисел, які мають бути додатними.
• Дрейф енергії через несимплектичну інтеграцію: відстежуйте повну енергію з часом. Якщо вона монотонно зростає (а не коливається), у вас систематичне вливання енергії через інтегратор. Перейдіть на симплектичний метод.

// Захисний цикл симуляції
function step(dt) {
  const newState = integrate(state, dt);
  if (!isFinite(newState.energy) || newState.energy > MAX_ENERGY) {
    console.warn('Instability detected, halving timestep');
    dt /= 2;
    return step(dt);  // повторна спроба з меншим кроком
  }
  state = newState;
}

Практичні техніки забезпечення стійкості

Для симуляцій рідини адвекція поля швидкості напівлагранжевими методами (простеження шляхів частинок назад у часі) є безумовно стійкою за будь-якого кроку за часом — ціною певної числової дифузії. Умова Куранта-Фрідріхса-Леві (CFL) обмежує крок за часом для явної адвекції: Δt ≤ CFL · Δx / max(|u|). З CFL = 0,9 комірка рідини, перетнута за один крок за часом, залишається стійкою; CFL > 1 означає, що інформація поширюється більш ніж на одну комірку за крок, і настає нестійкість.

Для систем пружина-маса (тканина, м'які тіла) критичний крок за часом масштабується як Δt_crit ~ √(m/k): жорсткіші пружини (більше k) вимагають менших кроків за часом. Динаміка на основі позицій (PBD) уникає цього, проєктуючи обмеження після кожного кроку замість прямої інтеграції сил — жертвуючи фізичною точністю заради безумовної стійкості за будь-якого кроку за часом, що й пояснює її широке використання в ігровій фізиці.