Диференціальні рівняння в симуляціях

Що таке диференціальне рівняння?

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

ЗДР (звичайне диференціальне рівняння) містить похідні за однією змінною, зазвичай часом. ДРЧП (диференціальне рівняння з частинними похідними) містить похідні за кількома змінними — часом і простором одночасно.

Частина 1 — Звичайні диференціальні рівняння (ЗДР)

Найпростіше ЗДР: експоненційне зростання і спад

ЗДР dy/dt = ky має точний розв'язок y(t) = y₀ · e^kt. Коли k < 0, воно описує радіоактивний розпад, виведення ліків з організму та спад популяції. Коли k > 0, воно описує складні відсотки та зростання бактерій. Кожен чисельний метод — це лише наближення до цього точного розв'язку.

Метод Ейлера — швидкий, але неточний

Метод Ейлера просуває стан на один малий крок часу Δt, використовуючи лише поточну похідну (нахил). Це найпростіший з можливих інтеграторів:

          // Метод Ейлера — глобальна похибка O(Δt)
          function
          eulerStep(state, dt) {
          const deriv =
          derivatives(state);
          // наприклад { vel: acc, pos: vel }
          return { pos: state.pos + deriv.vel * dt, vel:
          state.vel + deriv.acc * dt }; }
        

Проблема: похибка накопичується як O(Δt) на крок. Для коливальних систем на кшталт маятника метод Ейлера повільно додає енергію на кожному кроці — маятник розгойдується дедалі ширше, поки не "втече". Це можна побачити в симуляції маятника, увімкнувши режим "Ейлер".

Метод Рунге-Кутти 4-го порядку (RK4) — робочий кінь інтеграції

RK4 обчислює похідні у чотирьох точках всередині кроку часу і комбінує їх зі зваженим середнім. Похибка падає до O(Δt⁴) — у 10 000 разів точніше за метод Ейлера при тому самому розмірі кроку.

          // RK4 — глобальна похибка O(Δt⁴), приблизно у 4 рази дорожчий за
            метод Ейлера
          function
          rk4Step(state, dt) {
          const k1 =
          derivatives(state);
          const k2 =
          derivatives(add(state,
          scale(k1, dt/2)));
          const k3 =
          derivatives(add(state,
          scale(k2, dt/2)));
          const k4 =
          derivatives(add(state,
          scale(k3, dt)));

          // Зважене середнє: 1/6 k1 + 1/3 k2 + 1/3 k3 + 1/6 k4
          return add(state,
          scale( addAll(k1,
          scale(k2, 2),
          scale(k3, 2), k4), dt
          / 6
          )); }
        

Інтегрування Верле — збереження енергії для механіки

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

          // Швидкісний метод Верле — симплектичний, O(Δt²), але стабільний
            за енергією
          pos_new = pos + vel * dt + 0.5 * acc * dt *
          dt; vel_new = vel + 0.5 * (acc + acc_new) *
          dt;
        

Частина 2 — Рівняння з частинними похідними (ДРЧП)

Рівняння теплопровідності

Рівняння теплопровідності описує, як температура дифундує крізь середовище з часом: ∂T/∂t = α ∇²T. Просторовий лапласіан ∇²T (друга похідна у просторі) дискретизується методом скінченних різниць на сітці.

2D-рівняння теплопровідності — явна схема скінченних різниць

T[i][j] += α · dt / dx² · (T[i+1][j] + T[i-1][j] + T[i][j+1] + T[i][j-1] − 4·T[i][j])

Умова стабільності: α · dt / dx² ≤ 0.25

Хвильове рівняння

Звук, водяні хвилі та струни підкоряються рівнянню ∂²u/∂t² = c² ∇²u. На відміну від рівняння теплопровідності (перша похідна за часом — "без пам'яті"), друга похідна за часом у хвильовому рівнянні означає, що збурення подорожують, а не дифундують. Симуляція фігур Хладні розв'язує двовимірне хвильове рівняння на квадратній пластині.

          // 2D хвильове рівняння — явні скінченні різниці (хвилі у ванні,
            фігури Хладні)
          for (let i =
          1; i < N-1; i++)
          { for (let j =
          1; j < N-1; j++)
          { const laplacian = u_prev[i+1][j] + u_prev[i-1][j] + u_prev[i][j+1] + u_prev[i][j-1] -
          4 * u_prev[i][j]; u_next[i][j] =
          2 * u_curr[i][j] - u_old[i][j] + c2 * dt2 *
          laplacian; } }
        

ДРЧП реакції-дифузії

Система Грея-Скотта (показана у Хімічному фокусі) — це система ДРЧП, у якій дифузія (доданки лапласіана) поєднується з нелінійними реакційними доданками. Просторовий лапласіан обчислюється так само, як і в рівнянні теплопровідності — саме реакційний доданок створює нестійкість Тюрінга та формування візерунків.

Вибір правильного методу

Метод	Порядок похибки	Вартість	Найкраще для
Ейлер	O(Δt)	1 обчислення/крок	Прототипування; ніколи не для промислової фізики
Верле	O(Δt²)	1–2 обчислення/крок	Системи частинок, молекулярна динаміка, тканина
RK4	O(Δt⁴)	4 обчислення/крок	Орбітальна механіка, ЗДР, що потребують точності
МСР (явний)	O(Δt, Δx²)	O(N²) на крок	Рівняння теплопровідності, хвильове рівняння на фіксованих сітках
МСР (неявний)	O(Δt², Δx²)	Розв'язання лінійної системи	Жорсткі ДРЧП, коли потрібні великі кроки часу

Де ці методи з'являються на сайті

Подвійні зорі / Сонячний вітрило

RK4 інтегрує гравітаційні ЗДР. Дрейф енергії за 10 000 орбіт становить менш ніж 0.01% — неможливо для методу Ейлера.

Симуляція тканини

Інтегрування Верле над сіткою пружина-маса. Обмеження на основі позиції обробляють зіткнення без корекцій швидкості.

Хвилі у ванні

Явний МСР хвильового рівняння на сітці 256×256. Стабільність забезпечується умовою Куранта-Фрідріхса-Леві: c·Δt/Δx ≤ 1.

Реакція-дифузія

Пов'язані ДРЧП Грея-Скотта на текстурі 512×512. Скінченні різниці обчислюються у фрагментному шейдері GLSL для паралелізму на GPU.

Дифузія ліків

Другий закон Фіка (варіант рівняння теплопровідності) розв'язаний явним МСР. Умова стабільності: D·Δt/Δx² ≤ 0.5.

Фігури Хладні

Задача власних значень двовимірного хвильового рівняння, розв'язана аналітично для форм мод, потім відтворена через відкладення частинок на вузлових лініях.

Ключова ідея: Кожен чисельний метод — це наближення до одного й того самого базового диференціального рівняння. Питання завжди одне: скільки точності вам потрібно і скільки обчислень ви можете собі дозволити? RK4 не завжди є відповіддю — для тривалих симуляцій частинок симплектична властивість методу Верле переважає точність RK4.