📐 Euler czy Runge-Kutta: Którą metodę całkowania wybrać do swojej symulacji?

Problem całkowania

Większość symulacji fizycznych opisują równania różniczkowe zwyczajne (RRZ). Kąt wahadła θ spełnia d²θ/dt² = −(g/L)·sin(θ). Grawitacyjny układ N ciał spełnia drugie prawo Newtona dla każdego ciała. Płyn rządzi się równaniami Naviera-Stokesa. Aby symulować te układy, musimy całkować RRZ w czasie.

Podstawowe wyzwanie: komputery mogą reprezentować jedynie stany dyskretne. Wybieramy krok czasowy Δt i pytamy: „mając stan w chwili t, jaki jest stan w chwili t + Δt?" Odpowiedź zawsze jest przybliżeniem — a jakość tego przybliżenia zależy od metody całkowania. Złe wybory prowadzą do dryfu energii, niestabilności lub jawnej eksplozji numerycznej.

Metoda Eulera: prosta, ale niebezpieczna

Metoda Leonharda Eulera, sięgająca 1768 roku, jest najprostszym podejściem. Dla RRZ pierwszego rzędu dy/dt = f(y, t) krok naprzód wygląda tak:

y(t + Δt) ≈ y(t) + Δt · f(y(t), t)

To rozwinięcie Taylora pierwszego rzędu — wykorzystuje tylko pochodną w bieżącym punkcie do ekstrapolacji. Lokalny błąd obcięcia (błąd na krok) jest proporcjonalny do Δt², a globalny błąd zakumulowany w ustalonym przedziale czasu wynosi O(Δt) — stąd nazwa „metoda pierwszego rzędu".

Praktycznym problemem jest dryf energii. Dla układów zachowawczych, takich jak wahadła czy orbity, metoda Eulera systematycznie dodaje energię do układu. Wahadło symulowane metodą Eulera będzie powoli rozkręcać się na zewnątrz, nieskończenie zwiększając amplitudę. Zmniejszenie Δt spowalnia dryf — ale nigdy go nie zatrzymuje. To fundamentalna właściwość metody, a nie błąd, który można naprawić mniejszym krokiem czasowym.

Runge-Kutta 4. rzędu: koń roboczy

Metoda Rungego-Kutty czwartego rzędu (RK4) oblicza pochodną w czterech punktach wewnątrz każdego kroku czasowego i łączy je za pomocą średniej ważonej:

k1 = f(y,       t)
k2 = f(y + k1·Δt/2,  t + Δt/2)
k3 = f(y + k2·Δt/2,  t + Δt/2)
k4 = f(y + k3·Δt,    t + Δt)

y(t+Δt) = y(t) + (Δt/6)·(k1 + 2k2 + 2k3 + k4)

Lokalny błąd obcięcia wynosi O(Δt⁵), a błąd globalny O(Δt⁴). Podwojenie kroku czasowego mnoży błąd przez 16. Zmniejszenie go o połowę dzieli błąd przez 16. Ta dramatyczna wrażliwość na rozmiar kroku jest tym, co czyni RK4 tak potężną metodą — umiarkowane kroki czasowe zapewniają wysoką dokładność.

RK4 wymaga czterech obliczeń funkcji na krok (w porównaniu do jednego dla Eulera), więc kosztuje mniej więcej 4 razy więcej na krok. Ale ponieważ osiąga dokładność przy większych krokach czasowych, całkowity koszt osiągnięcia zadanego poziomu dokładności jest znacznie niższy. Dla większości gładkich problemów fizycznych RK4 jest domyślnym wyborem.

Problem dryfu energii: całkowania symplektyczne

RK4 jest zdecydowanie lepsza od Eulera, ale wciąż ma fundamentalną wadę dla zachowawczych układów mechanicznych: nie jest symplektyczna. Całkowanie symplektyczne zachowuje strukturę geometryczną mechaniki hamiltonowskiej — konkretnie zachowuje objętość przestrzeni fazowej (twierdzenie Liouville'a). Metody niesymplektyczne łamią tę własność i powodują wiekowy dryf energii.

Najprostszym całkowaniem symplektycznym jest metoda żabich skoków (leapfrog, znana też jako metoda Störmera-Verleta). Dla układu z pozycją x i prędkością v:

// Aktualizacja prędkości na pół kroku
v_half = v(t) + (Δt/2) · a(x(t))

// Aktualizacja pozycji na pełny krok
x(t+Δt) = x(t) + Δt · v_half

// Obliczenie przyspieszenia w nowej pozycji
a_new = f(x(t+Δt))

// Dokończenie kroku prędkości
v(t+Δt) = v_half + (Δt/2) · a_new

Leapfrog ma dokładność zaledwie drugiego rzędu — gorszą niż RK4 na krok. Ale ponieważ jest symplektyczna, zachowuje energię średnio na długich skalach czasowych. Energia oscyluje wokół prawdziwej wartości zamiast monotonicznie dryfować. Dla długotrwałych symulacji układów zachowawczych — orbit planetarnych, dynamiki molekularnej, akceleratorów cząstek — całkowania symplektyczne są niezbędne.

Symplektyczne metody wyższego rzędu

Forest i Ruth (1990) oraz Yoshida (1990) pokazali, jak konstruować całkowania symplektyczne wyższego rzędu, łącząc kroki leapfrog ze starannie dobranymi współczynnikami. Symplektyczne całkowanie czwartego rzędu wymaga w przybliżeniu 4–7 podkroków leapfrog na pełny krok, osiągając dokładność na poziomie RK4 przy zachowaniu struktury symplektycznej. Metody te są standardem w długoterminowych kodach dynamiki molekularnej i astronomii N ciał.

Macierz kompromisów:

• Euler: Szybka, pierwszego rzędu, z dryfem energii. Używaj tylko do szybkich prototypów lub układów niezachowawczych.
• RK4: Cztery obliczenia na krok, czwartego rzędu, z niewielkim dryfem energii. Najlepsza do dokładnych symulacji krótko- i średnioterminowych gładkich układów.
• Leapfrog/Verlet: Dwa obliczenia na krok, drugiego rzędu, zachowuje energię średnio. Najlepsza do długotrwałych symulacji układów zachowawczych.
• Symplektyczna 4. rzędu: 5–7 obliczeń na krok, czwartego rzędu, symplektyczna. Łączy zalety obu światów w długich, wymagających symulacjach.

Adaptacyjny rozmiar kroku

Metody o stałym kroku czasowym mają trudności, gdy rozwiązanie zmienia się szybko w jednych obszarach, a wolno w innych. Metody adaptacyjne szacują błąd lokalny, porównując rozwiązania różnych rzędów (np. RK4 kontra RK5), i dynamicznie dostosowują Δt. Metoda Dormanda-Prince'a (używana w funkcji ode45 MATLAB-a i solve_ivp SciPy) osiąga to zaledwie sześcioma obliczeniami funkcji. Rozmiar kroku jest zwiększany, gdy błąd jest mały (dla zysku prędkości), i zmniejszany, gdy błąd jest duży (dla zachowania dokładności).