Programiranje u fizici
- 3. DEO -
U ovom delu kursa će u programskom jeziku python biti napisani programi za primenu Ojlerovog metoda za rešavanje jednostavne diferencijalne jednačine.
1. Jednostavna diferencijalna jednačina – “ručni metod”
Pre nego što napišemo python kod za rešavanje jednostavne diferencijalne jednačine, podsetićemo se kako se ručno primenjuje ovaj metod. Dakle, primenom Ojlerovog metoda, treba pronaći vrednost diferencijalne jednačine:
\begin{equation}
\frac{dy}{dx}=2 \cdot x
\label{eq:prva_dj}
\end{equation}
kada \(x\) uzima vrednost 1,5. Uzeti da su početne vrednosti \(x_0 = 1\) i \(y_0 = 1\), a da je korak \(h=0,1\). Egzaktno rešenje ove diferencijalne jednačine je \(x^2\).
Naravno, za rešenje ovog zadatka primenjuje se Ojlerova formula, čiji je najopštiji oblik:
\begin{equation}
y_{i+1}=y_i + h \cdot f(x_i,y_i)
\label{eq:ojler}
\end{equation}
Ručnom metodom, dati zadatak je najlakše rešiti pravljenjem tabele, u kojoj se najpre definišu sve vrednosti \(x\) koje će biti iskorišćene u datoj formuli.
| \(n\) | \(x_n\) | \(y_n\) | \(y_{ex}\) |
|---|---|---|---|
| 0 | 1 | 1 | 1 |
| 1 | 1,1 | 1,20 | 1,21 |
| 2 | 1,2 | 1,42 | 1,44 |
| 3 | 1,3 | 1,66 | 1,69 |
| 4 | 1,4 | 1,92 | 1,96 |
| 5 | 1,5 | 2,20 | 2,25 |
Iteracije koje dovode do rešenja – vrednosti diferencijalne jednačine u tački \(x=1,5\)
\begin{equation*}
y_1=y_0+0,1 \cdot (2 \cdot x_0)=1,00+0,1\cdot(2\cdot1)=1+0,2=1,2 \\
y_2=y_1+0,1\cdot(2 \cdot x_1)=1,20+0,1\cdot(2\cdot1,1)=1,2+0,22=1,42 \\
y_3=y_2+0,1 \cdot (2 \cdot x_2 )=1,42+0,1 \cdot (2 \cdot 1,2)=1,42+0,24=1,66 \\
y_4=y_3+0,1 \cdot (2 \cdot x_3 )=1,66+0,1 \cdot (2 \cdot 1,3)=1,66+0,26=1,92 \\
y_5=y_4+0,1 \cdot (2 \cdot x_4 )=1,92+0,1 \cdot (2 \cdot 1,4)=1,92+0,28=2,20
\end{equation*}
U prvoj koloni su numerisane iteracije. Vrednosti koje uzima \(x\) su date u drugoj koloni. Pošto je početna vrednost \(x_0=0\), a korak \(h=0,1\), u drugoj koloni će se naći 6 vrednosti, od 1 do 1,5. U trećoj koloni su predstavljene vrednosti za \(y\) dobijene primenom Ojlerove formule. Sve iteracije primene ove formule su date sa desne strane tabele. Na kraju, u četvrtoj koloni su date egzaktne vrednosti (kvadrat svake od vrednosti \(x\) iz druge kolone). U tabeli su rezultati poslednje iteracije dati masnim fontom.
Vidi se da primena Ojlerovog metoda dovodi do rešenja koje je vrlo blisko egzaktnom rešenju. Naravno, ovo je vrlo jednostavan primer, pa dobro slaganje ne čudi. Međutim, bez obzira na jednostavnost ovog primera, svejedno je neophodno pravljenje tabele i pet iteracija proračuna, što ostavlja dosta prostora za grešku. Da bi se mogućnost greške svela na minimum, pomenutu proceduru je potrebno automatizovati nekim programskim jezikom – u slučaju ovog kursa programskim jezikom python.
Dalje će biti predstavljen program kojim se automatizuje procedura opisana u ovom delu.
2. Jednostavna diferencijalna jednačina – “programerski” pristup
Suštinski gledano, program za rešavanje diferencijalne jednačine \eqref{eq:prva_dj} se oslanja na način razmišljanja kao i ručna metoda. Najpre će biti definisane sve vrednosti koje promenljiva \(x\) može da uzima, a zatim će se kroz jednu for petlju sprovesti iterativno računanje vrednosti promenljive \(x\). Python program je sledeći:
Primer 1.1 – Program za rešavanje diferencijalne jednačine \eqref{eq:prva_dj}
U datom programu su najpre definisane početne vrednosti. U liniji 4 je definisana lista koja sadrži sve vrednostim koje \(x\) može da uzima. U liniji 5 je definisana lista koja sadrži vrednosti promenljive \(y\) i za početak je stavljeno da su sve vrednosti ove liste jednake 0. Ideja je da ove vrednosti kasnije u for petlji zamenjujemo vrednostima koje se dobijaju sukcesivno primenom Ojlerove metode. U liniji 7 je definisana vrednost koraka.
Od linije 10 startuje praktično glavni deo programa – for petlja zahvaljujući kojoj se 5 puta primenjuje Ojlerova formula \eqref{eq:prva_dj}. U liniji 10 se definiše niz vrednosti koje će uzimati brojač i (u ovom konkretnom slučaju to su vrednosti 0,1,2,3 i 4 – range() uvek ide do N-1 broja, pošto je ovde range(0,5) znači da ide do 4). Zatim se u linijama 11 i 12 prvi elementi listi \(x\) i \(y\) zamenjuju početnim vrednostima. U liniji 13 se sprovodi Ojlerova formula.
DOMAĆI 1:
- Dodati liniju koda koja bi štampala sva numerička rešenja.
- Dodati linije koda koje bi izračunale i sva egzaktna rešenja za sve vrednosti iz liste \(x\)
REŠENJE je dato u sledećem prozoru. Klikom na dugme “Code” otvoriće se kod koji rešava domaći. Probajte prvo sami da dođete do rešenja.
Prethodni program demonstrira koliko je jednostavno implementirati Ojlerov metod za numeričko rešavanje diferencijalne jednačine. Naravno, u realnosti se srećemo sa dosta komplikovanijim jednačinama, a jedan od parametara koji se najčešće menja jeste korak \(h\).
Ukoliko bi korak \(h\) uzimao npr. vrednost 0,01, to bi efektivno značilo da broj svih vrednosti u listi \(x\) skače sa 6 na 50
3. Jednostavna diferencijalna jednačina – Opštiji “programerski” pristup
Prethodni program demonstrira kako izgleda implementirati Ojlerov metod za numeričko rešavanje diferencijalne jednačine. U konkretnom primeru, ukoliko bi se želelo dobiti tačnije numeričko rešenje, neophodno je smanjiti vrednost koraka. Međutim, smanjenje koraka bi dovelo do mnogo većeg broja vrednosti koje \(x\) može da uzima. Konkretno, ukoliko bi korak \(h\) uzimao npr. vrednost 0,01, to bi efektivno značilo da broj svih vrednosti u listi \(x\) skače sa 6 na 51 vrednosti. Ručno unošenje svake pojedinačne vrednosti bi bila velika muka, a verovatnoća za pravljenje greške bi se povećavala svakim dodatnim smanjenjem vrednosti koraka.
Iz tog razloga je neophodno uopštiti prethodni kod, korišćenjem mogućnosti programskog jezika python predstavljenih u prethodnim delovima kursa. Najpre će biti naveden opštiji kod, a zatim sledi detaljno objašnjenje svake linije koda.
Primer 1.2 – Uopštenje koda – prvi
Za uopštenje koda predstavljenog u primeru 1.1 na ovoj stranici koristiće se mogućnosti numpy biblioteke. Konkretno, ključna će biti primena funkcije np.linspace(min,max,n), koja deli brojevni prostor između brojeva min i max na n tačaka. Prethodno je definisana konačna vrednost promenljive \(x_f\), tj. njena vrednost za koju se traži vrednost diferencijalne jednačine Primena funkcije linspace je ilustrovana na slici 1.
Sa Slike 1 se uočava sledeće. Ukoliko je interval između dva broja podeljen na \(n\) tačaka, onda postoji \(n-1\) koraka (tj. podintervala) između svake vrednosti. Vrednost koraka se tada može računati po sledećoj formuli:
\begin{equation}
h=\frac{x_f-x_0}{n-1} .
\label{eq:korak}
\end{equation}
Jasno je da što je veća vrednost broja \(n\), biće i niža vrednost koraka \(h\), a time i tačniji numerički rezultat.
U liniji 10 su upotrebom funkcije linspace definisane sve vrednosti koje \(x\) uzima tokom primene Ojlerove metode. Ova jedna linija koda omogućava jednostavno definisanje ogromnog broja tih vrednosti. U liniji 11 koristimo zeros funkciju numpy biblioteke da definišemo niz y koja sadrži \(n\) elemenata od koji svaki ima vrednost 0.
Zatim dolazi na red for petlja, koja ima identičan oblik kao i prethodni program. U liniji 19 definišemo egzaktno rešenje za krajnju tačku \(x_f\), kako bi numeričko i analitičko rešenje mogli da se porede. U te svrhe se definiše promenljiva “razlika”, a u linijama 23-26 se daju instrukcije za štampanje rezultata.
Ukoliko je \(n=6\), dobija se identičan rezultat kao i u prethodnom programu. Svako dalje povećanje ovog broja dovodiće i do boljeg numeričkog rešenja. Čitaocima se ostavlja da probaju različite vrednosti \(n\) i istraže tačnost dobijenih rezultata.
4. Jednostavna diferencijalna jednačina – Dalje uopštenje koda
Još jedan važan korak se može primeniti kako bi se kod uopštio i olakšao za dalju upotrebu. Naime, za aproksimativno rešavanje diferencijalne jednačine, neophodno je diferencijalni deo staviti na jednu strane jednakosti, a ostatak na drugu stranu jednakosti. Tada, u Ojlerovoj formuli za numeričko rešenje, figuriše strana jednakosti koja ne sadrži diferencijalni deo jednačine (npr. desna strana jednačine \eqref{eq:prva_dj}. Kako bi se uopštio program, pogodno je prvo definisati pomenutu stranu jednakosti kao neku novu funkciju, a zatim je tokom for petlje lako pozivati. Na ovaj način, jedan te isti program možemo koristiti iznova za rešavanje drugih diferencijalnih jednačina, naravno uz modifikacju funkcije koju definišemo, početnih uslova i koraka. Takav program je dat u sledećem primeru:
Primer 1.3 – Dalje uopštenje koda za primenu Ojlerovog metoda
U poslednjem primeru u linijama 3 i 4 definisana je funkcija, koja u stvari predstavlja desnu stranu jednačine \eqref{eq:prva_dj}. Pomenuta jednačina se posle poziva u for petlji, tj. računa se za svaku vrednost iz niza \(x\).