Porada: Wizualizacja ODE i PDE w Przeglądarce

Równania różniczkowe opisują, jak rzeczy się zmieniają — prędkość planety, temperaturę na metalowej płycie, wir płynu. Przeglądarka nie musi ich rozwiązywać symbolicznie w czasie rzeczywistym, ale nie musi tego robić. Wystarczy, że przesunie stan o mały krok czasowy, narysuje wynik i powtórzy. Zrób tę pętlę dobrze, a statyczne równanie zamienia się w żywy, przesuwalny, regulowalny obraz. Zrób to źle, a twoja symulacja albo pełznie, albo eksploduje do NaN w ciągu kilku sekund.

Poniżej znajduje się przepływ pracy, którego używamy w całym katalogu, od zwyczajnych równań różniczkowych (ODE) śledzących garstkę zmiennych po równania różniczkowe cząstkowe (PDE) ewoluujące całą siatkę wartości.

1. Wybierz Odpowiedni Integrator: Euler vs RK4

Solver ODE odpowiada wielokrotnie na jedno pytanie: mając bieżący stan i tempo zmiany, jaki będzie stan chwilę później? Najprostszą odpowiedzią jest metoda Eulera do przodu — weź bieżące nachylenie i zrób krok prosto wzdłuż niego.

// state = {x, v}, derivs(s) zwraca {dx, dv}
function eulerStep(s, dt) {
  const d = derivs(s);
  return { x: s.x + d.dx * dt, v: s.v + d.dv * dt };
}

Euler jest tani i całkowicie wystarczający dla łagodnych, wolnych systemów lub chmur cząstek, gdzie niewielki dryf jest niewidoczny. Ale szybko kumuluje błąd, a dla systemów oscylujących lub chaotycznych wstrzykuje fałszywą energię, aż orbita rozkręci się na zewnątrz. Tu właśnie zarabia na siebie metoda Rungego–Kutty czwartego rzędu (RK4). RK4 próbkuje nachylenie cztery razy na krok — raz na początku, dwa razy w środku, raz na końcu — i miesza je w znacznie dokładniejszą estymatę.

function rk4Step(s, dt) {
  const k1 = derivs(s);
  const k2 = derivs(add(s, scale(k1, dt / 2)));
  const k3 = derivs(add(s, scale(k2, dt / 2)));
  const k4 = derivs(add(s, scale(k3, dt)));
  return add(s, scale(add4(k1, k2, k2, k3, k3, k4), dt / 6));
}

RK4 kosztuje cztery ewaluacje pochodnej zamiast jednej, ale pozwala robić znacznie większe kroki czasowe przy tej samej dokładności — zazwyczaj to opłacalny kompromis. Nasz atraktor Lorenza w 3D polega na RK4: kształt motyla pozostaje wyraźny tylko dlatego, że integrator nie traci energii. Dla systemów sztywnych lub zachowujących energię wariant symplektyczny, jak velocity Verlet, jest często lepszym wyborem, ale RK4 jest bezpiecznym domyślnym wyborem, gdy nie masz pewności.

2. Siatki Różnic Skończonych dla PDE

PDE ewoluują pole — wartość w każdym punkcie przestrzeni — a nie kilka skalarów. Standardowym, przyjaznym przeglądarce podejściem jest metoda różnic skończonych: przechowuj pole jako płaski Float32Array indeksowany jako i = y * width + x, i przybliżaj pochodne przestrzenne, porównując każdą komórkę z jej sąsiadami.

Koniem roboczym jest dyskretny Laplasjan, który napędza dyfuzję, przepływ ciepła i propagację fal:

// 5-punktowy Laplasjan na siatce width×height
const i = y * width + x;
const lap = field[i - 1] + field[i + 1]
          + field[i - width] + field[i + width]
          - 4 * field[i];
next[i] = field[i] + alpha * lap * dt;

Zawsze pisz do osobnego bufora wyjściowego i zamień oba tablice po pełnym przejściu — aktualizowanie w miejscu psuje sąsiadów, których jeszcze nie odwiedziłeś. Ten wzorzec podwójnego buforowania jest dokładnie tym, co napędza symulacje siatkowe, takie jak niestabilność plazmy, gdzie pole naładowanego płynu jest krokowane komórka po komórce, oraz wygładzanie w stylu dyfuzji stojące za ogrzewaniem granularnym, gdzie lokalne zderzenia rozprowadzają energię przez upakowane złoże ziaren.

3. Szanuj Warunek Stabilności

Oto pułapka, która łapie każdego za pierwszym razem: jawne schematy różnic skończonych są stabilne tylko wtedy, gdy krok czasowy jest wystarczająco mały w stosunku do odstępu siatki. Dla równania dyfuzji w dwóch wymiarach warunek CFL (Couranta–Friedrichsa–Lewy'ego) z grubsza wymaga alpha * dt / dx² < 0,25. Przekrocz go, a pole nie tylko traci na dokładności — oscyluje dziko i przekracza zakres do nieskończoności w ciągu kilku klatek.

Dwie praktyczne poprawki chronią cię bez spowalniania widocznej animacji. Po pierwsze, uruchamiaj kilka małych podkroków fizyki na wyrenderowaną klatkę — matematyka pozostaje stabilna, a ekran wciąż aktualizuje się 60 razy na sekundę. Po drugie, przycinaj lub delikatnie tłum pole przy każdym kroku, by przypadkowy skok nie narastał lawinowo. PDE w stylu fal, jak dyspersyjny wzór fal w naszej symulacji kilwateru Kelvina, są tutaj szczególnie wrażliwe: V-kształtny kilwater utrzymuje swój kształt tylko wtedy, gdy prędkość fali i krok czasowy pozostają w równowadze.

4. Oddziel Fizykę od Renderowania

Powszechnym błędem jest wiązanie kroku całkowania bezpośrednio z requestAnimationFrame i przekazywanie rzeczywistego upłyniętego czasu jako dt. Gdy karta zawiśnie, ta delta rozdyma się do setek milisekund, a pojedynczy gigantyczny krok detonuje twój stabilny schemat. Rozwiązaniem jest akumulator o stałym kroku czasowym: przesuwaj fizykę stałymi wycinkami i renderuj tylko raz na klatkę.

let acc = 0;
const FIXED = 1 / 120; // sekundy fizyki na podkrok

function frame(now) {
  acc += Math.min((now - last) / 1000, 0.1); // ogranicz skok
  last = now;
  while (acc >= FIXED) { step(FIXED); acc -= FIXED; }
  render();
  requestAnimationFrame(frame);
}

Ograniczenie zakumulowanego czasu (Math.min(..., 0.1)) zapobiega osławionej „spirali śmierci", gdzie wolna klatka wymaga więcej podkroków, niż może sobie pozwolić następna klatka. Ten pojedynczy wzorzec to różnica między symulacją, która przetrwa kartę w tle, a taką, która wraca do ekranu pełnego śmieci.

5. Renderuj Szybko: Grupuj Wywołania Canvas

Nawet idealna matematyka wydaje się zepsuta przy 20 FPS. Canvas 2D jest wystarczająco szybki dla tych wizualizacji, jeśli przestaniesz z nim walczyć:

• Rysuj pola za pomocą putImageData. Dla pełnej siatki zapisz bajty RGBA bezpośrednio do bufora ImageData i wyślij całą klatkę jednym wywołaniem — nigdy jeden fillRect na komórkę.
• Grupuj ścieżki cząstek. Dla rojów ODE otwórz pojedynczy beginPath(), wypchnij każde moveTo/lineTo i wydaj jedno stroke(). Każde osobne obrysowanie to kosztowna zmiana stanu.
• Unikaj alokacji na klatkę. Ponownie wykorzystuj swoje tablice typowane i obiekty stanu. Pauzy odśmiecania pamięci są najczęstszą przyczyną widocznego zacinania przy 60 FPS.
• Dopasuj współczynnik pikseli urządzenia raz. Skaluj bufor canvas do devicePixelRatio przy inicjalizacji, nie w każdej klatce, dla ostrych linii na wyświetlaczach retina bez narzutu na każdą klatkę.

Razem zamieniają one powolny prototyp w coś, co wydaje się natychmiastowe, nawet na telefonie ze średniej półki.

Zbierając Wszystko Razem

Przepis jest ten sam, czy ewoluujesz trzy sprzężone ODE, czy pole PDE z milionem komórek: wybierz integrator pasujący do wrażliwości twojego systemu (Euler dla wybaczających, RK4 dla reszty), ułóż swoje pole jako płaską tablicę typowaną z podwójnie buforowaną aktualizacją różnic skończonych, utrzymuj krok czasowy poniżej granicy stabilności, oddziel fizykę od renderowania stałym akumulatorem i grupuj każde wywołanie Canvas, jakie możesz. Opanuj te pięć nawyków, a niemal każde równanie różniczkowe z podręcznika stanie się symulacją, którą możesz wypuścić.

Chcesz zobaczyć te techniki w akcji? Odkryj symulacje atraktora Lorenza, niestabilności plazmy, kilwateru Kelvina i ogrzewania granularnego — każda z nich to równanie różniczkowe uruchamiające opisaną powyżej pętlę. Albo przeglądaj pełny katalog symulacji, by znaleźć równanie stojące za twoim ulubionym zjawiskiem.

Powiązane Wpisy