6. To chce gule vzorové riešenie

Úloha vyzerá celkom priamočiaro – a aj je, akurát si treba dať pozor a nenechať sa oklamať čiastkovým úspechom. Simulovať môžeme buď v nejakom programovacom jazyku, alebo sa uspokojíme s nejakým tabuľkovým procesorom, napríklad Excelom. Simulácia bude v oboch prípadoch rovnaká, akurát za \(C_d\) prvý raz dosadíme konštantu, druhý raz si doň musíme v každom kroku priradiť výsledok pomerne zložitého výpočtu. Ten však počítač urobí za nás, stačí mu raz povedať, ako sa to robí.

Keďže kroky si môžeme voliť ako chceme¹, postačí nám jednoduchá Eulerova metóda: v každom kroku si analyticky spočítame zrýchlenie \(a\). Tu našťastie závisí iba od rýchlosti ako \[ a(v) = \frac{F(v)}{m} = -\frac{C_d S \rho_\mathrm{a} v^2}{2 m} = -\frac{C_d \pi r^2 p M_{\mathrm{a}} v^2}{2 m R T}, \qquad(1)\] kde sme hustotu vyjadrili zo stavovej rovnice ideálneho plynu, pričom sa však nič hrozné nestane, ak za ňu dosadíme tabuľkovú hodnotu.

Teraz v každom kroku rýchlosť zväčšíme o \(a \Delta t\)², polohu zväčšíme o \(v \Delta t\) a tak až donekonečna – alebo teda kým neuvidíme, že sa delová guľa zastavila, resp. jej rýchlosť klesla pod nejakú kritickú hodnotu, povedzme \(\SI{1}{\micro\metre\per\second}\). Teda máme predpis \[ \begin{aligned} x_0 &= \SI{0}{\metre}, \qquad &v_0 &= \SI{100}{\metre\per\second}, \\ x_{n + 1} &\leftarrow x_n + v_n \cdot \Delta t, \qquad &v_{n + 1} &\leftarrow v_n + a(v_n) \cdot \Delta t \end{aligned} \qquad(2)\] a s tým si už program alebo tabuľkový procesor hravo poradia.

Jednoduché? Vyzerá to tak. Lenže pri numerickom riešení diferenciálnych rovníc často natrafíme na nepríjemný jav. Presnosť totiž závisí od veľkosti kroku.

• Ak zvolíme priveľký krok, pri hodnotách, kde sa zrýchlenie prudko mení (v našom prípade na začiatku pohybu, keď je odporová sila veľká), stratíme presnosť. Ak by sme zvolili veľmi veľký krok, guľa by v našom výpočte zastavila prakticky okamžite, alebo sa dokonca začala pohybovať naspäť.
• Ak zvolíme primalý krok, výpočet síce bude veľmi presný… akurát jeho konca sa nemusíme dožiť.

Na tejto rovnici je trochu nepríjemné práve to, že na dosiahnutie dostatočnej presnosti potrebujeme veľmi malý krok³. Ak sa uspokojíme s prvým tipom, ktorý zvolíme, určite nedostaneme dobrý výsledok. Ak \(\Delta t\) trochu zväčšíme, zrazu nám vyjde niečo úplne iné. Aby nám malá zmena \(\Delta t\) dala iba malú zmenu výsledku, musíme použiť krok natoľko malý, že nás pred ukončením výpočtu zastihne dôchodkový vek. Rozumná presnosť sa dá dosiahnuť pre \(\Delta t = \SI{0.001}{\second}\). Pri sekundovom kroku už tratíme veľa presnosti na začiatku pohybu.

S trochou zamyslenia sa ale vieme vynájsť: čo keby sme použili malé časové kroky vtedy, keď sa zrýchlenie prudko mení (teda zo začiatku), a veľké kroky neskôr, keď budú rýchlosť a odporová sila menšie? Napríklad po každom kroku zväčšíme \(\Delta t\) o , čiže ho vynásobíme \(\num{1.01}\). Takto vieme dosiahnuť obstojnú presnosť a zároveň stihnúť termín odovzdania. Rovnica 2 samozrejme prejde na \[ \begin{aligned} x_0 &=& \SI{0}{\metre}, \qquad v_0 &=& \SI{100}{\metre\per\second}, \qquad \Delta t_0 &=& \SI{0.001}{\second}, \\ x_{n + 1} &\leftarrow& x_n + v_n \cdot \Delta t_n, \qquad v_{n + 1} &\leftarrow& v_n + a(v_n) \cdot \Delta t_n, \quad \Delta t_{n + 1} &\leftarrow& \Delta t_n \cdot \num{1.01}. \end{aligned} \qquad(3)\]

Prvý prípad

V prvom prípade by sme mali zistiť, že nech robíme, čo chceme, vzdialenosť nám vždy závisí od počtu krokov, ktoré necháme počítač odsimulovať. Dokonca pekne exponenciálne: ak celkovú dobu simulácie zväčšíme \(e\)-krát, guľa zaletí o skoro presne \(\SI{2500}{\metre}\) ďalej. Guľa sa teda nezastaví nikdy a preletená vzdialenosť môže byť ľubovoľne veľká, ak jej dáme dosť času.

Od vás sme to samozrejme nechceli, ale táto rovnica sa dá vyriešiť aj analyticky: \[ \frac{\mathrm{d} v}{\mathrm{d} t} = -\alpha v^2 \quad\text{s počiatočnou podmienkou}\quad v(0) = \SI{100}{\metre\per\second} \quad\Rightarrow\quad v(t) = \frac{1}{100 \alpha t + 1} \cdot \SI{100}{\metre\per\second}. \]

Pre \(t \to \infty\) síce rýchlosť klesá k nule, ale nie dosť rýchlo. Prejdená dráha je rovná \[ x(t) = \int\limits_{0}^{t} v(\tau) \mathrm{d} \tau = \frac{\ln\left(100 \alpha t + 1\right)}{\alpha}, \] a teda rastie síce veľmi pomaly, ale nad všetky medze.

Druhý prípad

Popísaný vzťah je nechutný, ale zato veľmi dobre sedí s experimentálnymi dátami. Ani v tomto prípade sa guľa nikdy úplne nezastaví, pokles rýchlosti je ale dostatočne rýchly na to, aby celková vzdialenosť bola konečná. Žiadnym nastavením presnosti alebo dĺžky výpočtu nedosiahneme, aby preletela viac ako \(\SI{32.639}{\kilo\metre}\).

Aj keď kompletné analytické riešenie je v tomto prípade nedosiahnuteľné⁴, pre malé \(\mathrm{Re}\) sa funkcia \(C_d\) rovná zhruba \(\frac{24}{\mathrm{Re}}\) a odporová sila sa približuje k Stokesovmu odporu \[ F = - \underbrace{6 \pi \mu r}_{\beta} v. \]

Ako guľa brzdí, skôr či neskôr musí spomaliť natoľko, aby viskozita vzduchu začala hrať dôležitú rolu. Pohyb je potom popísaný rovnicou \[ \frac{\mathrm{d} v}{\mathrm{d} t} = -\beta v \quad\Rightarrow\quad v(t) = v(t_0) \mathrm{e}^{-\beta t}, \]

Rýchlosť klesá k nule, aj keď ani teraz ju nikdy nedosiahne. Prejdená dráha je však rovná \[ x(\infty) = \int\limits_{t_0}^{\infty} v(\tau) \mathrm{d} \tau = \frac{v(t_0)}{\beta}, \] čo celkom veselo konverguje ku konečnej hodnote.

Komentár k riešeniam

Úloha mala dve podstatné podčasti: správne naprogramovať samotné riešenie, a potom si zvoliť správnu dĺžku kroku. Pravda, nešlo o to, či nájdete konkrétnu hodnotu alebo predpis, ktorý dá ako-tak správny výsledok, ale či si všimnete, čo to robí.