Devlog #87 – Fala 66: Epicykle Fouriera, wahadło magnetyczne, zwilżanie powierzchni, przekładnia zębata, plecak & Huffman

Nowe symulacje

🌀 Epicykle Fouriera — rysowanie za pomocą DFT

Od narysowanej ścieżki do złożonej DFT

Gdy rysujesz zamkniętą krzywą na płótnie, symulator próbkuje ją przy N równomiernie rozłożonych wartościach parametru, traktując każdy punkt próbki jako liczbę zespoloną z_n = x_n + i·y_n. Obliczana jest złożona Dyskretna Transformata Fouriera:

X[k] = Σ_{n=0}^{N-1}  z_n · exp(−2πi·k·n / N)    dla k = 0, 1, …, N−1

Każdy koszyk k reprezentuje jeden epicykl:
  częstotliwość środkowa : k  obrotów na pełny okres rysowania
  amplituda              : |X[k]| / N
  faza początkowa        : arg(X[k])

Rekonstrukcja (suma częściowa M epicykli):
  z(t) = (1/N) · Σ_{k=0}^{M-1}  X[k] · exp(2πi·k·t / N)

Symulacja sortuje epicykle malejąco wg amplitudy, aby największe okręgi rysowały się jako pierwsze. Gdy M rośnie od 1 do N/2, ścieżka sumy częściowej zbiega się do oryginalnego rysunku. Suwak pozwala obserwować poprawę przybliżenia okrąg po okręgu.

Czym różni się od symulatora Szeregu Fouriera

Istniejący na platformie symulator fourier-series stosuje 1D DFT do okresowych przebiegów (kwadratowy, piłokształtny, trójkątny) i wizualizuje amplitudy harmonicznych na wykresie słupkowym. Symulator fourier-epicycles używa 2D złożonej DFT na dowolnej narysowanej ścieżce — zarówno x(t), jak i y(t) są kodowane razem jako pojedynczy sygnał zespolony. Oznacza to, że każdy koszyk częstotliwości wnosi 2D ruch kołowy (epicykl), a nie tylko oscylację pionową. Zbiór obracających się okręgów jest wizualnym dowodem, że dowolna ciągła zamknięta krzywa może zostać rozłożona na ruchy kołowe.

Wbudowane presety i ich widma

Pięć presetów zostało wybranych, aby zilustrować różne profile widmowe:

• Koło — pojedynczy niezerowy koszyk przy k = 1; jeden epicykl wystarcza do idealnej rekonstrukcji.
• Gwiazda (5-ramienna) — dominujące harmoniczne przy k = 1, 4, 6, 9, 11, … (wielokrotności 5 ± 1); daje charakterystyczne pięcioramienne rozbłyski.
• Serce — asymetryczny kształt wymaga wielu harmonicznych; powolna zbieżność ujawnia ostrze na dole.
• Ósemka — symetryczna; w składowej x obecne są tylko nieparzyste harmoniczne, dając czyste widmo.
• Lemniskata (∞) — sparametryzowana jako x = cos t / (1 + sin²t), y = sin t cos t / (1 + sin²t); zdominowana przez k = 2 z powodu podwójnej pętli.

Renderowanie i wydajność

Obliczenie pełnej DFT przy N = 512 kosztuje O(N²) mnożeń. Symulator używa iteracyjnego algorytmu FFT Cooleya-Tukeya radix-2 (dopełniając do najbliższej potęgi dwójki), aby zredukować to do O(N log N):

// Iteracyjna FFT Cooleya-Tukeya w miejscu
function fft(x):   // x: tablica zespolona, długość musi być potęgą 2
    N = x.length
    permutacja bit-reversal x
    for s = 1 to log2(N):
        m = 2^s
        w_m = exp(-2πi / m)
        for k = 0 to N-1 step m:
            w = 1
            for j = 0 to m/2 - 1:
                t = w * x[k + j + m/2]
                u = x[k + j]
                x[k + j]       = u + t
                x[k + j + m/2] = u - t
                w = w * w_m

Animacja epicykli renderuje się przy 60 fps, przeliczając co klatkę tylko pozycje końców obracających się ramion (bez FFT na klatkę — DFT jest obliczana raz po puszczeniu przycisku myszy). Ślad śledzący jest rysowany na płótnie poza ekranem i komponowany, dzięki czemu główna pętla pozostaje lekka.

🧲 Wahadło magnetyczne — chaos i fraktalne baseny

Równania ruchu

Kulka wahadła jest modelowana jako masa punktowa m zawieszona na sznurku o długości L, ograniczona do ruchu w płaszczyźnie poziomej bezpośrednio nad trzema magnesami. Wysokość z kulki nad płaszczyzną magnesów jest utrzymywana na stałym poziomie h (wysokość spoczynkowa). Równania ruchu we współrzędnych kartezjańskich (x, y) to:

Siły działające na kulkę:
  Grawitacja przywracająca:  F_g = −(mg/L) · (x, y)   (przybliżenie małych kątów)
  Tłumienie:                 F_d = −b · (ẋ, ẏ)
  Magnetyczna (magnes i w pozycji r_i = (x_i, y_i)):
    d_i = sqrt((x − x_i)² + (y − y_i)² + h²)
    F_m_i = μ · (r_i − r) / d_i³

Pełne ODE:
  ẍ = −(g/L)·x  −  b·ẋ  +  Σ_i  μ·(x_i − x) / d_i³
  ÿ = −(g/L)·y  −  b·ẏ  +  Σ_i  μ·(y_i − y) / d_i³

Parametry domyślne: L = 0,5 m, h = 0,05 m, b = 0,2 (tłumienie), μ = 1,0 (siła magnetyczna). Trzy magnesy leżą pod kątami 0°, 120°, 240° na okręgu o promieniu 0,3 m.

Integracja RK4

Każda trajektoria jest całkowana klasyczną metodą Rungego-Kutty 4. rzędu ze stałym krokiem czasowym Δt = 0,01 s. Wektor stanu to s = (x, y, ẋ, ẏ):

k1 = f(t,       s)
k2 = f(t + h/2, s + h·k1/2)
k3 = f(t + h/2, s + h·k2/2)
k4 = f(t + h,   s + h·k3)
s_{n+1} = s_n + (h/6)·(k1 + 2k2 + 2k3 + k4)

Symulacja zatrzymuje się, gdy prędkość kulki spada poniżej 0,001 m/s przez 50 kolejnych kroków, rejestrując najbliższy magnes jako „przechwytujący". Dla mapy basenów siatka pozycji startowych (do 400×400 = 160 000 trajektorii) jest obliczana przy użyciu Web Workerów dla podzielonego na fragmenty renderowania, aby interfejs pozostał responsywny.

Fraktalne granice basenów

Wynikowa mapa basenów jest kanonicznym przykładem fraktala. Na każdej pozornej granicy między dwoma obszarami kolorystycznymi znajduje się trzeci region zagnieżdżony nieskończenie głęboko — granice mają wymiar Hausdorffa ściśle większy niż 1. Jest to bezpośrednia konsekwencja wrażliwości na warunki początkowe (dodatnie wykładniki Lapunowa) w pobliżu separatrys. Powiększenie granicy ujawnia samopodobną drobną strukturę na każdej skali, co czyni to jedną z najbardziej wizualnie przystępnych demonstracji chaosu na platformie.

Dostępne są dwa interaktywne suwaki: „siła magnetyczna", który osłabia lub wzmacnia magnesy względem grawitacji, oraz „tłumienie". Przy bardzo niskim tłumieniu kulka wędruje znacznie dłużej, zanim się ustabilizuje, dając ostrzejszą i głębszą strukturę fraktalną. Przy bardzo wysokim tłumieniu prawie od razu wpada do najbliższego magnesu, upraszczając granice basenów do prostych linii Woronoja.

💧 Zwilżanie powierzchni — od równania Younga do efektu lotosu

Równanie Younga (idealna płaska powierzchnia)

Gdy kropla cieczy spoczywa na płaskiej powierzchni stałej w środowisku pary, na linii styku działają trzy napięcia międzyfazowe: ciało stałe-para (γ_SV), ciało stałe-ciecz (γ_SL) i ciecz-para (γ_LV). Równowaga mechaniczna wzdłuż powierzchni wymaga:

Równanie Younga:
  γ_SV = γ_SL + γ_LV · cos(θ_Y)

  → cos(θ_Y) = (γ_SV − γ_SL) / γ_LV

Zakresy kąta zwilżania:
  θ_Y = 0°         : całkowite zwilżanie (ciecz rozprzestrzenia się bez końca)
  0° < θ_Y < 90°   : hydrofilowa (woda zwilża powierzchnię)
  90° < θ_Y < 150° : hydrofobowa (woda tworzy kropelki)
  θ_Y > 150°        : superhydrofobowa (efekt lotosu)

Animowany przekrój kropli zmienia się w sposób ciągły wraz z suwakiem kąta zwilżania. Trzy strzałki interfejsu aktualizują kierunek i etykietę, pokazując równowagę sił.

Model Wenzla (powierzchnia chropowata)

Prawdziwe powierzchnie nigdy nie są idealnie płaskie. Wenzel (1936) wprowadził współczynnik chropowatości r (stosunek rzeczywistej powierzchni do powierzchni rzutowanej, zawsze r ≥ 1) do bilansu energetycznego:

Równanie Wenzla:
  cos(θ_W) = r · cos(θ_Y)

  Konsekwencja:
    Jeśli θ_Y < 90°  (hydrofilowa): chropowatość zwiększa zwilżanie  → θ_W < θ_Y
    Jeśli θ_Y > 90°  (hydrofobowa): chropowatość zmniejsza zwilżanie → θ_W > θ_Y

Przykłady współczynnika chropowatości:
  Polerowany metal:  r ≈ 1,05
  Piaskowany:        r ≈ 1,5 – 2,0
  Liść lotosu:       r ≈ 2,5 (podwójna struktura mikro- i nanoskali)

Wizualizacja chropowatości powierzchni pokazuje powiększony przekrój podłoża z regulowaną amplitudą i częstotliwością przestrzenną nierówności. Kąt zwilżania kropli aktualizuje się w czasie rzeczywistym wraz ze zmianą zarówno θ_Y, jak i r.

Model Cassiego-Baxtera i efekt lotosu

Gdy powierzchnia ma uwięzione kieszenie powietrza (interfejs kompozytowy), kropla spoczywa częściowo na ciele stałym, a częściowo na powietrzu. Cassie i Baxter (1944) wyprowadzili kąt zwilżania dla tej heterogenicznej powierzchni:

Równanie Cassiego-Baxtera:
  cos(θ_CB) = f₁ · cos(θ₁) + f₂ · cos(θ₂)

Dla kompozytu ciało stałe-powietrze (f₂ = frakcja powietrza, θ₂ = 180°):
  cos(θ_CB) = f_s · cos(θ_Y) − (1 − f_s)

  gdzie f_s = frakcyjny obszar styku z ciałem stałym

Liść lotosu: f_s ≈ 0,05 (tylko 5% styku z ciałem stałym)
  cos(θ_CB) ≈ 0,05 · cos(105°) − 0,95 ≈ −0,963
  θ_CB ≈ 164°  (superhydrofobowa)

Symulator udostępnia przełącznik między stanem Wenzla (kropla penetruje rowki chropowatości) a stanem Cassiego-Baxtera (kropla siedzi na szczytach nierówności z powietrzem pod spodem). Przejście z Cassiego do Wenzla (przejście zwilżające) można wywołać, wciskając kroplę, modelowane jako przekroczenie bariery energetycznej.

⚙️ Przekładnia zębata — przełożenia, moment obrotowy i stopnie złożone

Podstawy przełożenia przekładni

Dwa zazębiające się zewnętrzne koła zębate walcowe o liczbie zębów N₁ i N₂ obracają się z prędkościami kątowymi ω₁ i ω₂. Ponieważ ich koła podziałowe toczą się bez poślizgu, prędkość styczna w punkcie zazębienia jest równa:

Promienie kół podziałowych: r₁ = N₁·m/2,  r₂ = N₂·m/2  (m = moduł = rozmiar zęba)

Warunek zazębienia bez poślizgu: r₁·ω₁ = r₂·ω₂

  → przełożenie: ω₂/ω₁ = N₁/N₂  (wartość bezwzględna)

Kierunek: zazębienie zewnętrzne odwraca obrót → ω₂/ω₁ = −N₁/N₂
          zazębienie wewnętrzne (koło koronowe) zachowuje obrót

Przełożenie momentu obrotowego (idealne, bez tarcia):
  T₂/T₁ = N₂/N₁   (moment obrotowy mnożony o ten sam czynnik, przez który dzielona jest prędkość)

Koła pośrednie

Koło pośrednie (zwane też kołem pasożytniczym) jest wstawiane między koło napędzające a napędzane. Zazębia się z obydwoma, ale nie jest połączone z żadnym obciążeniem ani źródłem mocy. Ponieważ wprowadza dwa odwrócenia zazębienia, kierunek netto koła napędzanego jest taki sam jak koła napędzającego — jakby koła pośredniego wcale nie było. Ogólne przełożenie między pierwszym a ostatnim kołem nie jest zależne od koła pośredniego:

Napędzające (N₁) → Pośrednie (N_i) → Napędzane (N₂):
  ω_i / ω₁ = −N₁/N_i       (pierwsze zazębienie, odwrócone)
  ω₂ / ω_i = −N_i/N₂       (drugie zazębienie, ponownie odwrócone)
  ω₂ / ω₁  = N₁/N₂         (liczba zębów koła pośredniego się redukuje)

Praktyczne zastosowanie: odwrócenie kierunku wału wyjściowego bez zmiany przełożenia

Złożone przekładnie zębate

Złożona przekładnia zębata ma dwa lub więcej kół dzielących wał (sprzężonych, obracających się z tą samą prędkością kątową). Całkowite przełożenie jest iloczynem przełożeń poszczególnych stopni:

3-stopniowa złożona przekładnia:
  Stopień 1: napędzające N₁ zazębia się z kołem N₂ (na wale A)
  Stopień 2: koło N₃ (na wale A) zazębia się z kołem N₄ (na wale B)
  Stopień 3: koło N₅ (na wale B) zazębia się z wyjściowym N₆

  Przełożenie całkowite: ω₆/ω₁ = (N₁·N₃·N₅) / (N₂·N₄·N₆)

Przykład: N₁=10, N₂=40, N₃=12, N₄=60, N₅=15, N₆=75
  Przełożenie = (10·12·15)/(40·60·75) = 1800/180000 = 1/100
  Prędkość wejściowa 1000 obr/min → wyjściowa 10 obr/min, moment ×100

Animacja renderuje każde koło z uproszczonym profilem ewolwentowym zębów. Zazębienie zębów jest śledzone przez monitorowanie punktu styku wzdłuż linii działania, a zęby animują się płynnie, wchodząc i wychodząc z zazębienia. Panel wydajności pokazuje przybliżoną stratę mocy z tarcia zazębienia (typowo 1–3% na stopień dla dobrze naoliwionych kół walcowych).

🎒 Plecak — odsłonięta tabela DP 0/1

Definicja problemu i złożoność brute-force

Mając n przedmiotów, każdy o wadze w_i i wartości v_i, oraz plecak o pojemności W, znajdź podzbiór przedmiotów o maksymalnej łącznej wartości, przy sumie wag ≤ W. Każdy przedmiot może być użyty co najwyżej raz (ograniczenie „0/1"). Podejście brute-force próbuje wszystkich 2ⁿ podzbiorów — wykonalne tylko dla n ≤ ok. 20.

Programowanie dynamiczne od dołu do góry

Rekurencja DP wypełnia tabelę (n+1) × (W+1), gdzie K[i][w] = maksymalna wartość osiągalna przy użyciu pierwszych i przedmiotów w plecaku o pojemności w:

Przypadek bazowy:
  K[0][w] = 0   dla każdego w   (brak przedmiotów → zerowa wartość)

Rekurencja (dla i = 1…n, w = 0…W):
  jeśli w_i > w:
      K[i][w] = K[i-1][w]            // przedmiot i za ciężki: pomiń go
  w przeciwnym razie:
      K[i][w] = max(
          K[i-1][w],                 // nie bierz przedmiotu i
          v_i + K[i-1][w - w_i]     // weź przedmiot i
      )

Odpowiedź: K[n][W]
Złożoność: czas O(n·W), pamięć O(n·W) (redukowalna do O(W) z tablicą przetaczaną)

Symulacja renderuje pełną tabelę, kolorowaną wg wielkości wartości. Animacja podczas wypełniania każdej komórki podświetla dwie porównywane komórki (K[i-1][w] i K[i-1][w−w_i]) krótkim błyskiem, dzięki czemu struktura rekurencji jest natychmiast widoczna.

Odtwarzanie wyboru przedmiotów

Po ukończeniu tabeli faktycznie wybrane przedmioty są odzyskiwane przez odtwarzanie wstecz od K[n][W]:

i = n, w = W
while i > 0 and w > 0:
    if K[i][w] != K[i-1][w]:
        include item i           // komórka się zmieniła → przedmiot i był wzięty
        w = w − w_i
    i = i − 1

Wybrane przedmioty są podświetlane złotą obwódką na liście przedmiotów. Panel porównania obok uruchamia algorytm zachłanny (sortowanie wg stosunku wartość/waga, wypełnianie zachłanne) i oznacza wszelkie rozbieżności — dla plecaka 0/1 podejście zachłanne często nie osiąga optimum.

Dlaczego DP nie jest naprawdę wielomianowe

Złożoność O(nW) jest czasem nazywana „pseudo-wielomianową", ponieważ W może być wykładniczo duże w liczbie bitów użytych do jego zakodowania. Jeśli W = 2^k, tabela ma n·2^k komórek — wykładniczo wobec rozmiaru wejścia. Prawdziwie wielomianowe algorytmy dla plecaka 0/1 są nieznane (jest to NP-trudne). Symulacja czyni to konkretnym za pomocą wyświetlacza „długości bitowej": przesuń W od 100 do 100 000 i obserwuj, jak rozmiar tabeli (i czas renderowania) skaluje się liniowo z W, ale wykładniczo z log W.

🌳 Kodowanie Huffmana — optymalna kompresja bezprefiksowa

Budowa drzewa za pomocą kopca min.

Mając alfabet symboli o znanych częstotliwościach (lub prawdopodobieństwach), algorytm Huffmana buduje optymalny kod binarny bezprefiksowy:

Wejście: symbole s_1…s_k o częstotliwościach f_1…f_k

Algorytm:
1. Wstaw wszystkie symbole do kolejki priorytetowej min. (kopiec min.), kluczowanej wg częstotliwości.
2. Dopóki rozmiar kopca > 1:
     a. Wyciągnij dwa węzły o minimalnej częstotliwości: lewy (f_L), prawy (f_R)
     b. Utwórz węzeł wewnętrzny o częstotliwości f_L + f_R
     c. Wstaw węzeł wewnętrzny z powrotem do kopca
3. Pozostały węzeł to korzeń drzewa Huffmana.

Przypisanie kodu:
  Przejdź drzewo od korzenia: w lewo dopisuje '0', w prawo dopisuje '1'.
  Kod każdego liścia = ciąg bitów na ścieżce od korzenia do tego liścia.

Symulacja animuje każdy krok wyciągania i łączenia. Rosnące drzewo jest rysowane z etykietami częstotliwości na każdym węźle, a dwa wyciągnięte węzły świecą podczas łączenia.

Entropia i średnia długość kodu

Twierdzenie Shannona o kodowaniu źródłowym mówi, że żaden jednoznacznie dekodowalny kod nie może osiągnąć średniej długości kodu krótszej niż entropia H źródła:

Entropia Shannona:
  H = −Σ_{i=1}^{k}  p_i · log₂(p_i)    (bity na symbol)

Średnia długość kodu Huffmana:
  L̄ = Σ_{i=1}^{k}  p_i · len(code_i)

Gwarancja:  H ≤ L̄ < H + 1

Dla optymalnego kodu blokowego o długości n symboli: L̄ → H gdy n → ∞

Wykres na żywo w panelu po prawej stronie przedstawia H (linia odniesienia pozioma) i L̄ (słupek) obok siebie. Podczas edycji tekstu wejściowego obie wartości się aktualizują, a różnica między nimi jest podpisana w milibitach. Dla naturalnego tekstu angielskiego H ≈ 4,1 bitów/symbol, a Huffman osiąga L̄ ≈ 4,2 bitów/symbol — w granicach ok. 2% teoretycznego minimum przy użyciu kodów jednosymbolowych.

Dowód własności bezprefiksowości

Każdy kod uzyskany przez odczyt ścieżek drzewa jest automatycznie bezprefiksowy: żadne słowo kodowe nie jest prefiksem innego, ponieważ słowa kodowe pojawiają się tylko w liściach, a żaden liść nie jest przodkiem innego liścia. Gwarantuje to jednoznaczne dekodowanie bez żadnych znaków separujących. Symulator zawiera panel dekodera — wpisz ciąg bitów, a on parsuje go symbol po symbolu, podświetlając ścieżkę drzewa dla każdego zdekodowanego symbolu.

Co dalej

Fala 67 jest już zbudowana i przynosi bibliotece przełomowy kamień milowy: granicę 600 symulacji. Sześć dodatków obejmuje teorię chaosu (diagram bifurkacji odwzorowania logistycznego ze stałą Feigenbauma), optykę (wzory mory i aliasing Nyquista), geometrię obliczeniową (triangulację Delaunaya), modelowanie transportu (automat komórkowy ruchu Nagela-Schreckenberga), algorytmy (podział przestrzenny quadtree) oraz mechanikę orbitalną (siły pływowe z cyklami przypływów wiosennych i martwych). Przeczytaj o tym wszystkim w Devlogu #88.