Devlog #81 – Fala 61: Zespół podwójnych wahadeł, płyn lepkosprężysty i wizualizator DFT/STFT

Nowe symulacje

⚖️ Zespół podwójnych wahadeł — ilościowe ujęcie efektu motyla

Po co zespół?

Pojedyncza symulacja podwójnego wahadła pokazuje chaos wizualnie — nieprzewidywalny ruch — ale nie mówi nic o tym, jak szybko rozbiegają się bliskie sobie trajektorie. Zespół N=30 wahadeł, wszystkich startujących z tego samego kąta z wyjątkiem maleńkiego przesunięcia Δθ w θ₁, czyni tę szybkość rozbieżności mierzalną. To numeryczny odpowiednik motyla Lorenza: zmień stan początkowy o 10⁻⁶ rad, a obie trajektorie staną się całkowicie nieskorelowane po ~15 sekundach.

Równania Lagrange'a (m₁=m₂=l₁=l₂=1, g=9,81)

denom = 3 - cos(2θ₁ - 2θ₂)

θ₁'' = [ -3g·sinθ₁ - g·sin(θ₁-2θ₂) - 2sin(θ₁-θ₂)·(θ₂'² + θ₁'²cos(θ₁-θ₂)) ] / denom

θ₂'' = [ 2sin(θ₁-θ₂)·(2θ₁'² + 2g·cosθ₁ + θ₂'²cos(θ₁-θ₂)) ] / denom

Każdy element jest całkowany niezależnie metodą RK4 z krokiem Δt = 0,005 s. Element referencyjny (k=0) używa dokładnie θ₁ = θ₁₀. Element k używa θ₁₀ + k·Δθ. Rozrzut zespołu Δθ₁(t) = max|θ₁(k,t) − θ₁(0,t)| jest śledzony w każdej klatce.

Szacowanie wykładnika Lapunowa

Dla układu chaotycznego średnia szybkość wzrostu infinitezymalnych różnic wynosi δ(t) ~ δ(0) · e^λt, więc log δ rośnie liniowo z nachyleniem λ (bity/s przy log₂, naty/s dla ln). Symulacja gromadzi log₁₀|Δθ|_max w czasie, a następnie dopasowuje regresję liniową przez środkowe 75% historii (odrzucając początkowy stan przejściowy i plateau nasycenia). Nachylenie × log₂(10) daje λ w bitach/s. Dodatnie λ potwierdza chaos; dla regularnych warunków początkowych zbiega ono do zera.

Kodowanie kolorem i statystyki

• Wachlarz zespołu — Każdy element jest rysowany z hsl(270°…60°), stopniowo przechodząc od fioletu (referencja) do żółci (najbardziej przesunięty). Szerokość wachlarza to bezpośredni wizualny wskaźnik separacji.
• Liczba spójnych — Liczone są elementy, których |Δθ₁| jest poniżej 0,1 rad; gdy spadnie to do 1, układ znajduje się w reżimie w pełni chaotycznym.
• Etykieta reżimu — Spójny (wszystkie się pokrywają), Rozbieżny (wachlarz się rozszerza), Chaotyczny (nasycona separacja).
• Ustawienia gotowe — Mały rozrzut (Δθ = 10⁻⁴), duży rozrzut (Δθ = 0,2), kąty krzyżujące się (θ₂ = 2,5), ultra-precyzyjny (Δθ = 10⁻⁷).

💧 Płyn lepkosprężysty — reologia materiałów polimerowych

Lepkosprężystość: stan pośredni

Większość rzeczywistych materiałów nie jest ani czysto sprężysta (sprężyny pamiętające swój kształt), ani czysto lepka (płyny, które go zapominają). Polimery, hydrożele, tkanki biologiczne, Silly Putty i farba wykazują wszystkie zachowanie lepkosprężyste: reagują sprężyście na szybkie odkształcenia (wysoka liczba Deborah) i lepko na wolne (niska liczba Deborah). Liczba Deborah De = τ_R/t_obs to kluczowy stosunek bezwymiarowy, gdzie τ_R = η/G₀ to czas relaksacji Maxwella.

Trzy modele i ich równania

Wszystkie trzy modele opisują element materiału poddany naprężeniu σ i odkształceniu ε:

• Model Maxwella (sprężyna + tłumik szeregowo): naprężenie relaksuje wykładniczo przy stałym odkształceniu. G(t) = G₀ · e^−t/τᴿ. Czysto lepki przy długich czasach: brak granicy pełzania — odkształcenie rośnie nieograniczenie. G′(ω) = G₀(ωτ)²/[1+(ωτ)²], G″(ω) = G₀ωτ/[1+(ωτ)²].
• Model Kelvina-Voigta (sprężyna + tłumik równolegle): brak relaksacji naprężeń, ale pełzanie się nasyca. J(t) = (1 − e^−t/τᴿ) / G₀. Czysto sprężysty przy wysokich częstotliwościach: G′(ω) = G₀, G″(ω) = ηω.
• Standardowe Ciało Liniowe (SLS / model Zenera): dodaje trzeci element (sprężynę równoległą do ramienia Maxwella). Łączy ograniczone pełzanie z prawidłową relaksacją naprężeń. G(t) = G₀[α + (1−α)e^−t/τᴿ], gdzie α to stosunek modułu równowagowego. Najdokładniejszy dla polimerów i tkanek biologicznych.

Symulacja cząstek

Lewe płótno pokazuje siatkę 10×10 punktów masy połączonych sprężynami Maxwella (sprężyna G₀ szeregowo z tłumikiem η). Każde wiązanie śledzi swoje wewnętrzne wydłużenie lepkie L_v:

F = G₀_bond × (L − L₀ − Lv)       // siła sprężyny
dLv/dt = F / η_bond                 // relaksacja tłumika (τR = η/G₀)
overdamped: ζ × v = ΣF_springs      // brak bezwładności; natychmiastowa prędkość

Gdy przeciągasz cząstkę, sąsiednie wiązania rozciągają się (od niebieskiego do czerwonego wg wielkości naprężenia), a następnie powoli relaksują z powrotem po zwolnieniu. Ogólne zachowanie bryły przypomina model Maxwella: im szybciej ją odkształcasz, tym sztywniejsza się wydaje.

🎵 Wizualizator DFT i STFT — od sinusoid do spektrogramów

Dyskretna transformata Fouriera

Mając N próbek w dziedzinie czasu x[n], DFT oblicza N zespolonych współczynników:

X[k] = Σₙ x[n] · w[n] · e^{−j2πkn/N}    k = 0, 1, …, N/2

gdzie w[n] to funkcja okna. |X[k]| to amplituda k-tego przedziału częstotliwości, o częstotliwości f_k = k · f_s/N Hz. Wizualizator oblicza to bezpośrednio (O(N²)) dla N ≤ 1024 przy 60 fps — wystarczająco szybko dla wyświetlania w czasie rzeczywistym.

Funkcje okna i przeciek widmowy

Niecałkowita liczba cykli w oknie DFT powoduje „przeciek" energii z ostrych szczytów widmowych do sąsiednich przedziałów. Wybór okna to kompromis między szerokością listka głównego (rozdzielczość częstotliwościowa) a poziomem listków bocznych (przeciek):

• Prostokątne — Najwyższa rozdzielczość, najgorszy przeciek (−13 dB pierwszy listek boczny). Najlepsze, gdy sygnał idealnie wypełnia okno.
• Hanning — w[n] = 0,5(1 − cos(2πn/N)). Pierwszy listek boczny przy −31 dB. Standardowy wybór dla audio.
• Hamming — w[n] = 0,54 − 0,46cos(2πn/N). Pierwszy listek boczny przy −41 dB, ale niezerowe punkty końcowe → niewielka nieciągłość.
• Blackman — Suma trzech wyrazów cosinusowych, pierwszy listek boczny przy −57 dB. Najszerszy listek główny; najlepszy do wykrywania słabych tonów blisko silnych.

Spektrogram STFT

Krótkoczasowa transformata Fouriera stosuje DFT do kolejnych nakładających się ramek: skok = N − nakładanie próbek. Wynikiem jest macierz 2D: kolumny to ramki czasowe, wiersze to przedziały częstotliwości. Symulacja utrzymuje bufor pierścieniowy 460 kolumn (odpowiadający szerokości płótna) kolorowany paletą podobną do viridis, mapującą wartość dB (−80…0) na ciemny niebieski → cyjan → żółty.

Wypróbuj ustawienie chirp (częstotliwość zamiatająca 100 → 800 Hz w ciągu 1 sekundy) i obserwuj przesuwający się w górę ukośny żółty grzbiet na spektrogramie. Następnie przełącz na ustawienie wieloczęstotliwościowe i zobacz trzy poziome linie odpowiadające podstawowej częstotliwości i harmonicznym. Zmiana typu okna podczas oglądania czystej sinusoidy z oknem prostokątnym dramatycznie pokazuje przeciek widmowy: wstęgi boczne częstotliwości pojawiają się i znikają w miarę przełączania okien.

Typy sygnałów

• Sinus — Pojedyncza częstotliwość; jeden szczyt w DFT. Idealny do demonstrowania przecieku okna.
• Prostokąt — Szereg Fouriera 4/π Σ sin((2k−1)2πft)/(2k−1); nieparzyste harmoniczne widoczne jako malejące szczyty.
• Piła — Wszystkie harmoniczne o amplitudzie 1/k; bogatsza zawartość harmoniczna niż fala prostokątna.
• Chirp — Liniowe zamiatanie częstotliwości od f₀ do 8f₀ na sekundę; ukośny pasek na spektrogramie.
• Sygnał AM — Modulacja amplitudy: (1 + 0,8sin(2πf_mt))cos(2πf_ct); dwie wstęgi boczne widoczne wokół fali nośnej.
• Szum biały — Płaskie widmo; wszystkie przedziały częstotliwości mają równą (losową) amplitudę.
• Wieloczęstotliwościowy — Częstotliwość podstawowa f₀ plus harmoniczne przy 3f₀, 5f₀, 7f₀ o malejącej amplitudzie.

Co dalej

Fala 62 będzie kontynuować rozbudowę biblioteki. Rozważane kandydatury obejmują symulacje dynamiki płynów, wizualizacje mechaniki statystycznej oraz dodatkowe interaktywne demonstracje fizyczne. Bądźcie czujni.