Learning #32 – Równania różniczkowe w fizyce: oscylatory, chaos, PDE i metody numeryczne

Równanie różniczkowe wiąże funkcję z jej pochodnymi. Rozwiązania opisują, jak systemy fizyczne ewoluują: jak rośnie populacja, jak propaguje się sygnał elektryczny, jak orbituje planeta, jak dyfunduje ciepło, jak kołysze się wahadło. Każda gałąź fizyki jest ostatecznie zbiorem równań różniczkowych i systematycznych technik ich rozwiązywania. Ten wpis buduje te techniki od podstaw, kończąc dyskusją o systemach chaotycznych, w których długoterminowe przewidywanie staje się niemożliwe mimo w pełni deterministycznych równań.

1. ODE pierwszego rzędu: zanik, wzrost i linie fazowe

ODE pierwszego rzędu ma postać dy/dt = f(t, y). Podstawowe techniki rozwiązywania zależą od struktury f.

Separacja zmiennych

Jeśli f(t, y) = g(t) h(y), równanie jest rozdzielne i można je scałkować bezpośrednio: ∫dy/h(y) = ∫g(t) dt. Klasyczne przykłady: rozpad promieniotwórczy (h(y)=y, g(t)=−λ), dający y(t)=y₀e^−λt; oraz równanie logistyczne.

dN/dt = r N (1 − N/K)

N = populacja,  r = wewnetrzne tempo wzrostu,  K = pojemnosc srodowiska

Rozwiazanie dokladne:  N(t) = K / (1 + ((K−N₀)/N₀) e^{−rt})

Punkty stale: N* = 0 (niestabilny),  N* = K (stabilny)
  → wyznaczone z linii fazowej dN/dt = f(N)

Liniowe ODE pierwszego rzędu i czynnik całkujący

Dla dy/dt + p(t)y = q(t), pomnóż obie strony przez czynnik całkujący μ(t) = exp(∫p(t) dt), by lewa strona stała się pochodną dokładną: d(μy)/dt = μ q(t). Ta technika pojawia się w rozładowaniu obwodu RC szeregowego, prawie stygnięcia Newtona i równaniu Langevina dla ruchu Browna.

2. Liniowe ODE drugiego rzędu: oscylatory i rezonans

Prototypem jest tłumiony oscylator harmoniczny: m&ddot;x + b˙x + kx = F(t). Rozwiązanie jednorodne (F=0) zależy od dyskryminanty wielomianu charakterystycznego mλ² + bλ + k = 0.

m x'' + b x' + k x = F cos(ωt)

ω₀ = √(k/m)     (czestosc wlasna)
ζ  = b/(2mω₀)  (wspolczynnik tlumienia)

Pierwiastki charakterystyczne: λ = −ζω₀ ± ω₀√(ζ²−1)

Niedotlumiony (ζ < 1): x(t) = A e^{−ζω₀t} cos(ω_d t + ϕ)
  ω_d = ω₀ √(1−ζ²)   (tlumiona czestosc wlasna)

Krytycznie tlumiony (ζ = 1): x(t) = (A + Bt) e^{−ω₀t}

Nadtlumiony (ζ > 1): x(t) = A e^{λ₁t} + B e^{λ₂t},  λ₁,λ₂ < 0

Dobroc:  Q = ω₀/2ζω₀ = 1/(2ζ)
  Q ≫ 1  → system rezonansowy (dzwon, obwod LC, wneka optyczna)
  Q ≈ 0,5  → krytycznie tlumiony zamykacz drzwi

Oscylator wymuszony i rezonans

Gdy F(t) = F₀cos(ωt), stacjonarne rozwiązanie szczególne ma amplitudę:

A(ω) = F₀/m / √((ω₀² − ω²)² + (2ζω₀ω)²)

Maksymalna amplituda przy:  ω_rez = ω₀ √(1 − 2ζ²)  (dla ζ < 1/√2)

Amplituda szczytowa:  A_max = F₀/(2mζω₀²√(1−ζ²)) ≈ F₀Q/(mω₀²)  (male ζ)

Opoznienie fazowe:  φ(ω) = arctan(2ζωω₀ / (ω₀² − ω²))
  φ → 0    gdy ω → 0  (w fazie z sila)
  φ = π/2 przy ω = ω₀  (kwadratura przy rezonansie)
  φ → π    gdy ω → ∞  (przeciwfaza 180°)

3. Systemy ODE i analiza płaszczyzny fazowej

Wiele systemów fizycznych opisują autonomiczne systemy ODE: dx/dt = f(x,y), dy/dt = g(x,y). Płaszczyzna fazowa (x,y) wizualizuje trajektorie bez jawnego rozwiązywania. Kluczowe narzędzie to linearyzacja wokół punktów stałych.

Punkty stale (x*, y*): f(x*,y*) = g(x*,y*) = 0

Jacobian w (x*,y*):
  J = | ∂f/∂x   ∂f/∂y |
      | ∂g/∂x   ∂g/∂y |

Wartosci wlasne λ₁, λ₂ macierzy J okreslaja lokalna stabilnosc:
  Re(λ) < 0 obie: stabilny wezel/spirala (atraktor)
  Re(λ) > 0 obie: niestabilny wezel/spirala
  Re(λ) przeciwne znaki: punkt siodlowy (niestabilny)
  Re(λ) = 0: centrum (neutralny, potrzebna prawdziwa analiza nieliniowa)

Slad τ = λ₁+λ₂ = tr(J),   Wyznacznik Δ = λ₁λ₂ = det(J)
Wezel stabilny:   τ < 0, Δ > 0, τ² > 4Δ
Spirala stabilna: τ < 0, Δ > 0, τ² < 4Δ

System drapieżnik-ofiara Lotki–Volterry

dx/dt = αx − βxy, dy/dt = δxy − γy (x = ofiara, y = drapieżnik). Dwa punkty stałe: (0,0) (siodło) i (γ/δ, α/β) (centrum w modelu nieliniowym — orbity zamknięte). Wielkość zachowana to V(x,y) = δx − γln x + βy − αln y.

4. Równania różniczkowe cząstkowe: fala, ciepło i Laplace

Gdy funkcja zależy zarówno od przestrzeni, jak i czasu (lub wielu współrzędnych przestrzennych), jej ewolucja jest rządzona przez PDE. Trzy kanoniczne liniowe PDE drugiego rzędu reprezentują różne archetypy fizyczne.

Równanie fali

∂²u/∂t² = c² ∂²u/∂x²

Rozwiazanie d'Alemberta:  u(x,t) = f(x − ct) + g(x + ct)
  f: fala biegnaca w prawo, g: fala biegnaca w lewo
  majac u(x,0) = φ(x), u_t(x,0) = ψ(x):
  f(ξ) = [φ(ξ) − Κ(ξ)]/2,  g(ξ) = [φ(ξ) + Κ(ξ)]/2
  gdzie Κ(x) = (1/c) ∫ ψ dx

Fala stojaca na strunie [0,L] z u(0,t)=u(L,t)=0:
  u_n(x,t) = sin(nπx/L) [A_n cos(ω_n t) + B_n sin(ω_n t)]
  ω_n = nπc/L  (czestosci modow normalnych)

Równanie ciepła

∂u/∂t = α ∂²u/∂x²       α = dyfuzyjnosc cieplna [m²/s]

Rozwiazanie szeregiem Fouriera na [0,L] z u(0,t)=u(L,t)=0:
  u(x,t) = ∑_n B_n sin(nπx/L) exp(−α(nπ/L)² t)

B_n = (2/L) ∫₀ᴸ u(x,0) sin(nπx/L) dx

Charakterystyczna dlugosc dyfuzji:  ℓ ≈ √(αt)
  Dyfuzyjnosc cieplna: Cu ≈ 1,17×10⁻⁴ m²/s, Si ≈ 8×10⁻⁵ m²/s

Równanie Laplace'a

∇²u = 0 rządzi potencjałem elektrostatycznym w wolnej przestrzeni, ustalonym rozkładem temperatury, nieściśliwym bezwirowym przepływem płynu i funkcjami analitycznymi zespolonymi. Rozwiązania (funkcje harmoniczne) spełniają własność wartości średniej i nie mają lokalnych ekstremów wewnątrz dziedziny. Separacja zmiennych we współrzędnych kartezjańskich, biegunowych i sferycznych daje rodziny rozwiązań (wielomiany Legendre'a, harmoniki sferyczne), które są podstawą rozwinięć multipolowych w elektrostatyce i grawitacji.

5. Metody numeryczne: Euler i Runge–Kutta

Większości ODE spotykanych w praktyce nie da się rozwiązać analitycznie. Numeryczne metody całkowania dyskretyzują czas i krok po kroku posuwają rozwiązanie naprzód. Kluczowy kompromis to dokładność (lokalny błąd obcięcia) kontra koszt obliczeniowy.

Majac dy/dt = f(t,y),  krok h,  biezace (t_n, y_n):

k₁ = f(t_n,          y_n)
k₂ = f(t_n + h/2,   y_n + h k₁/2)
k₃ = f(t_n + h/2,   y_n + h k₂/2)
k₄ = f(t_n + h,     y_n + h k₃)

y_{n+1} = y_n + (h/6)(k₁ + 2k₂ + 2k₃ + k₄)

Lokalny blad obciecia: O(h⁵)   Blad globalny: O(h⁴)

Porownanie:
  Euler:        LBO O(h²), BG O(h)   [1 wywolanie f/krok]
  Trapezow:     LBO O(h³), BG O(h²)  [2 wywolania f/krok]
  RK4:          LBO O(h⁵), BG O(h⁴)  [4 wywolania f/krok]
  Dormand-Prince RK45: adaptacyjny, uzywany przez scipy.integrate.solve_ivp

Sztywne ODE (gdzie niektóre skale czasowe są dużo szybsze niż inne) wymagają metod niejawnych (np. wsteczny Euler, Crank–Nicolson) lub niejawnych metod Runge–Kutty, by uniknąć konieczności stosowania bardzo małych kroków czasowych. Klasycznym przykładem sztywnym jest oscylator Van der Pola z dużym parametrem tłumienia μ.

6. Chaos i dynamika nieliniowa

Chaos deterministyczny nie jest losowością — równania są dokładne i całkowicie deterministyczne. To, co czyni system chaotycznym, to czuła zależność od warunków początkowych: bliskie trajektorie rozchodzą się wykładniczo szybko, czyniąc długoterminowe przewidywanie niemożliwym w praktyce, mimo że jest ono możliwe w zasadzie.

dx/dt = σ(y − x)         σ = 10
dy/dt = x(ρ − z) − y   ρ = 28
dz/dt = xy − βz          β = 8/3

Maks. wykladnik Lapunowa λ₁ ≈ +0,906 (dodatni ⇒ chaotyczny)
  ⇒ rozbieznosc: |Δ(t)| ~ |Δ(0)| e^{λ₁ t}
  Horyzont przewidywalnosci: T ≈ ln(Δ_max/|Δ(0)|) / λ₁
  Przy |Δ(0)|=10⁻10, Δ_max=1: T ≈ 23/0,906 ≈ 25 jednostek czasu Lorenza

Wymiar atraktora (Kaplana-Yorke'a): D_KY ≈ 2,06

Diagramy bifurkacyjne i stała Feigenbauma

Gdy parametr kontrolny się zmienia, system nieliniowy może doświadczyć nagłych jakościowych zmian zachowania (bifurkacji). Odwzorowanie logistyczne x_n+1 = r x_n(1−x_n) pokazuje kaskadę podwajania okresu, gdy r rośnie od 3 do 4: stabilny punkt stały → cykl okresu-2 → okresu-4 → … → chaos przy r≈3,57. Stosunek kolejnych przedziałów bifurkacyjnych zbiega do stałej Feigenbauma δ ≈ 4,6692… — uniwersalnej dla wszystkich jednowymiarowych odwzorowań z kwadratowym maksimum (uniwersalność podwajania okresu).

Równania różniczkowe jako intuicja fizyczna

Progresja w tym wpisie — od separacji zmiennych przez całkowanie RK4 po atraktor Lorenza — odzwierciedla historyczny rozwój dynamiki na przestrzeni trzech stuleci. Druga zasada dynamiki Newtona zapoczątkowała program zapisywania fizyki jako równań różniczkowych. Równanie ciepła Fouriera i rozwój metod spektralnych rozszerzyły go na systemy rozłożone. Geometryczne podejście Poincarégo do równań różniczkowych (lata 80. XIX wieku) wprowadziło płaszczyzny fazowe, punkty stałe i pierwsze przebłyski zachowania chaotycznego.

Interaktywne symulacje pozwalają regulować tłumienie i częstotliwość wymuszenia, by śledzić krzywą rezonansową oscylatora harmonicznego w czasie rzeczywistym, obserwować cykl graniczny drapieżnik–ofiara w przestrzeni fazowej i obserwować czułą zależność atraktora Lorenza od warunków początkowych, uruchamiając dwie niemal identyczne trajektorie i obserwując, jak się rozchodzą. To nie są tylko ćwiczenia obliczeniowe; to te same zjawiska fizyczne, z którymi na co dzień mierzą się inżynierowie, biolodzy, klimatolodzy i fizycy.