Programiranje u fizici
- 5. DEO -
U ovom delu će biti upotrebljen Ojlerov metod za rešavanje nekoliko diferencijalnih jednačina iz mehanike.
1. Slobodan pad u vakuumu
Jedan od tipičnih problema kojim se demonstrira upotreba Ojlerovog metoda na realnim sistemima jeste slobdan pad. Prvo će se obraditi slobodan pad u vakuumu, a zatim i u vazduhu.
Slobodan pad u vakuumu je trivijalna problem opisan sledećom diferencijalnom jednačinom:
$$\frac{dv}{dt}=g$$
Ova diferencijalna jednačina se lako rešava i daje rešenje:
$$v=v_0+gt$$
Pretpostavimo da diferencijalnu jednačinu želimo da rešimo tako da dobijemo vrednost brzine u 5. sekundi. Primenjujući gradivo iz 3. i 4. dela kursa, program koji upotrebom Ojlerovog metoda daje rešenje je dat u sledećem primeru.
Primer 5.1. – Slobodan pad u vakuumu rešen Ojlerovim metodom
Naravno, poklapanje između aproksimativnog i egzaktnog rešenja će biti potpuno, jer se radi o trivijalnoj diferencijalnoj jednačini u kojoj sa desne strane jednalosti figuriše samo konstanta.
2. Slobodan pad u vazduhu
Situacija se značajno otežava ukoliko se slobodan pad odigrava u vazduhu, tj. kada postoji sila otpora sredine. Neka se posmatra situacija u kojoj padobranac od oko 80 kg skače iz aviona. U tom slučaju važi sledeća jednakost:
$$F=mg+F_o$$
Za male brzine, sila otpora sredine, \(F_o\) je srazmerna prvom stepenu brzine, dok je za veće brzine srazmerna drugom stepenu brzine. U slučaju slobodnog npr. pada padobranca iz aviona, brzina spada u veće vrednosti, pa sila otpora sredine ima sledeći oblik:
$$F_o=k\cdot v^2$$
Jednačina koja opisuje slobodan pad u vazduhu tada postaje:
$$F=mg-k\cdot v^2 $$
dok je njen odgovarajući diferencijalni oblik:
$$m\frac{dv}{t}=mg-k\cdot v^2 $$
odnosno nakon jednostavnog sređivanja izraza:
\begin{equation}
\frac{dv}{dt}=g+\frac{k}{m} \cdot v^2
\label{eq:dj_slpad_vazduh}
\end{equation}
Suštinski, modelovanje slobodnog pada se svodi na rešavanje diferencijalne jednačine \eqref{eq:dj_slpad_vazduh}, međutim, doći do njenog egzaktnog rešenja nije trivijalno. Njeno egzaktno rešenje je:
\begin{equation}
v=\sqrt{\frac{mg}{k}}\cdot \tanh(gt\cdot \sqrt{\frac{k}{mg}})
\label{eq:dj_slpad_vazduh_egz}
\end{equation}
Brzina padobranca nakon skoka iz aviona se brzo povećava, međutim ona se zbog sile otpora vazduha povećava do izvesne vrednosti nakon čega se više ne menja. Drugim rečima, vrednost brzine u saturira i iz izraza \eqref{eq:dj_slpad_vazduh} lako se može doći do te vrednosti. Kada se sila otpora sredine izjednači sa gravitacionom silom, ubrzanje je nula. Brzina se od tog momenta više ne povećava, a dostignuta vrednost se naziva terminalna brzina, \(v_t\):
$$mg=kv_t^2$$
odakle sledi izraz za \(v_t\):
\begin{equation}
v_t=\sqrt{\frac{mg}{k}}
\label{eq:slpad_vazduh_term}
\end{equation}
Sada se egzaktno rešenje može izraziti preko terminalne brzine na sledeći način:
\begin{equation}
v=v_t \cdot \tanh(\frac{gt}{v_t})
\label{eq:dj_slpad_vazduh_egz_term}
\end{equation}
Dalje će biti dat primer programa gde je iskorišćen Ojlerov metod kako bi se grafički prikazala zavisnost brzine od vremena, tokom prvih 25 sekundi slobodnog pada padobranca.
Primer 5.2. – Slobodan pad u vazduhu
Forma programa je identična kao što je to predstavljeno u delu 2 ovog kursa. Dakle, prvo se definiše nova funkcija (linije 4 i 5) koja je zapravo desna strana jednakosti diferencijalne jednačine (u ovom slučaju jednačine \eqref{eq:dj_slpad_vazduh}). Ova funkcija se posle poziva u petlji kojom se računaju nove vrednosti brzine u skladu sa Ojlerovom formulom. Pre petlje su definisane vrednosti konstanti, uz napomenu da koeficijent otpora za oblik čoveka od oko 80 kg iznosi 0,24.
Za početak, vremenski interval će biti podeljen na 21 tačku, odnosno 20 podintervala (\(n=21\), linija 14). Funkcija np.linspace() je iskorišćena u liniji 16, kako bi se definisale vrednosti vremena za koje će se računati vrednosti brzine, a u sledećoj liniji je definisano kako se računa korak. Linija 24 u for petlji je praktično glavna, jer određuje kako se računaju nove vrednosti brzine prema Ojlerovom metodu. U linijama 26-31 su instrukcije za crtanje grafika.
Malo bliži pogled na dobijenu krivu otkriva da je malo izlomljena, a glađa kriva se može dobiti prostim povećanjem vrednosti \(n\). Dalja brza analiza dobijenog grafika ukazuje na to da posle 10. sekunde vrednost brzine saturira na vrednosti od oko 55 km/s. I zaista, kada se iskoristi izraz \eqref{eq:slpad_vazduh_term}, dobija se vrednost \(v_t=57,2 \: \mathrm{km/s}\).
Iako grafički prikaz jeste zgodan za inicijalnu analizu, sa njega je malo problematično očitati tačne vrednosti, što je nekad od velikog interesa. Međutim, dodavanjem dve linije koda, komentarisane linije 33 i 34 u poslednjem primeru, lako se dolazi do vrednosti \(v_t\). Naime, ugrađena funkcija max() daje maksimalnu vrednost iz tražene liste. Vrednost terminalne brzine je definisana u liniji 33, dok je štampanje ove linije instruisano linijom 34. Čitalac se upućuje da “odkomentariše” ove linije kako bi se dobila vrednost terminalne brzine. Čak i veoma niska vrednost za parametar \(n=21\) će dati praktično identičnu vrednost terminalne brzine kao izraz \eqref{eq:slpad_vazduh_term}.
Grafički prikaz je od velikog značaja, ali je prethodno neophodno utvrditi da li se i koliko se slažu aproksimativno i egzaktno rešenje. U te svrhe korisno je nacrtati oba rešenja i vizuelno ih uporediti. Program kojim se mogu uporediti ta dva rešenja je dat u sledećem primeru.
Primer 5.3. – Upoređivanje aproksimativnog i egzaktnog rezultata
U odnosu na prethodni primer, poslednji primer sadrži novu definisanu funkciju (linije 7 i 8). To je funkcija kojom će se računati vrednosti brzine prema egzaktnom rešenju \eqref{eq:dj_slpad_vazduh_egz}. Skoro sve ostalo je isto, uz dodatak linije 30, gde se instruiše da se crta i kriva koja predstavlja egzaktno rešenje. Naravno, dodati su i parametri za crtanje legende.