Spadanie liczone krok po kroku: Euler, błąd i stabilność

Potrzebne wiadomości: Druga zasada Newtona, pochodna i funkcja wykładnicza. Pomaga wcześniejsza lekcja „Opór ośrodka i prędkość graniczna”.

Najpierw siły i skala czasu

Oś z kierujemy pionowo w dół. Ciało o stałej masie m spada w nieruchomym ośrodku; przyjmujemy opór −bv, stałe g oraz b > 0. Pomijamy wypór. Model liniowy musi być uzasadniony dla danego zakresu prędkości; dokładniejszy rachunek nie naprawi niewłaściwego prawa oporu. Jednostki b to kg/s, ponieważ bv ma być siłą w niutonach.

\frac{dz}{dt}=v,\qquad \frac{dv}{dt}=g-\frac{v}{\tau},\qquad \tau=\frac{m}{b},\qquad v_g=g\tau.

τ jest czasem charakterystycznym, a v_g prędkością graniczną ze znakiem zgodnym z osią. Dla z(0) = 0 i v(0) = 0 rozwiązanie ma postać:

v(t)=g\tau(1-e^{-t/\tau}),\qquad z(t)=g\tau\left[t-\tau(1-e^{-t/\tau})\right].

Pochodna z jest równa v, a pochodna v jest równa g e^−t/τ, czyli g − v/τ. Oba warunki początkowe też są spełnione. To niezależny wzorzec do kontroli obliczeń, a nie wynik dopasowany do naszej symulacji.

Co znaczy jeden krok jawnego Eulera?

Niech h > 0 oznacza krok w sekundach, t_n = nh, a z_n i v_n przybliżenia w tej samej chwili. Z rozwinięcia v(t+h) = v(t) + h v′(t) + O(h²) zatrzymujemy dwa pierwsze składniki. Tak samo postępujemy z położeniem:

z_{n+1}=z_n+h v_n,\qquad v_{n+1}=v_n+h\left(g-\frac{v_n}{\tau}\right).

Obie prawe strony korzystają ze starej prędkości. Najpierw obliczamy nowe wartości w zmiennych pomocniczych, potem zastępujemy nimi stary stan. Wstawienie już zaktualizowanego v_n+1 do wzoru na z daje inny algorytm. Nie jest równoważnym przestawieniem dwóch poleceń.

Przyspieszenie ma jednostkę m/s², więc ha jest zmianą prędkości, a hv zmianą położenia. Pierwszy krok od spoczynku daje z₁ = 0. Nie znaczy to, że ciało naprawdę nie poruszało się przez h sekund: prostokątna aproksymacja pola pod rosnącym wykresem v(t) pominęła przemieszczenie wewnątrz kroku.

Policzmy i porównajmy ten sam czas

Przyjmijmy m = 0,20 kg, b = 0,10 kg/s, g = 10 m/s². Wtedy τ = 2 s i v_g = 20 m/s. Dla h = 0,50 s prędkości po kolejnych krokach wynoszą 5; 8,75; 11,5625; 13,671875 m/s. Na przykład drugi krok daje v₂ = 5 + 0,50(10 − 5/2) = 8,75 m/s. Położenia wynoszą 0; 2,50; 6,875; 12,65625 m.

Po 2 s rozwiązanie dokładne daje v = 12,642411 m/s i z = 14,715178 m. Euler zawyża tutaj prędkość, a zaniża przebytą drogę. Nie ma sprzeczności: położenie jest sumowane ze starych prędkości, a błędy obu współrzędnych nie muszą mieć tego samego znaku.

Stabilność nie jest dokładnością

Odejmijmy v_g od obu stron rekurencji. Odchylenie d_n = v_n − v_g spełnia:

d_{n+1}=q d_n,\qquad q=1-\frac{h}{\tau},\qquad d_n=q^n d_0.

Aby odchylenie zanikało, potrzebujemy |q| < 1, czyli 0 < h < 2τ. Dla 0 < h ≤ τ i startu poniżej v_g wynik nie przekracza granicy. Dla τ < h < 2τ zbiega do niej naprzemiennie od góry i od dołu: stabilnie, lecz z niefizycznymi oscylacjami. Przy h = 2τ amplituda odchylenia nie maleje; dla większego kroku rośnie. Są to własności algorytmu, bo dokładne d(t) = d(0)e^−t/τ nie zmienia znaku.

Przy h = τ Euler trafia dokładnie w v_g po jednym kroku, choć model ciągły zbliża się do niej asymptotycznie. To dobry przykład stabilnego wyniku, który nadal może być bardzo niedokładny w interesującej nas chwili.

Jak uczciwie zmniejszać krok?

Porównujemy h, h/2 i h/4 przy tym samym czasie końcowym i tych samych danych. W gładkim zadaniu na ustalonym skończonym przedziale globalny błąd Eulera jest rzędu O(h): po dwukrotnym zmniejszeniu dostatecznie małego h powinien w przybliżeniu zmaleć dwukrotnie. Jednokrokowe pominięcie ma rząd O(h²), lecz kroków jest około t/h. Nie nazywamy zatem całej metody metodą drugiego rzędu.

Dochodząc do zadanego T, wybieramy h = T/N albo skracamy wyłącznie ostatni krok do T − t. Nie porównujemy wyników z różnych chwil. Należy przechować więcej cyfr niż w tabeli, a zaokrąglić dopiero wynik. Błąd kroku, zaokrąglenia i błąd fizycznego modelu to trzy różne rzeczy.

Zadania do tego tematu

• Dwie aktualizacje, jedna stara prędkość
• Sztuczne odbijanie się od prędkości granicznej
• Ile daje podzielenie kroku przez dwa?

Dalsza lektura

Materiał referencyjny