Równania różniczkowe w symulacjach

Czym jest równanie różniczkowe?

Równanie różniczkowe to równanie wiążące funkcję z jej własnymi pochodnymi. W fizyce „funkcją" jest zwykle położenie, prędkość, temperatura lub stężenie substancji chemicznej, a „pochodna" to tempo jej zmiany w czasie lub przestrzeni.

ODE (równanie różniczkowe zwyczajne) zawiera pochodne względem jednej zmiennej, zwykle czasu. PDE (równanie różniczkowe cząstkowe) zawiera pochodne względem wielu zmiennych — czasu i przestrzeni jednocześnie.

Część 1 — Równania różniczkowe zwyczajne (ODE)

Najprostsze ODE: wzrost i zanik wykładniczy

Równanie ODE dy/dt = ky ma dokładne rozwiązanie y(t) = y₀ · e^kt. Gdy k < 0, opisuje ono rozpad promieniotwórczy, eliminację leku z organizmu i spadek populacji. Gdy k > 0, opisuje procent składany i wzrost bakterii. Każda metoda numeryczna jest jedynie przybliżeniem tego dokładnego rozwiązania.

Metoda Eulera — szybka, ale niedokładna

Metoda Eulera przesuwa stan o jeden mały krok czasowy Δt, wykorzystując wyłącznie bieżącą pochodną (nachylenie). To najprostszy możliwy integrator:

          // Metoda Eulera — błąd globalny O(Δt)
          function
          eulerStep(state, dt) {
          const deriv =
          derivatives(state);
          // np. { vel: acc, pos: vel }
          return { pos: state.pos + deriv.vel * dt, vel:
          state.vel + deriv.acc * dt }; }
        

Problem: błąd narasta jako O(Δt) na krok. Dla układów oscylujących, takich jak wahadło, metoda Eulera powoli dodaje energię przy każdym kroku — wahadło rozkołysuje się coraz szerzej, aż w końcu „ucieka". Można to zaobserwować w symulacji wahadła po włączeniu trybu „Euler".

Metoda Rungego-Kutty 4. rzędu (RK4) — koń roboczy integracji

RK4 oblicza pochodne w czterech punktach wewnątrz kroku czasowego i łączy je za pomocą średniej ważonej. Błąd spada do O(Δt⁴) — 10 000 razy dokładniej niż metoda Eulera przy tym samym rozmiarze kroku.

          // RK4 — błąd globalny O(Δt⁴), około 4× droższy niż metoda
            Eulera
          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)));

          // Średnia ważona: 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
          )); }
        

Całkowanie Verleta — zachowujące energię dla mechaniki

Całkowanie Verleta jest preferowaną metodą dla symulacji cząstek. Jest symplektyczne — zachowuje zmodyfikowaną energię przy długich całkowaniach — dlatego symulacje tkanin, ciał miękkich i dynamiki molekularnej używają go zamiast RK4:

          // Prędkościowy Verlet — symplektyczny, O(Δt²), ale stabilny
            energetycznie
          pos_new = pos + vel * dt + 0.5 * acc * dt *
          dt; vel_new = vel + 0.5 * (acc + acc_new) *
          dt;
        

Część 2 — Równania różniczkowe cząstkowe (PDE)

Równanie przewodnictwa cieplnego

Równanie przewodnictwa cieplnego opisuje, jak temperatura dyfunduje przez ośrodek w czasie: ∂T/∂t = α ∇²T. Przestrzenny laplasjan ∇²T (druga pochodna w przestrzeni) jest dyskretyzowany metodą różnic skończonych na siatce.

2D równanie przewodnictwa cieplnego — jawny schemat różnic skończonych

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])

Warunek stabilności: α · dt / dx² ≤ 0.25

Równanie falowe

Dźwięk, fale wodne i struny podlegają równaniu ∂²u/∂t² = c² ∇²u. W przeciwieństwie do równania przewodnictwa cieplnego (pierwsza pochodna czasowa — „bez pamięci"), druga pochodna czasowa w równaniu falowym oznacza, że zaburzenia przemieszczają się, zamiast dyfundować. Symulacja figur Chladniego rozwiązuje dwuwymiarowe równanie falowe na kwadratowej płycie.

          // 2D równanie falowe — jawne różnice skończone (fale w wannie,
            figury Chladniego)
          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; } }
        

PDE reakcji-dyfuzji

Układ Graya-Scotta (widoczny w materiale „W centrum uwagi: chemia") to układ PDE, w którym dyfuzja (człony laplasjanu) łączy się z nieliniowymi członami reakcji. Przestrzenny laplasjan obliczany jest tak samo jak w równaniu przewodnictwa cieplnego — to właśnie człon reakcji tworzy niestabilność Turinga i powstawanie wzorów.

Wybór odpowiedniej metody

Metoda	Rząd błędu	Koszt	Najlepsza do
Euler	O(Δt)	1 obliczenie/krok	Prototypowanie; nigdy do fizyki produkcyjnej
Verlet	O(Δt²)	1–2 obliczenia/krok	Systemy cząstek, dynamika molekularna, tkaniny
RK4	O(Δt⁴)	4 obliczenia/krok	Mechanika orbitalna, ODE wymagające dokładności
MRS (jawna)	O(Δt, Δx²)	O(N²) na krok	Równanie ciepła, równanie falowe na stałych siatkach
MRS (niejawna)	O(Δt², Δx²)	Rozwiązanie układu liniowego	Sztywne PDE, gdy potrzebne są duże kroki czasowe

Gdzie te metody pojawiają się na stronie

Podwójne gwiazdy / Żagiel słoneczny

RK4 całkuje grawitacyjne ODE. Dryf energii po 10 000 orbitach wynosi <0,01% — niemożliwe przy metodzie Eulera.

Symulacja tkaniny

Całkowanie Verleta na siatce sprężyna-masa. Ograniczenia oparte na pozycji obsługują kolizje bez korekt prędkości.

Fale w wannie

Jawna metoda różnic skończonych dla równania falowego na siatce 256×256. Stabilność wymuszona warunkiem CFL: c·Δt/Δx ≤ 1.

Reakcja-dyfuzja

Sprzężone PDE Graya-Scotta na teksturze 512×512. Różnice skończone obliczane są w shaderze fragmentów GLSL dla równoległości na GPU.

Dyfuzja leku

Drugie prawo Ficka (wariant równania ciepła) rozwiązane jawną metodą różnic skończonych. Warunek stabilności: D·Δt/Δx² ≤ 0,5.

Figury Chladniego

Zagadnienie własne dwuwymiarowego równania falowego rozwiązane analitycznie dla kształtów modów, następnie wizualizowane poprzez osadzanie cząstek na liniach węzłowych.

Kluczowy wniosek: Każda metoda numeryczna jest przybliżeniem tego samego, leżącego u podstaw równania różniczkowego. Pytanie zawsze brzmi: ile dokładności potrzebujesz i na ile obliczeń możesz sobie pozwolić? RK4 nie zawsze jest odpowiedzią — dla długich symulacji cząstek symplektyczna właściwość metody Verleta przewyższa dokładność RK4.