🔬 Stabilność numeryczna: dlaczego symulacje eksplodują i jak temu zapobiec

Liczby zmiennoprzecinkowe: to nie liczby rzeczywiste

64-bitowe liczby zmiennoprzecinkowe (podwójna precyzja IEEE 754) reprezentują liczby z około 15–17 znaczącymi cyframi dziesiętnymi precyzji. Nie są rozmieszczone równomiernie: odstęp między liczbami reprezentowalnymi w pobliżu 1,0 wynosi około 2,2 × 10⁻¹⁶ (epsilon maszynowy), podczas gdy w pobliżu 10⁶ odstęp wynosi około 2,2 × 10⁻¹⁰. Dlatego kolejność wykonywania działań ma znaczenie numeryczne nawet w przypadku „dokładnej” arytmetyki.

Klasyczną usterką jest katastrofalne skracanie: odejmowanie dwóch niemal równych liczb powoduje utratę znaczących cyfr. Rozważmy obliczenie (1 + x) - 1 dla x = 1e-16. W arytmetyce dokładnej wynik ten równa się x. W arytmetyce zmiennoprzecinkowej 1 + 1e-16 zaokrągla się dokładnie do 1,0 (poniżej epsilonu maszynowego), więc wynikiem jest 0,0 — błąd rzędu 100%. Ten problem pojawia się w wielu wzorach fizycznych:

// Źle: katastrofalne skracanie dla małych x
const bad = (1 - Math.cos(x)) / (x * x);

// Dobrze: użyj tożsamości trygonometrycznej
// (1 - cos x) = 2 sin²(x/2)
const good = 2 * Math.sin(x/2)**2 / (x * x);

Inne zagrożenia zmiennoprzecinkowe: sumowanie wielu liczb w naiwnej pętli kumuluje błędy zaokrąglenia rzędu O(n · ε). Sumowanie Kahana koryguje ten problem, utrzymując bieżący wyraz kompensujący, redukując skumulowany błąd do O(ε) niezależnie od n. W symulacjach fizycznych sumujących siły od tysięcy cząstek sumowanie Kahana może zadecydować o różnicy między stabilną a dryfującą energią.

Liczby uwarunkowania: zadania źle postawione

Liczba uwarunkowania zadania mierzy, jak bardzo zmienia się wynik względem niewielkich zmian danych wejściowych. Liczba uwarunkowania równa 100 oznacza, że względny błąd wejściowy 10⁻¹⁵ może wywołać względny błąd wyjściowy 10⁻¹³ — wciąż w porządku. Liczba uwarunkowania równa 10¹⁵ oznacza, że sama precyzja zmiennoprzecinkowa czyni wynik bezsensownym.

Rozwiązanie układu liniowego Ax = b ma liczbę uwarunkowania κ(A) = ||A|| · ||A⁻¹||. Macierze niemal osobliwe (gdzie dwa wiersze są prawie równoległe) mają ogromne liczby uwarunkowania. Dlatego bezpośrednie odwracanie macierzy jest niewskazane w obliczeniach numerycznych — zamiast tego stosuje się rozkład LU, faktoryzację QR lub metody iteracyjne wykorzystujące strukturę zadania.

W symulacjach złe uwarunkowanie często wynika ze skrajnych proporcji parametrów. Symulacja płynu sprzęgająca szybkie fale akustyczne (rozchodzące się z prędkością 340 m/s) z wolnym przepływem (0,1 m/s) ma liczbę Macha równą 0,0003 i stosunek skali czasowej akustyka-do-przepływu równy 3400. Jawne całkowanie czasowe musi rozdzielić szybkie fale akustyczne, nawet jeśli interesuje nas tylko wolny przepływ — co jest ogromnie marnotrawne. Liczba uwarunkowania sprzężonego układu jest proporcjonalna do tego stosunku.

Sztywność: kiedy metody jawne zawodzą

Sztywne równanie różniczkowe zawiera procesy o skrajnie różnych skalach czasowych. Kinetyka chemiczna to klasyczny przykład: szybka reakcja o skali czasowej 10⁻⁹ s sprzężona jest z wolną reakcją o skali czasowej 10³ s. Jawny integrator musi wykonywać kroki czasowe < 10⁻⁹ s, by pozostać stabilnym — miliardy kroków, by zasymulować milisekundę, mimo że szybki proces osiąga niemal natychmiast quasi-równowagę i nie wnosi żadnej użytecznej dynamiki w wolnej skali czasowej.

Obszar stabilności metody jawnej jest ograniczonym obszarem na płaszczyźnie zespolonej. Dla metody Eulera stabilność wymaga |1 + Δt·λ| ≤ 1 dla wszystkich wartości własnych λ jakobianu. Dla sztywnego układu z λ = −10⁶ stabilność wymaga Δt ≤ 2/10⁶ = 2 μs. Jawna metoda RK4 ma większy obszar stabilności, ale to samo fundamentalne ograniczenie.

Metody uwikłane uwzględniają przyszły stan w obliczeniu pochodnej, co czyni warunki stabilności znacznie mniej restrykcyjnymi. Metoda Eulera wstecz y(t+Δt) = y(t) + Δt · f(y(t+Δt)) wymaga rozwiązania nieliniowego układu na każdym kroku (przy użyciu metody Newtona), ale jest A-stabilna: stabilna dla dowolnego Δt w przypadku zadań z wartościami własnymi o ujemnej części rzeczywistej. Pozwala to na ogromne kroki czasowe dla zadań sztywnych, co często czyni metody uwikłane 10–1000 razy szybszymi mimo wyższego kosztu na krok.

Diagnozowanie niestabilnej symulacji

Kiedy symulacja eksploduje, przyczyna niemal zawsze jest jedną z trzech rzeczy:

• Zbyt duży krok czasowy: sprawdź Δt · max(|wartości własne|). Jeśli przekracza to obszar stabilności twojego integratora, zmniejsz Δt lub przełącz się na metodę uwikłaną.
• Katastrofa zmiennoprzecinkowa: loguj wartości pośrednie i szukaj nagłych skoków lub propagacji NaN. Częste sprawcy: dzielenie przez wartość bliską zeru (odległości kolizji cząstek), pierwiastki kwadratowe z liczb ujemnych (szum zmiennoprzecinkowy powodujący, że wartość, która powinna być zerem, jest nieznacznie ujemna) lub logarytmy liczb, które powinny być dodatnie.
• Dryf energii wskutek niesymplektycznego całkowania: monitoruj całkowitą energię w czasie. Jeśli rośnie monotonicznie (bez oscylacji), masz systematyczny wtrysk energii z twojego integratora. Przełącz się na metodę symplektyczną.

// Defensywna pętla symulacji
function step(dt) {
  const newState = integrate(state, dt);
  if (!isFinite(newState.energy) || newState.energy > MAX_ENERGY) {
    console.warn('Wykryto niestabilność, zmniejszanie kroku czasowego o połowę');
    dt /= 2;
    return step(dt);  // ponowna próba z mniejszym krokiem
  }
  state = newState;
}

Praktyczne techniki stabilizacji

W symulacjach płynów adwekcja pola prędkości metodami semi-Lagrange'a (śledzenie ścieżek cząstek wstecz w czasie) jest bezwarunkowo stabilna dla dowolnego kroku czasowego — kosztem pewnej dyfuzji numerycznej. Warunek Couranta-Friedrichsa-Lewy'ego (CFL) ogranicza krok czasowy dla jawnej adwekcji: Δt ≤ CFL · Δx / max(|u|). Przy CFL = 0,9 komórka płynu przekraczana w jednym kroku czasowym pozostaje stabilna; CFL > 1 oznacza, że informacja propaguje się przez więcej niż jedną komórkę na krok, co prowadzi do niestabilności.

W układach masa-sprężyna (tkaniny, ciała miękkie) krytyczny krok czasowy skaluje się jako Δt_krytyczne ~ √(m/k): sztywniejsze sprężyny (wyższe k) wymagają mniejszych kroków czasowych. Dynamika oparta na pozycjach (PBD) unika tego, projektując więzy po każdym kroku zamiast bezpośrednio całkować siły — poświęcając dokładność fizyczną na rzecz bezwarunkowej stabilności przy dowolnym kroku czasowym, dlatego jest powszechnie stosowana w fizyce gier.