Rzut z oporem: metoda Heuna i lądowanie między krokami

Potrzebne wiadomości: Rzut ukośny, składowe wektora i druga zasada Newtona. Metoda Eulera jest wprowadzona w pierwszej lekcji pakietu.

Siła musi być przeciwna do całego wektora prędkości

Oś x jest pozioma, y pionowa w górę. Ciało traktujemy jako punkt materialny o stałej masie m. Przyjmujemy stałe g, nieruchome powietrze i model oporu F_o = −c|v|v, c > 0, bez siły nośnej, obrotu i wyporu. Jednostką c jest kg/m. Dla rzeczywistego przedmiotu stałość tego współczynnika w całym locie wymaga weryfikacji; tutaj definiuje model.

s=\sqrt{v_x^2+v_y^2},\qquad \beta=\frac{c}{m},\qquad a_x=-\beta s v_x,\qquad a_y=-g-\beta s v_y.

s oznacza szybkość, nie drogę. β ma jednostkę m⁻¹, zatem βsv_x ma jednostkę przyspieszenia. Podczas opadania v_y jest ujemne, więc pionowa składowa oporu jest dodatnia. Zapis −cv_y² dawałby błędnie siłę w dół także przy opadaniu. Osobne −cv_x|v_x| oraz −cv_y|v_y| z kolei zwykle nie tworzą wektora dokładnie przeciwnego do v.

Cztery liczby opisują jeden stan

Stan S = (x; y; v_x; v_y) ma pochodną F(S) = (v_x; v_y; a_x; a_y). Składowe wektora stanu mają różne jednostki; dodawanie hF do S jest poprawne, bo w każdej odpowiadającej sobie składowej jednostki się zgadzają.

Metoda Heuna, czyli jawna metoda trapezów, najpierw przewiduje cały następny stan Eulerem, a następnie uśrednia pochodną na początku i w przewidywanym końcu:

K_1=F(S_n),\qquad S^*=S_n+hK_1,\qquad K_2=F(S^*),

S_{n+1}=S_n+\frac{h}{2}(K_1+K_2).

W K₂ trzeba obliczyć nową szybkość z obu przewidzianych składowych prędkości. Nie mieszamy przewidzianego v_x ze starym v_y. Wzór przypomina całkowanie trapezem, lecz używa przewidywanego końca; nie rozwiązuje niejawnego równania trapezów. Dla gładkiego przebiegu, dostatecznie małego kroku i ustalonego przedziału czasu ma globalny błąd O(h²).

Jeden krok z kompletnymi danymi

Niech m = 0,20 kg, c = 0,010 kg/m, g = 10 m/s² oraz S₀ = (0 m; 0 m; 12 m/s; 16 m/s). Początkowa szybkość 20 m/s i β = 0,050 m⁻¹ dają a₀ = (−12; −26) m/s². Dla h = 0,10 s predyktor to S* = (1,20 m; 1,60 m; 10,8 m/s; 13,4 m/s).

W predyktorze s* = √(10,8² + 13,4²) = 17,210462 m/s. Stąd a* = (−9,293649; −21,531009) m/s². Uśrednienie daje S₁ = (1,14 m; 1,47 m; 10,935318 m/s; 13,623450 m/s). Są to przybliżone dane w t = 0,10 s; drugi etap nie oznacza wykonania dwóch fizycznych kroków po 0,10 s.

Lądowanie to zdarzenie, a nie numer wiersza

Start następuje z y = 0 i v_y > 0; nie kończymy symulacji już w tej chwili. Kolejne pełne kroki wykonujemy aż do pierwszego opadającego przejścia z y_n > 0 do y_n+1 ≤ 0. Punkt końcowy może być pod poziomem gruntu. Zwrócenie jego x jako zasięgu wprowadza dodatkowy błąd, zwykle rzędu h, mimo metody drugiego rzędu.

Przy małym kroku i pojedynczym przejściu przez poziom przybliżamy końcowy odcinek linią. Ułamek kroku i zdarzenie obliczamy jako:

\alpha=\frac{y_n}{y_n-y_{n+1}},\qquad t_{hit}=t_n+\alpha h,\qquad R=x_n+\alpha(x_{n+1}-x_n).

Wstawienie α do interpolacji y daje dokładnie zero. Jeśli chcemy również oszacować prędkość w chwili trafienia, interpolujemy jej obie składowe z tym samym α. Jest to przybliżenie, nie dokładne rozwiązanie równania y(t) = 0. Przy gładkim, poprzecznym przejściu z niezerową prędkością pionową interpolacja zachowuje drugi rząd błędu zdarzenia; styczne dotknięcie, odbicie albo zbyt duży krok wymagają osobnej analizy.

Dla danych przykładu otrzymujemy: h = 0,20 s — czas 2,405904 s i zasięg 16,053901 m; h = 0,10 s — 2,414253 s i 16,105305 m; h = 0,05 s — 2,416051 s i 16,116081 m; h = 0,025 s — 2,416512 s i 16,118912 m. Same cyfry po przecinku nie świadczą o dokładności. Ich porównanie pokazuje, które cyfry przestają zależeć od siatki. Iloraz różnic nie musi być idealnie równy 4, bo miejsce przecięcia wewnątrz kroku zmienia się po zagęszczeniu.

Dwa testy fizyczne przed zaufaniem trajektorii

Przy c = 0 metoda daje dokładne położenia i prędkości w węzłach dla stałego g, z dokładnością arytmetyki. Wzorzec bez oporu ma T = 2v_y0/g = 3,2 s i R = v_x0T = 38,4 m. Interpolacja lądowania między węzłami nadal ma własny błąd, chyba że dokładne trafienie wypadnie w węźle.

Druga kontrola dotyczy energii mechanicznej E = m(v_x² + v_y²)/2 + mgy. W tym modelu:

\frac{dE}{dt}=\mathbf F_o\cdot\mathbf v=-c s^3\leq 0.

Praca grawitacji znosi zmianę energii potencjalnej, zostaje ujemna moc oporu. Metoda Heuna nie gwarantuje dyskretnego spadku energii dla dowolnego h. Duże wzrosty wskazują na błąd implementacji albo zbyt gruby krok; sam spadek energii nie dowodzi jednak, że kierunek siły lub zasięg są poprawne.

Zadania do tego tematu

• Dwa etapy Heuna nie są dwoma krokami czasu
• Zasięg między dwoma punktami po przeciwnych stronach gruntu
• Obie składowe oporu muszą mieć wspólną szybkość

Dalsza lektura

Materiał referencyjny