Порівняння методів симуляції — інтегратори, моделі рідин і Монте-Карло

1. Таблиця порівняння інтеграторів

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

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

2. Явний Ейлер

Найпростіший можливий інтегратор: обчислити похідну у поточному стані, а потім зробити крок вперед, використовуючи цю одну оцінку нахилу.

// Явний (прямий) Ейлер — НЕ використовуйте для будь-чого коливального
function stepEuler(body, dt) {
  const a = body.force.scale(body.invMass);
  body.pos = body.pos.add(body.vel.scale(dt)); // використовує СТАРУ швидкість
  body.vel = body.vel.add(a.scale(dt));
}

Оскільки позиція використовує швидкість до застосування сили, явний Ейлер систематично додає енергію до будь-якої коливальної системи (пружини, орбіти, маятники). Пружина без тертя, симульована явним Ейлером, буде помітно розкручуватися назовні й зрештою «вибухне», незалежно від того, наскільки малим зроблено dt — похибка лише зменшується, але якісно не зникає. Цей метод наведено тут суто як базовий рівень, який покращують усі інші.

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

Виправлення в один рядок: спочатку оновити швидкість, а потім використати нову швидкість для оновлення позиції. Та сама вартість, що й у явного Ейлера, але значно краща поведінка енергії.

function stepSemiImplicit(body, dt) {
  const a = body.force.scale(body.invMass);
  body.vel = body.vel.add(a.scale(dt));      // спочатку швидкість
  body.pos = body.pos.add(body.vel.scale(dt)); // потім позиція, з новою v
}

4. Верле та швидкісний Верле

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

function stepVerlet(p, dt) {
  const a = p.force.scale(p.invMass);
  const next = p.pos
    .scale(2)
    .sub(p.prevPos)
    .add(a.scale(dt * dt));
  p.prevPos = p.pos;
  p.pos = next;
}

Швидкісний Верле зберігає явну швидкість, даючи ту саму точність другого порядку, але з придатною для використання v(t) для демпфування, опору чи сил, залежних від швидкості:

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

RK4 обчислює похідну чотири рази за крок — на початку, двічі в середині та один раз наприкінці — і комбінує їх у зважене середнє. Це метод четвертого порядку точності: зменшення dt вдвічі зменшує локальну похибку приблизно у 16 разів, порівняно з 4 разами для Верле та 2 разами для Ейлера.

function rk4Step(state, dt, derivative) {
  const k1 = derivative(state);
  const k2 = derivative(addScaled(state, k1, dt / 2));
  const k3 = derivative(addScaled(state, k2, dt / 2));
  const k4 = derivative(addScaled(state, k3, dt));
  return addScaled(
    state,
    combine(k1, k2, k3, k4), // k1 + 2k2 + 2k3 + k4
    dt / 6
  );
}

Симуляції цього сайту Атрактор Лоренца та Подвійний маятник використовують RK4, бо хаотичні системи експоненційно посилюють малі числові похибки — похибка інтегратора низького порядку призводить до візуально неправильних траєкторій уже за кілька секунд. Недолік RK4: він не є симплектичним, тому орбітальні симуляції з RK4 повільно набирають або втрачають енергію протягом тисяч орбіт (для N-тіл, які мають працювати нескінченно довго, натомість використовуйте Leapfrog чи Верле).

6. Неявний Ейлер і жорсткі системи

Неявний (зворотний) Ейлер обчислює силу в майбутньому стані замість поточного, що вимагає розв'язання рівняння (часто розрідженої лінійної системи) на кожному кроці, а не просто підстановки чисел.

Ця додаткова робота дає безумовну стійкість: навіть з дуже великим dt неявний інтегратор не «вибухне» — він просто гасить високочастотний рух, що зазвичай є перевагою. Ця властивість найбільш важлива для жорстких систем, де сили змінюються за дуже різними часовими масштабами (наприклад, дуже натягнута пружина тканини поруч із розхлябаною). Стаття Барафа та Віткіна 1998 року «Large Steps in Cloth Simulation» популяризувала неявне інтегрування для тканини саме тому, що явні методи потребували неприпустимо малих кроків часу для збереження стійкості на жорстких пружинах.

7. Лагранжева vs ейлерова симуляція рідини

Окрім інтеграторів окремих частинок, суцільні середовища (рідини, дим, ерозія рельєфу) можна моделювати з двох протилежних точок зору.

Лагранжевий підхід: стежити за матеріалом

Представити рідину як набір рухомих частинок, кожна з яких несе масу, швидкість і густину. Smoothed Particle Hydrodynamics (SPH), що використовується у симуляції цього сайту SPH-рідина, є класичним прикладом: властивості кожної частинки обчислюються шляхом згладжування (зважування ядром) за сусідніми частинками.

SPH природно зберігає масу (частинки Є масою), обробляє вільні поверхні та бризки без додаткового обліку, а його вартість масштабується з кількістю частинок, а не з розміром порожнього простору. Недолік — пошук сусідів O(N·k) на кожному кроці та шумніші поверхні при малій кількості частинок.

Ейлерів підхід: зафіксувати сітку, дозволити матеріалу текти крізь неї

Дискретизувати простір у фіксовану сітку та зберігати швидкість/густину/тиск у кожній комірці, розв'язуючи рівняння Нав'є-Стокса методом скінченних різниць. Метод «Stable Fluids» Йоса Стама (напівлагранжева адвекція + проєкція тиску) є стандартною технікою реального часу:

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

Сіткові методи дають гладкі поля без артефактів і чудово підходять для диму, газів і патернів реакції-дифузії (див. Реакція-дифузія на цьому сайті, яка використовує фіксовану сітку з пінг-понг GPU текстурою). Їхня слабкість — вартість, прив'язана до роздільної здатності: сітка 256×256 витрачає однакові обчислення на порожнє повітря та на цікаву зону рідини — а також втрата дрібних деталей (маленьких крапель, тонких плівок) менших за розмір однієї комірки.

8. Методи Монте-Карло

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

Ця швидкість збіжності 1/√N є водночас найбільшою силою та найбільшою слабкістю Монте-Карло: вона однакова, скільки б вимірів не мала задача, тому Монте-Карло перевершує сіткову квадратуру для інтегралів високої розмірності (трасування шляхів світла, ціноутворення фінансових опціонів), але сходиться повільніше за гарний детермінований метод при низькій розмірності.

• Пряма вибірка: генерувати рівномірні випадкові точки та усереднювати — використовується у симуляції Генетичний алгоритм цього сайту для мутації та вибору схрещування.
• Марковський ланцюг Монте-Карло (MCMC): вибірка з складних розподілів ймовірності шляхом випадкового блукання, яке переважно відвідує зони високої ймовірності (Метрополіс- Гастінгс). Використовується у статистичній механіці та баєсівському висновуванні.
• Трасування шляхів методом Монте-Карло: у рендерингу випадково вибирають шляхи світла, що відбиваються по сцені, для розв'язання рівняння рендерингу — характерний «зашумлений, потім чіткий» вигляд трасованого рендеру — це саме видима збіжність 1/√N.

9. Як обрати метод

• Ігрова/частинкова фізика реального часу: напівнеявний Ейлер, за бажанням із 2–4 підкроками за кадр для додаткової стійкості — найдешевший варіант із прийнятною поведінкою енергії.
• Тканина, канати, м'які тіла: позиційний Верле з релаксацією обмежень, або неявні/PBD-солвери, якщо матеріал жорсткий і не повинен «вибухати» при великих кроках часу.
• Орбіти, N тіл, усе, що має працювати нескінченно довго без набору енергії: Leapfrog чи швидкісний Верле — симплектичні інтегратори тут правильний вибір за замовчуванням, а не RK4.
• Хаотичні системи, де точність траєкторії важливіша за довгостроковий баланс енергії (Лоренц, подвійний маятник): RK4, або адаптивний RK45, якщо потрібен контроль похибки.
• Бризки рідини, вільні поверхні: SPH (лагранжевий підхід).
• Дим, газ, патерни реакції-дифузії: сіткові ейлерові солвери (Stable Fluids, пінг-понг шейдери).
• Оцінка високої розмірності, рендеринг, пошук оптимізації: Монте-Карло / MCMC / генетичні алгоритми.