Devlog #88 – Fala 67: Odwzorowanie logistyczne, wzór mory, triangulacja Delaunaya, korek uliczny, quadtree i pływy

Nowe symulacje

📉 Odwzorowanie logistyczne — okno w uniwersalny chaos

Odwzorowanie i jego bifurkacje

Odwzorowanie logistyczne to najprostszy jednowymiarowy system, który wykazuje pełną drogę podwajania okresu do chaosu:

x_{n+1} = r · x_n · (1 − x_n)      x ∈ [0, 1], r ∈ [0, 4]

Punkt stały (orbita okresu 1):     x* = 1 − 1/r     stabilny dla r < 3
Bifurkacja okresu 2 przy r₁ ≈ 3.000
Bifurkacja okresu 4 przy r₂ ≈ 3.449
Bifurkacja okresu 8 przy r₃ ≈ 3.544
...
Punkt akumulacji (początek chaosu) przy r∞ ≈ 3.56995

Aby wyrenderować diagram bifurkacji, symulacja iteruje odwzorowanie 300 razy (stan przejściowy), a następnie zapisuje kolejne 300 iteracji dla każdej wartości r, rysując pionową kolumnę punktów atraktora. Oś pozioma obejmuje r od 2.5 do 4.0 z rozdzielczością 1200 pikseli. Przeciąganie do powiększania rerenderuje wybrany prostokąt w pełnej rozdzielczości, ujawniając samopodobną strukturę wewnątrz każdego chaotycznego okna.

Stała Feigenbauma

Mitchell Feigenbaum odkrył w 1975 roku, że stosunek kolejnych interwałów bifurkacji zbiega do uniwersalnej stałej, niezależnej od konkretnego odwzorowania:

Stała Feigenbauma:
  δ = lim_{n→∞} (r_n − r_{n-1}) / (r_{n+1} − r_n) ≈ 4.66920160910299...

Zbieżność:
  (r₁ − r₀) / (r₂ − r₁) = (3.449 − 3.000) / (3.544 − 3.449) = 0.449 / 0.095 ≈ 4.73
  (r₂ − r₁) / (r₃ − r₂) ≈ 4.66  (już bardzo blisko granicy)

Ta sama δ pojawia się w odwzorowaniu kwadratowym, odwzorowaniu sinusoidalnym
i każdym unimodalnym odwzorowaniu z kwadratowym maksimum — klasa
uniwersalności systemów chaotycznych.

Symulacja opisuje każdy punkt bifurkacji podwajania okresu pionową przerywaną linią i oblicza bieżący stosunek między kolejnymi interwałami, pokazując zbieżność do δ na małym panelu z boku.

Diagram pajęczyny i wykładnik Lapunowa

Diagram pajęczyny (cobweb plot) to narzędzie geometryczne do wizualizacji dynamiki orbity na paraboli y = rx(1−x). Zaczynając od x₀, narysuj pionową linię do paraboli w punkcie (x₀, x₁), następnie poziomą linię do przekątnej y = x w punkcie (x₁, x₁), potem pionową do paraboli w punkcie (x₁, x₂) i tak dalej. Stabilny punkt stały tworzy spiralę zbiegającą do przecięcia; orbita okresu 2 tworzy prostokąt; chaos tworzy gęstą pajęczynę.

Wykładnik Lapunowa (oszacowanie skończonego czasu na N iteracji):
  λ = (1/N) · Σ_{n=0}^{N-1}  ln |df/dx w x_n|
    = (1/N) · Σ_{n=0}^{N-1}  ln |r(1 − 2x_n)|

  λ < 0 : stabilna orbita okresowa (sąsiednie trajektorie zbiegają)
  λ = 0 : punkt bifurkacji lub graniczna stabilność
  λ > 0 : chaos (sąsiednie trajektorie rozbiegają eksponencjalnie)

Wykres wykładnika Lapunowa jest rysowany bezpośrednio pod diagramem bifurkacji, dzieląc tę samą oś r. Pasma chaotyczne (λ > 0, pokazane na czerwono) i okna okresowe (λ < 0, pokazane na niebiesko) dokładnie pokrywają się ze strukturą bifurkacji powyżej — najbardziej uderzające jest duże okno okresu 3 przy r ≈ 3.83, przewidziane przez twierdzenie Li-Yorke'a „okres 3 implikuje chaos".

🌐 Wzór mory — interferencja, próbkowanie i aliasing

Częstotliwość przestrzenna i okres mory

Siatka liniowa o częstotliwości przestrzennej f (linii na piksel) może być opisana jako okresowa funkcja binarna g(x) = round(sin(2πfx)). Gdy dwie siatki o częstotliwościach f₁ i f₂ są mnożone (lub równoważnie, jedna jest oglądana przez drugą), iloczyn zawiera częstotliwości sumy i różnicy:

g₁(x) · g₂(x) = cos(2πf₁x) · cos(2πf₂x)
              = (1/2)[cos(2π(f₁+f₂)x) + cos(2π(f₁−f₂)x)]

Częstotliwość mory:    f_m = |f₁ − f₂|
Okres mory:       d_m = 1 / f_m    (pikseli na pasmo mory)

Przykład: f₁ = 0.10 linii/px, f₂ = 0.09 linii/px
  f_m = 0.01 linii/px  →  d_m = 100 px (szerokie, widoczne pasmo mory)

Symulacja renderuje każdą siatkę na osobnym płótnie poza ekranem i łączy je za pomocą globalCompositeOperation = 'multiply'. Odczyt na żywo pokazuje f₁, f₂, f_m i d_m aktualizujące się w czasie rzeczywistym podczas przeciągania suwaków częstotliwości.

Obrócone siatki

Gdy dwie siatki o identycznej częstotliwości f są obrócone o kąt θ względem siebie, wzór mory ma okres określony przez obrót:

Okres mory dla siatek o równej częstotliwości pod kątem θ:
  d_m = d / (2 · sin(θ/2))

  Przy θ = 5°:  d_m = d / (2 · sin(2.5°)) ≈ 11.5 d   (szerokie pasma)
  Przy θ = 30°: d_m = d / (2 · sin(15°))  ≈ 1.93 d   (wąskie pasma)
  Przy θ = 0°:  d_m → ∞  (siatki równoległe: brak mory, jednolite pole)

Tarcza obrotu pozwala przemiatać od 0° do 90°. Przy małych kątach pasma mory są szerokie i wolno oscylują; przy 45° otrzymuje się wzór siatki diamentowej; przy 90° dwie siatki są ortogonalne i tworzą kwadratową siatkę. Etykieta okresu mory aktualizuje się w sposób ciągły.

Związek z aliasingiem Nyquista

Mora w obrazach cyfrowych jest bezpośrednim przejawem aliasingu: gdy sensor obrazu próbkuje drobny wzór z częstotliwością poniżej granicy Nyquista, wysokie częstotliwości przestrzenne „zawijają się" do niższych częstotliwości — pojawiając się jako grubsze wzory, których nie było w oryginalnej scenie. Demonstracja aliasingu w symulacji renderuje drobny sinusoidalny cel o częstotliwości f, a następnie próbkuje go w dół z częstotliwością f_s:

Kryterium Nyquista: f_s ≥ 2 · f   (bez aliasingu)

Jeśli f_s < 2f (niedopróbkowanie):
  Częstotliwość aliasu = |f − round(f / f_s) · f_s|

Przykład: f = 0.45 cykli/px, f_s = 1 próbka/px
  Alias = |0.45 − 0| = 0.45  (brak aliasu, ledwo w granicach Nyquista)

  f = 0.55 cykli/px, f_s = 1 próbka/px
  Alias = |0.55 − 1| = 0.45  (alias do 0.45 — to samo co powyżej!)

Wyświetlacz obok pokazuje oryginalny wzór (drobna siatka), wynik próbkowania (mora) i etykietę częstotliwości aliasu. Daje to najbardziej intuicyjne wyjaśnienie, dlaczego aparaty potrzebują optycznych filtrów dolnoprzepustowych przeciw aliasingowi przed sensorami.

△ Triangulacja Delaunaya — optymalne trójkąty i dual Woronoja

Kryterium Delaunaya

Triangulacja zbioru punktów jest triangulacją Delaunaya wtedy i tylko wtedy, gdy dla każdego trójkąta okrąg opisany (circumcircle) nie zawiera żadnego innego punktu zbioru w swoim wnętrzu. To kryterium „pustego okręgu opisanego" jest równoważne maksymalizacji minimalnego kąta wśród wszystkich trójkątów triangulacji — czyniąc triangulacje Delaunaya możliwie najbardziej „równobocznymi" dla danego zbioru punktów.

Okrąg opisany trójkąta (A, B, C):
  Środek D można znaleźć na przecięciu symetralnych AB i BC.

  W praktyce (obliczenia stabilne numerycznie):
    ax, ay = A.x − C.x,  A.y − C.y
    bx, by = B.x − C.x,  B.y − C.y
    D_denom = 2 · (ax·by − ay·bx)
    D.x = (by·(ax²+ay²) − ay·(bx²+by²)) / D_denom + C.x
    D.y = (ax·(bx²+by²) − bx·(ax²+ay²)) / D_denom + C.y
    R   = distance(D, A)

Inspektor przy najechaniu rysuje okrąg opisany trójkąta pod kursorem, zaznaczając go na czerwono, jeśli (błędnie) zawiera inny punkt — przydatny tryb debugowania, który weryfikuje poprawność algorytmu w czasie rzeczywistym.

Inkrementalny algorytm Bowyera-Watsona

Symulacja wykorzystuje algorytm Bowyera-Watsona, który inkrementalnie wstawia punkty jeden po drugim:

Inicjalizacja „supertrójkątem", który zawiera wszystkie punkty.

Dla każdego nowego punktu P:
  1. Znajdź wszystkie „złe" trójkąty: trójkąty, których okrąg opisany zawiera P.
  2. Znajdź wielokąt graniczny „złych" trójkątów:
       krawędzie niewspólne między dwoma złymi trójkątami.
  3. Usuń wszystkie złe trójkąty.
  4. Utwórz nowe trójkąty od P do każdej krawędzi granicznej.
  5. Sprawdź i napraw wszelkie pozostałe naruszenia okręgu opisanego (odwrócenia krawędzi).

Po wstawieniu wszystkich punktów:
  Usuń wszystkie trójkąty mające wspólny wierzchołek z supertrójkątem.

Złożoność: O(n log n) oczekiwana dla jednorodnie losowych zbiorów punktów;
             O(n²) w najgorszym przypadku (konfiguracje zdegenerowane).

Każdy krok wstawiania jest animowany: złe trójkąty migają na czerwono, wielokąt wnęki rysowany jest na pomarańczowo, po czym nowe trójkąty wachlarzowo rozchodzą się od wstawionego punktu. Tryb krok po kroku pozwala przechodzić po jednym punkcie, aby prześledzić algorytm.

Diagram Woronoja jako dual

Diagram Woronoja jest geometrycznym dualem triangulacji Delaunaya. Dla każdego trójkąta Delaunaya umieść wierzchołek w jego środku okręgu opisanego; połącz środki okręgów opisanych sąsiednich trójkątów (mających wspólną krawędź). Wynikiem jest diagram Woronoja: każda komórka zawiera wszystkie punkty bliższe jednemu punktowi wejściowemu niż jakiemukolwiek innemu.

Dualność:
  Krawędź Delaunaya (A,B) ↔ krawędź Woronoja rozdzielająca komórki A i B
  Wierzchołek Delaunaya P   ↔ komórka Woronoja punktu P
  Trójkąt Delaunaya   ↔ wierzchołek Woronoja (środek okręgu opisanego)

Komórka Woronoja punktu P to wielokąt wypukły:
  { x : |x − P| ≤ |x − Q| dla wszystkich innych punktów wejściowych Q }

Przełącz przycisk „Pokaż Woronoja", aby nałożyć graf dualny. Obie struktury są rysowane jednocześnie różnymi kolorami, co czyni relację dualności wizualnie oczywistą: każda krawędź Delaunaya jest przecinana prostopadle przez odpowiadającą jej krawędź Woronoja.

🚗 Korek uliczny — jak fantomowe korki powstają z reguł

Model Nagela-Schreckenberga

Automat komórkowy Nagela-Schreckenberga (NaSch) (1992) modeluje ruch drogowy na jednopasmowej drodze kolistej podzielonej na L komórek. Każda komórka jest pusta lub zajęta przez jeden samochód o całkowitej prędkości v ∈ {0, 1, 2, …, v_max}. Jeden krok czasowy stosuje jednocześnie cztery reguły do wszystkich samochodów:

Niech d = odstęp z przodu (liczba pustych komórek przed następnym samochodem)
    v = obecna prędkość
    v_max = maksymalna prędkość (domyślnie 5)
    p = prawdopodobieństwo losowości (domyślnie 0.3)

Reguła 1 — Przyspieszanie:
    jeśli v < v_max:  v ← v + 1

Reguła 2 — Hamowanie (unikanie kolizji):
    jeśli v ≥ d:  v ← d − 1

Reguła 3 — Losowość (szum ludzkiej reakcji):
    jeśli v > 0 i random() < p:  v ← v − 1

Reguła 4 — Ruch:
    samochód przesuwa się o v komórek do przodu

Krok losowości jest kluczowy: bez niego (p = 0) ruch płynie swobodnie przy każdej gęstości. Przy p > 0 losowe zdarzenia hamowania rozprzestrzeniają się wstecz jako fale „stop-and-go" — fantomowy korek, którego kierowcy doświadczają na autostradach bez widocznej przyczyny.

Diagram czasoprzestrzenny i diagram fundamentalny

Diagram czasoprzestrzenny pokazuje pozycje samochodów (oś pozioma = pozycja na drodze, oś pionowa = czas, przewijana w dół). Każdy samochód to punkt; korki pojawiają się jako ciemne przekątne pasma nachylone w lewo — wskazując, że korek przemieszcza się wstecz z w przybliżeniu stałą prędkością. Prędkość korka w modelu NaSch wynosi zwykle −1 komórkę/krok (jedna komórka wstecz na krok czasowy) niezależnie od gęstości ruchu.

Zmienne przepływu ruchu:
  ρ = gęstość = N / L  (samochodów na komórkę, 0 ≤ ρ ≤ 1)
  q = przepływ    = ρ · ⟨v⟩  (samochodów przechodzących przez punkt na krok czasowy)

Diagram fundamentalny (q względem ρ) dla NaSch z v_max = 5, p = 0.3:
  Gałąź swobodnego przepływu  (ρ < ρ_crit ≈ 0.06):  q ≈ ρ · v_max
  Gałąź zatłoczona  (ρ > ρ_crit):          q maleje wraz z ρ
  Maksymalny przepływ przy  ρ_crit ≈ 0.06,  q_max ≈ 0.32 samochodów/krok/komórkę

Żywy diagram fundamentalny jest rysowany jako punkt rozrzutu (bieżąca gęstość, bieżący przepływ) na płaszczyźnie q-ρ, który kreśli charakterystyczny odwrócony kształt V w miarę zmiany gęstości suwakiem. Histereza między gałęzią swobodnego przepływu a zatłoczoną, znana z empirycznych danych autostradowych, jest jakościowo odtwarzana przez model.

Przejście fazowe i efekt „szybciej znaczy wolniej"

Przy krytycznej gęstości model przechodzi przejście fazowe pierwszego rodzaju od swobodnego przepływu do zatłoczonego. Powyżej ρ_crit zwiększanie v_max może paradoksalnie zmniejszyć przepustowość: większe indywidualne prędkości oznaczają wymagane większe bezpieczne odstępy, zmniejszając gęstość, co zmniejsza przepływ. Ten efekt „szybciej znaczy wolniej" obserwuje się także w modelach tłumów pieszych (Helbing) i wiąże się z szerszą fizyką napędzanych systemów dysypatywnych w pobliżu przejść do zatoru.

🌲 Quadtree — podział przestrzenny i szeroka faza kolizji

Struktura quadtree i wstawianie

Quadtree rekurencyjnie dzieli region 2D na cztery równe kwadranty. Każdy węzeł wewnętrzny reprezentuje prostokąt; każdy liść albo przechowuje jeden punkt, albo jest pusty. Wstawianie nowego punktu P:

function insert(node, P):
    jeśli node jest liściem i node jest pusty:
        zapisz P w node
        zwróć
    jeśli node jest liściem i node zawiera punkt Q:
        podziel node na 4 dzieci
        wstaw ponownie Q do odpowiedniego dziecka
        wstaw P do odpowiedniego dziecka
        zwróć
    // node jest wewnętrzny:
    wstaw P do kwadrantu dziecka zawierającego P

Głębokość dla n jednorodnie losowych punktów: O(log n) oczekiwana
Najgorszy przypadek: n punktów współliniowych → O(n) głębokość (przypadek zdegenerowany)

Wizualizacja płótna rysuje każdą granicę prostokąta w czasie rzeczywistym. W miarę dodawania punktów jeden po drugim, pojawiają się nowe podziały z krótką animacją błysku. Głębokość drzewa i liczba liści są pokazywane na panelu statystyk.

Zapytanie zakresowe z przycinaniem gałęzi

Zapytanie zakresowe pyta: „które punkty leżą wewnątrz tego prostokąta R?" Naiwne przeszukiwanie sprawdza wszystkie n punktów w czasie O(n). Quadtree przycina gałęzie, których prostokąt otaczający nie nachodzi na R:

function range_query(node, R):
    jeśli node.bounds nie nachodzi na R:
        zwróć []    // przytnij: pomiń całe poddrzewo
    jeśli node jest liściem:
        zwróć [P dla P w node, jeśli P w R]
    zwróć range_query(NW) + range_query(NE) +
           range_query(SW) + range_query(SE)

Złożoność: O(√n + k) dla jednorodnie losowych punktów,
            gdzie k = liczba zwróconych wyników.

Interaktywne zapytanie zakresowe pozwala narysować prostokąt zapytania kliknięciem i przeciągnięciem. Przycięte poddrzewa są natychmiast zacieniane na szaro; aktywne poddrzewa świecą. Licznik pokazuje „odwiedzone węzły" wobec „wszystkich węzłów", aby ilościowo określić korzyść z przycinania.

Szeroka faza wykrywania kolizji

Panel wykrywania kolizji symuluje n poruszających się okręgów. Naiwne podejście testuje wszystkie n(n−1)/2 par w każdej klatce. Szeroka faza quadtree odbudowuje drzewo w każdej klatce (O(n log n)) i testuje tylko okręgi w tej samej lub sąsiednich komórkach jako kandydatów szerokiej fazy:

Naiwne:     O(n²) sprawdzeń na klatkę
Quadtree:  O(n log n) budowa + O(n·k) sprawdzeń
           gdzie k = średnia liczba sąsiadów na komórkę ≈ stała

Współczynnik przyspieszenia przy n = 500:
  Naiwne: 124 750 sprawdzeń par na klatkę
  Quadtree: ~2500 kandydatów szerokiej fazy (≈ 50× mniej sprawdzeń)

Licznik na żywo pokazuje obie liczby sprawdzeń obok siebie, aktualizując się co klatkę. Współczynnik przyspieszenia jest wykreślany względem n jako wykres pomocniczy, potwierdzając asymptotyczną rozbieżność O(n²) wobec O(n log n).

🌊 Pływy — gradient grawitacji

Siła pływowa: gradient grawitacji

Siły pływowe nie są samą grawitacją, lecz różnicą w przyspieszeniu grawitacyjnym na przestrzeni ciała o skończonym rozmiarze. Dla układu Ziemia-Księżyc, przyspieszenie pływowe w punkcie przesuniętym o Δr od środka Ziemi w kierunku Księżyca (odległość D) wynosi:

Przyspieszenie grawitacyjne od Księżyca w środku Ziemi:
  a₀ = G·M_Księżyc / D²

W punkcie powierzchni przesuniętym o Δr w kierunku Księżyca:
  a = G·M_Księżyc / (D − Δr)²
  ≈ G·M_Księżyc / D² · (1 + 2Δr/D)   (rozwinięcie Taylora, Δr ≪ D)

Przyspieszenie pływowe (różnicowe):
  Δa = a − a₀ = 2·G·M_Księżyc·Δr / D³   ∝ 1/D³

Ta zależność 1/D³ oznacza, że siły pływowe maleją szybciej niż grawitacja (1/D²),
więc efekt pływowy Księżyca na Ziemię przewyższa efekt Słońca, mimo większego
przyciągania grawitacyjnego Słońca:
  (Siła pływowa)_Księżyc / (Siła pływowa)_Słońce ≈ 2.2

Wyświetlacz pola wektorowego rysuje strzałki przyspieszenia pływowego w punktach siatki wokół powierzchni Ziemi. Strzałki tworzą charakterystyczny wzór „dwóch płatów": rozciąganie na zewnątrz w kierunku i od Księżyca, ściskanie do wewnątrz prostopadle do linii Ziemia-Księżyc. Daje to dwa garby pływowe i wyjaśnia mniej więcej dwukrotny w ciągu doby cykl przypływów.

Pływy syzygijne i kwadraturowe

Księżyc i Słońce każde wytwarzają własną siłę pływową. Ich względne ustawienie określa typ pływu:

Syzygia (ustawienie Księżyc-Ziemia-Słońce lub Słońce-Ziemia-Księżyc):
  Siły pływowe sumują się konstruktywnie → pływy syzygijne
  Występuje przy nowiu (koniunkcja) i pełni (opozycja)
  Amplituda pływu syzygijnego ≈ 1 + 0.46 = 1.46 × średnia

Kwadratura (Księżyc 90° od Słońca widziany z Ziemi):
  Siły pływowe częściowo się znoszą → pływy kwadraturowe
  Występuje w pierwszej i ostatniej kwadrze
  Amplituda pływu kwadraturowego ≈ 1 − 0.46 = 0.54 × średnia

  (gdzie 0.46 to stosunek sił pływowych Słońce/Księżyc wynikający z masy i odległości)

Tarcza fazy Księżyca przechodzi przez 29,5-dniowy miesiąc synodyczny. W miarę jej obrotu superpozycja pola wektorowego aktualizuje się w czasie rzeczywistym, a mareograf poniżej rejestruje wynikowy zakres pływów. Szczyty pływów syzygijnych pokrywają się z pozycjami nowiu i pełni; minima kwadraturowe przypadają na kwadry.

Blokada pływowa

Księżyc zawsze zwraca do Ziemi tę samą stronę — konsekwencja blokady pływowej. Grawitacja Ziemi wywołuje garb pływowy w skale Księżyca. We wczesnej historii Księżyca, gdy obracał się szybciej, garb ten był przesuwany do przodu względem linii Ziemia-Księżyc przez bezwładność obrotową. Moment grawitacyjny od Ziemi na tym niewspółosiowym garbie hamował obrót Księżyca, aż jego okres obrotu zrównał się z jego okresem orbitalnym. Symulacja pokazuje animację tego procesu z oznaczonym kątem garbu pływowego i hamującym momentem obrotowym.

Skala czasowa blokady pływowej:
  t_lock ∝ a⁶ · ω₀ · Q / (M_planeta · R_księżyc⁵)

  gdzie a = wielka półoś orbity, Q = współczynnik jakości pływowej,
        ω₀ = początkowa prędkość obrotowa, R_księżyc = promień ciała zablokowanego pływowo

Księżyc-Ziemia: t_lock ≈ kilkaset milionów lat (już zablokowany)
Ziemia-Słońce:  t_lock ≈ 50 miliardów lat (nie zablokuje się w czasie życia Słońca)

Co dalej

Fala 68 kontynuuje rozpęd kolejnymi sześcioma symulacjami: grawitacja N ciał Barnesa-Huta O(n log n), konturowanie izoliniowe maszerującymi kwadratami z metakulami, próbkowanie niebiesko-szumowe dysków Poissona metodą Bridsona, eksperyment nierówności Bella CHSH osiągający granicę Tsirelsona, atraktor podwójnej spirali obwodu Chuy oraz model siły społecznej Helbinga tłumu pieszych z formowaniem pasów ruchu. Ta fala i jej devlog są już dostępne — zobacz Devlog #89.