Devlog #93 – Fala 72: Filtr Kalmana, Propagacja Wsteczna, Wahadło 3D, Eutrofizacja, Soczewka Akustyczna i Pierścienie Planetarne

Nowe symulacje

Filtr Kalmana: optymalna estymacja stanu

Filtr Kalmana (1960) to jeden z najbardziej udanych algorytmów matematyki stosowanej. Optymalnie łączy przewidywanie modelu dynamicznego z zaszumionymi pomiarami z czujników, tworząc nieobciążony estymator stanu o minimalnej wariancji — i, co kluczowe, dokładnie śledzi, jak bardzo jest niepewny, za pomocą macierzy kowariancji P.

Stan 4D i pełne operacje macierzowe

Wektor stanu to [x, y, v_x, v_y] — pozycja i prędkość w 2D. Macierz kowariancji P 4×4 śledzi niepewności i korelacje między wszystkimi składowymi stanu. Każdy krok predykcji oblicza:

x̂⁻ = F·x̂          (predykcja stanu przez macierz przejścia F)
P⁻  = F·P·Fᵀ + Q   (predykcja kowariancji, Q = szum procesu)

Każdy krok aktualizacji (przy otrzymaniu obserwacji pozycji 2D przez H = [I₂|0₂]):

K  = P⁻·Hᵀ·(H·P⁻·Hᵀ + R)⁻¹  (wzmocnienie Kalmana)
x̂  = x̂⁻ + K·(z − H·x̂⁻)       (aktualizacja stanu)
P  = (I − K·H)·P⁻·(I − K·H)ᵀ + K·R·Kᵀ  (forma Josepha)

Forma Josepha dla stabilności numerycznej

Standardowa aktualizacja P = (I−KH)P⁻ cierpi na asymetrię wynikającą z arytmetyki zmiennoprzecinkowej — po wielu iteracjach P może utracić dodatnią określoność. Forma Josepha P = (I−KH)P⁻(I−KH)ᵀ + KRKᵀ jest algebraicznie równoważna, ale numerycznie stabilna: gwarantuje, że P pozostaje symetryczna i dodatnio określona niezależnie od zaokrągleń zmiennoprzecinkowych.

Wizualizacja elipsy kowariancji

Niepewność pozycyjna 2D jest wizualizowana jako elipsa obliczona z rozkładu na wartości własne podmacierzy pozycyjnej 2×2 macierzy P. Półosie to √(λ₁) i √(λ₂) (odchylenia standardowe), a orientacja podąża za wektorami własnymi. Elipsa wizualnie puchnie podczas kroku predykcji i kurczy się, gdy nadchodzi obserwacja — piękna demonstracja cyklu informacyjnego filtra.

• 4 trajektorie: linia prosta, orbita kołowa, błądzenie losowe, manewr wymijający
• Porównanie RMSE: estymata Kalmana kontra surowe pomiary kontra sama predykcja modelu
• Suwaki szumu procesu Q i szumu pomiaru R
• Przełącznik częstotliwości obserwacji, by zobaczyć wzrost kowariancji w przerwach

Propagacja wsteczna: przepływ gradientu uwidoczniony

Propagacja wsteczna to po prostu reguła łańcuchowa zastosowana efektywnie do grafu obliczeniowego — ale zrozumienie jej na poziomie intuicyjnym wymaga zobaczenia, jak impulsy gradientu przepływają wstecz przez każdą warstwę. Ta symulacja czyni ten przepływ jawnym.

Ręczna implementacja reguły łańcuchowej

Architektura sieci jest konfigurowalna: 2 wejścia → H jednostek ukrytych → H jednostek ukrytych → 2 (klasyfikacja) lub 1 (regresja) wyjścia. Każdy przebieg w przód oblicza aktywacje warstwa po warstwie. Każdy przebieg wsteczny oblicza błędy δ: warstwa wyjściowa otrzymuje δ_L = ∂L/∂z_L, następnie propaguje: δ_l = (Wᵀ_{l+1} · δ_{l+1}) ⊙ σ'(z_l). Gradienty wag to ∂L/∂W_l = δ_l · aᵀ_{l-1}.

Wizualna animacja gradientu

Podczas propagacji wstecznej kolorowe impulsy przemieszczają się od prawej do lewej przez diagram sieci. Grubość i intensywność koloru krawędzi reprezentują |∂L/∂w| dla każdej wagi. Duże gradienty świecą jasno; bliskie zeru (zanikający gradient) pokazane są jako przygaszone linie. Przy ReLU martwe neurony (zawsze zwracające 0) są jawnie oznaczone — ich wychodzące gradienty są trwale zerowe.

• Funkcje aktywacji: ReLU, tanh, sigmoid; wyjście: softmax+CE lub liniowe+MSE
• Mapa cieplna granicy decyzyjnej aktualizowana co N kroków treningu
• Wykres paskowy straty: strata treningowa w czasie
• Tryb krok po kroku dla pojedynczej próbki do nauki algorytmu

Wahadło 3D: więz mnożnika Lagrange'a

Wahadło sferyczne — masa na sztywnym pręcie, swobodnie wychylająca się we wszystkich kierunkach — jest łatwe do opisania fizycznie, ale niewygodne do symulacji numerycznej. Więzem jest to, że odległość od zawieszenia do masy jest zawsze dokładnie równa L. Zamiast przechodzić na współrzędne sferyczne (które wprowadzają osobliwości gimbal lock), symulacja działa we współrzędnych kartezjańskich i wymusza więz za pomocą mnożnika Lagrange'a.

RK4 kartezjańskie z projekcją więzu

Równania ruchu we współrzędnych kartezjańskich (x,y,z) to: m·r̈ = F_grawitacja + λ·r̂, gdzie λ to wielkość siły więzu (naprężenie). Na każdym kroku RK4, po swobodnym przesunięciu r i v, pozycje są rzutowane z powrotem na sferę, a prędkości na płaszczyznę styczną: λ = (−|v|² − r·g_wek) / L², następnie stosowana jest korekta. To zapewnia dokładne spełnienie więzu bez osobliwości układu współrzędnych.

Trzy tryby fizyczne

• Płaska rozeta: warunki początkowe w jednej płaszczyźnie — tworzy wzór przypominający figurę Lissajous widziany z góry z powodu nieliniowego okresu
• Wahadło stożkowe: masa krąży pod stałym kątem θ z okresem T = 2π√(L cosθ / g) — krótszym niż okres wahadła prostego z powodu zredukowanej efektywnej grawitacji
• Wahadło Foucaulta: dodaje siłę Coriolisa F_c = −2mΩ × v z Ω = Ω_Ziemia·sin(szerokość geograficzna)·ẑ, powodując precesję płaszczyzny oscylacji z prędkością Ω·sin(szer.)

• Projekcja 3D z obrotem przeciąganiem myszy
• Suwak szerokości geograficznej dla trybu Foucaulta (0° = równik, 90° = biegun)
• Długość śladu i gradient koloru w czasie

Eutrofizacja: punkt krytyczny jeziora

Eutrofizacja — wzbogacenie w składniki odżywcze prowadzące do zakwitu glonów i deficytu tlenu — jest jednym z najważniejszych przykładów katastroficznego punktu krytycznego w naukach środowiskowych. Gdy jezioro przekroczy próg, przeskakuje do mętnego, ubogiego w tlen stanu, który może utrzymywać się przez dekady, nawet po zmniejszeniu dopływu fosforu.

Układ pięciu sprzężonych równań różniczkowych

Symulacja całkuje pięć sprzężonych równań różniczkowych reprezentujących P (rozpuszczony fosfor), A (biomasa glonów), O (rozpuszczony tlen), F (biomasa ryb) i D (detrytus). Kluczową nieliniowością jest człon recyklingu fosforu: w warunkach beztlenowych (O poniżej progu) osad uwalnia zmagazynowany fosfor — dodatnie sprzężenie zwrotne, które popycha system do alternatywnego stanu stabilnego. Integrator RK4 niezawodnie radzi sobie ze sztywnymi przejściami.

Histereza i wykres fazowy

Wykres fazowy śledzi (ładunek fosforu, biomasa glonów) w czasie. Gdy ładunek powoli rośnie, system pozostaje w stanie czystej wody, aż nagle przeskakuje do stanu zakwitu w górnym punkcie bifurkacji siodło-węzeł. Gdy ładunek jest następnie zmniejszany, system pozostaje w stanie zakwitu aż do dolnego punktu bifurkacji siodło-węzeł — to jest histereza. Pole zamknięte przez pętlę histerezy na wykresie fazowym reprezentuje „dług odbudowy" — ile dodatkowego wysiłku potrzeba, by przywrócić jezioro w porównaniu z zapobieżeniem zakwitowi od samego początku.

• Wizualizacja zbiornika: gradient koloru wody od czystego błękitu do mętnej zieleni, mata powierzchniowa zakwitu, ikony ryb, ciemność osadu
• Suwak ładunku fosforu z trybami powolnego wzrostu i powolnego spadku
• Odczyty poziomu tlenu i populacji ryb w czasie rzeczywistym

Soczewka akustyczna: symulacja fal metodą FDTD

Ogniskowanie akustyczne — koncentrowanie fal dźwiękowych w punkcie — leży u podstaw ultrasonografii medycznej, sonaru i akustyki architektonicznej. Ta symulacja rozwiązuje równanie fali 2D bezpośrednio metodą różnic skończonych w dziedzinie czasu (FDTD), pokazując zaginanie się czół fal w czasie rzeczywistym.

Leapfrog FDTD i granice Mura

Równanie fali ciśnienia 2D ∂²p/∂t² = c²(∂²p/∂x² + ∂²p/∂y²) jest dyskretyzowane na siatce 240×144 metodą leapfrog: ciśnienie w chwili n+1 zależy od chwil n i n-1 oraz przestrzennego laplasjanu w chwili n. Warunek stabilności CFL wymaga c·dt/dx ≤ 1/√2. Warunek brzegowy absorbujący Mura pierwszego rzędu eliminuje odbicia od krawędzi siatki, dopasowując prędkość fali wychodzącej: p_brzeg(t+dt) = p_{brzeg+1}(t) + (c·dt−dx)/(c·dt+dx)·(p_{brzeg+1}(t+dt) − p_brzeg(t)).

Cztery typy soczewek

• Soczewka wypukła: wyższa prędkość dźwięku w obszarze soczewki (materiał o niższej gęstości) skupia czoło fali
• Soczewka wklęsła: ten sam materiał, ale wklęsła geometria rozprasza wiązkę
• Soczewka Fresnela: konstrukcja płyty strefowej — naprzemienne pierścienie koncentryczne o wysokiej/niskiej prędkości, które sumują się w fazie w punkcie ogniskowym
• GRIN (paraboliczna): prędkość zmienia się jako c(r) = c₀ / (1 + (r/r₀)²) — soczewka o gradientowym współczynniku bez dyskretnych granic

• Rozbieżna mapa cieplna: niebiesko-czerwone pole ciśnienia z addytywnym poświatą
• Nakładka promieni wg prawa Snella: przybliżenie optyki geometrycznej nałożone na symulację falową
• Ciśnienie EMA w przewidywanym punkcie ogniskowym weryfikuje długość ogniskowej wobec wzoru paraosiowego

Pierścienie planetarne: cząstki keplerowskie i luki rezonansowe

Pierścienie Saturna są jedną z najbardziej fotogenicznych struktur Układu Słonecznego — a ich struktura jest kontrolowana przez mechanikę orbitalną. Podział Cassiniego, wyraźna luka między pierścieniami A i B, jest opróżniana przez rezonans orbitalny 2:1 z księżycem Mimas.

Cząstki keplerowskie i kopnięcia rezonansowe

Do 10 000 cząstek testowych umieszczonych jest na orbitach keplerowskich z półosią wielką a rozłożoną jak n ∝ a^(−3/2) (płaska gęstość powierzchniowa). Każda cząstka otrzymuje kopnięcie grawitacyjne od księżyca podobnego do Mimasa przy każdej koniunkcji. W promieniu rezonansu 2:1 (gdzie cząstka wykonuje dokładnie 2 orbity na jedną orbitę księżyca) kopnięcia te zawsze zachodzą w tej samej fazie orbitalnej — kumulują się spójnie i podbijają mimośród, aż orbita zaczyna się przecinać, po czym cząstka jest szybko rozpraszana. Po tysiącach orbit promień rezonansu wykazuje lukę gęstości odpowiadającą podziałowi Cassiniego.

Granica Roche'a i księżyce-pasterze

Granica Roche'a R_Roche = a_planeta · (2ρ_planeta/ρ_satelita)^(1/3) wyznacza miejsce, w którym siły pływowe przewyższają samograwitację — księżyce wewnątrz niej są rozrywane na materiał pierścieni. Symulacja obsługuje także pary księżyców-pasterzy: dwa małe księżyce okalające wąski pierścionek wymieniają moment pędu z cząstkami pierścienia, utrzymując pierścionek ściśle ograniczony — mechanizm stojący za strukturą pierścienia F Saturna.

• Suwak N: od 100 do 10 000 cząstek
• Regulowana masa i promień orbity księżyca — obserwuj powstawanie i przesuwanie się luki
• Kolor według mimośrodu, by wizualizować wzbudzenie rezonansowe
• Przełącznik pary księżyców-pasterzy dla ograniczenia wąskiego pierścionka

Co dalej

Fala 73 to specjalne wydanie WebGL/GLSL: wszystkie pięć symulacji wykorzystuje własne shadery, raymarching lub techniki GPGPU. Nadchodzi: raymarching fraktala 3D Mandelbulb z estymatorem odległości, gąbka Mengera SDF z wymiarem Hausdorffa, osadzenie paraboloidy Flamma tunelu czasoprzestrzennego Morrisa-Thorne'a, produkt sferyczny supershape Gielisa oraz symulacja GPGPU z 250 000 cząstek wykorzystująca cele renderowania ping-pong.