Oscylator w komputerze: Verlet, energia i uciekająca faza

Potrzebne wiadomości: Ruch harmoniczny, energia sprężyny i pochodna. Pierwsza lekcja pakietu wyjaśnia krok czasu oraz błąd globalny.

Układ zachowawczy ma własny test poprawności

Masa m porusza się poziomo na idealnej sprężynie o stałej k. x mierzymy od równowagi, pomijamy opór i wymuszenie. Z ma = −kx wynika a(x) = −ω²x, gdzie ω = √(k/m). Dla x(0) = A, v(0) = 0 dokładne rozwiązanie to x = A cos(ωt), v = −Aω sin(ωt), a energia E = mv²/2 + kx²/2 pozostaje równa kA²/2.

Jawny Euler używa jednocześnie x_n+1 = x_n + hv_n i v_n+1 = v_n − hω²x_n. Po podniesieniu obu do kwadratu wyrazy mieszane się znoszą:

E_{n+1}=(1+h^2\omega^2)E_n.

Dla każdego h > 0 energia takiego dyskretnego oscylatora rośnie. Mniejszy krok spowalnia wzrost w ustalonym przedziale czasu, ale nie zapewnia stabilnych drgań przez dowolnie wiele okresów. Otrzymana energia nie pochodzi z fizycznej siły; wytworzyło ją przybliżenie.

Dwie połowy zmiany prędkości

Dla przyspieszenia zależnego od położenia możemy w każdym kroku najpierw zmienić prędkość o połowę impulsu, potem przesunąć ciało i zakończyć drugą połową impulsu:

v_{n+1/2}=v_n+\frac{h}{2}a(x_n),

x_{n+1}=x_n+h v_{n+1/2},

v_{n+1}=v_{n+1/2}+\frac{h}{2}a(x_{n+1}).

To metoda Verleta z prędkością, zwana również velocity Verlet. Po wyeliminowaniu prędkości pośredniej dostajemy x_n+1 = x_n + hv_n + h²a(x_n)/2 oraz v_n+1 = v_n + h[a(x_n) + a(x_n+1)]/2. Nowe przyspieszenie obliczamy w nowym położeniu; zastąpienie go starym niszczy ten schemat.

Pośrednie v_n+1/2 jest etapem algorytmu, nie dokładnie znaną prędkością ciągłą w połowie czasu. Dla gładkiej siły błąd globalny x i v na ustalonym przedziale jest rzędu O(h²). Kolejność półimpuls–przesunięcie–półimpuls jest symetryczna: wykonanie kroku z −h odwraca te operacje w arytmetyce dokładnej. Opór zależny od prędkości wymaga innego schematu; nie wstawiamy go bezrefleksyjnie do tej jawnej postaci.

Pierwszy krok ze wszystkimi jednostkami

Dla m = 1 kg, k = 4 N/m, A = 0,50 m i h = 0,25 s mamy ω = 2 rad/s, a₀ = −2 m/s². Pierwszy półimpuls daje v_1/2 = −0,25 m/s. Przesunięcie: x₁ = 0,50 + 0,25(−0,25) = 0,4375 m. Nowe przyspieszenie jest −1,75 m/s², a drugi półimpuls daje v₁ = −0,46875 m/s. Energia wynosi 0,492675781 J wobec początkowych 0,50 J.

Wzorzec w t = 0,25 s daje x = 0,50 cos 0,5 = 0,438791 m i v = −sin 0,5 = −0,479426 m/s. Verlet nie jest dokładnym rozwiązaniem nawet wtedy, gdy wykres wygląda dobrze. Euler w tym samym kroku dałby x₁ = 0,50 m, v₁ = −0,50 m/s i E₁ = 0,625 J.

Granica stabilności i to, co rzeczywiście się zachowuje

Wprowadźmy q = hω. Eliminacja prędkości daje x_n+1 − 2x_n + x_n−1 = −q²x_n. Rozwiązanie oscylacyjne można zapisać z kątem przyrostu fazy θ takim, że cos θ = 1 − q²/2. Stabilny, nieosobliwy zakres oscylatora to 0 < q < 2. Sam punkt q = 2 nie jest bezpieczny dla ogólnych danych początkowych: zaburzenie może dać wzrost liniowy. Dla q > 2 pojawia się wzrost wykładniczy.

Dla stałego kroku i tego liniowego oscylatora dokładnie zachowuje się w arytmetyce dokładnej inna forma kwadratowa:

\widetilde E=\frac{m}{2}\left[v_n^2+\omega^2\left(1-\frac{q^2}{4}\right)x_n^2\right].

Można to sprawdzić przez podstawienie trzech aktualizacji. Ta pomocnicza wielkość nie jest fizyczną energią sprężyny. Dla startu z A i v₀ = 0 mamy x_n = A cos(nθ), stąd E_n/E₀ = 1 − (q²/4)sin²(nθ). Energia fizyczna oscyluje w ograniczonym przedziale [1 − q²/4; 1] razy E₀. Przy naszym q = 0,5 jest to od 0,46875 do 0,50 J.

Energia może wyglądać dobrze, kiedy faza jest już zła

Numeryczna częstość kołowa wynosi ω_num = θ/h = 2arcsin(q/2)/h. Jest większa od ω dla 0 < q < 2. Kolejne maksima pojawiają się więc za wcześnie. Ten błąd fazy kumuluje się z czasem, mimo ograniczonego błędu energii.

Przy ω = 2 rad/s i h = 0,25 s otrzymujemy ω_num = 2,021442 rad/s. Po 20 s przesunięcie fazy wynosi około 0,428841 rad. Dla h = 0,125 s spada do 0,104906 rad, blisko czterokrotnie. Porównujmy więc energię, położenie oraz fazę w interesującym nas przedziale; żadna pojedyncza liczba nie zastępuje wszystkich kontroli.

Zadania do tego tematu

• Drugi krok sprężyny: porównaj stan, nie tylko amplitudę
• Czy spadek energii w tabeli dowodzi tarcia?
• Dobra energia, zła godzina maksimum

Dalsza lektura

Materiał referencyjny