Vektorska polja

Opis dinamičkih sustava se obično izvodi diferencijalnim jednadžbama. Obično se takav sustav nalazi u nekom potencijalu opisanom skalarnom funkcijom, čiji negativni gradijent daje polje sile. Promotrimo prvo reprezentaciju takvih polja u Sageu.

Neka se čestica u ravnini nalazi u potencijalu $U(x,y)= -4 x^2 y^3$.  Odredimo odgovarajuću silu ${\bf F} = -{\bf \nabla} U$

Dakle, polje sile je ${\bf F} = 8 x y^3\, {\bf i} + 12 x^2 y^2\, {\bf j}$  i možemo ga skicirati funkcijom plot_vector_field():

Ukoliko potencijal ovisi i o nekim parametrima moramo eksplicitno deklarirati varijable koordinatnog sustava. Isto trebamo i kad potencijal ne ovisi o nekoj koordinati. Npr. ukoliko je gornji potencijal zapravo $U(x, y, z)$, dakle u 3D prostoru samo što ne ovisi o $z$, onda imamo:

Sage nema specijalizirane funkcije za divergenciju i rotaciju, ali nije ih teško napisati. Jedan primjer nalazi se u datoteci dostupnoj ovako:

Provjerimo da je divergencija rotacije proizvoljne funkcije nula:

Crtanje 3D vektorskog polja:

Diferencijalne jednadžbe

Simboličko rješavanje

Osnovna Sage funkcija koja traži analitička rješenja običnih diferencijalnih jednadžbi je desolve(). Riješimo pomoću nje jednadžbu gušenog harmoničkog oscilatora

$$ \frac{d^2 y}{d t^2} + 2 \frac{dy}{dt} + 4y = 0 \,.$$

Prilikom definiranja diferencijalne jednadžbe treba eksplicitno deklarirati funkciju funkcijom function()

Napomene:

1. Prisjetite se da je znak jednakosti u jednadžbama '=='
2. Kao što je kod običnih jednadžbi trebalo eksplicirati po kojoj varijabli se traži rješenje, tako se kod diferencijalnih specificira po kojoj funkciji se traži rješenje
3. Rješenja uključuju nepoznate konstante ($k_1$ i $k_2$) koje treba odrediti iz početnih uvjeta

Ukoliko želimo partikularno rješenje za neke konkretne početne uvjete, zadamo te uvjete u dodatnom argumentu ics (=initial conditions), a koji je lista oblika ics = [$t_0,\, y(t_0),\, y'(t_0)$]:

♦ Zadatak 2-6.1: Riješite simbolički diferencijalnu jednadžbu

$$ \frac{dy}{dt} - y^2 - 1 = 0 \;, $$

uz početni uvjet $y(0)=1/2$, provjerite rješenje uvrštavanjem te ga nacrtajte za t u rasponu $0<t<1$.

Numeričko rješavanje

Neke jednadžbe se ne mogu riješiti analitički, pa moramo pribjeći numeričkom integriranju, npr. funkcijom desolve_system_rk4(). (rk4 = Runge-Kutta metoda 4. reda, što je vrlo popularan algoritam). Npr. na kolegiju Diferencijalne jednadžbe i dinamički sustavi proučava se tzv. van der Polova jednadžba.

$$ \frac{d^2 x}{d t^2} - (1-x^2) \frac{dx}{dt} + x = 0 \,.$$

Jednadžbe višeg reda se za numeričku analizu treba prirediti tako da se pretvore u sustav jednadžbi prvog reda. To se u ovom slučaju radi uvođenjem nove zavisne varijable $y(t)=dx/dt$ pa je gornja jednadžba ekvivalentna sustavu

$$ \frac{dx}{dt} = y$$

$$ \frac{dy}{dt} = (1-x^2)y -x $$

Prvi argument funkcije desolve_system_rk4  je lista desnih strana ovih jednadžbi, drugi argument je odgovarajuća lista funkcija čije prve derivacije su na lijevoj strani, treći argument su početni uvjeti kao gore, zatim se eksplicitno navodi nezavisna varijabla (ivar=independent variable) i konačna točka (ili interval) end_points do koje provodimo numeričku integraciju:

Rezultat je lista točaka istog formata kao početni uvjet ics, koji je sam prva točka te liste, a slijedeće točke su evolucija sustava: res = [ics=[$t_0,\, x(t_0),\, y(t_0)$], [$t_1,\, x(t_1),\, y(t_1)$], ... ]. (Broj točaka, rezolucija itd. se može kontrolirati dodatnim opcijama.)

Prvih nekoliko točaka možemo očitati koristeći tzv slice sintaksu (vidi kasnije odjeljak 3.1):

Grafički prikaz ovog sustava možemo izvesti funkcijom list_plot(). Da bismo pretvorili gornju listu res u listu s elementima $[t, x(t)]$ (položaj sustava u vremenu) odnosno u listu s elementima $[t, y(t)]$ (brzina sustava u vremenu), koristimo tzv. obuhvaćanje liste (list comprehension), vrlo elegantan python idiom koji ćemo često sretati:

♦ Zadatak 2-6.2: Nacrtajte graf brzine gornjeg rješenja Van der Polove jednadžbe y(t)=dx(t)/dt, ali ne u ovisnosti o vremenu t, već  o položaju x(t)  (tzv. fazni dijagram).

♦ Zadatak 2-6.3: Proučite funkciju spline() koja listu točaka pretvara u funkciju koja interpolira te točke, te nacrtajte gornja dva dijagrama rješenja Van der Polove jednadžbe (obični i fazni) pomoću plot() i parametric_plot().

♦ Zadatak 2-6.4: Riješite diferencijalnu jednadžbu iz zadatka 2-6.1 numerički, nacrtajte rješenje, te usporedite sa simboličkim rješenjem. Smanjite korak integracije tako da krivulja bude glatka.

♦ Zadatak 2-6.5: Neka na česticu u ravnini djeluje  potencijal $U(x,y)= -4 x^2 y^3$ (potencijal s početka ovog odjeljka). Nacrtajte potencijal U(x,y) i pomoću contour_plot() i pomoću plot3d(). Nadalje, rješavajući numerički Newtonovu diferencijalnu jednadžbu gibanja, pronađite putanju čestice, uz početne uvjete $x(0)=y(0)=-1, v_x(0)=2, v_y(0)=-3$, od trenutka $t=0$ do trenutka $t=1.43$. Nacrtajte tu putanju, a onda je pokušajte superponirati na oba gornja dijagrama potencijala. Igrajte se malo s početnim uvjetima i uvjerite se da je sve OK. Izvrijednite ukupnu energiju za razna vremena i uvjerite se da je sustav zaista konzervativan tj. da je energija konstantna.

(*) Poboljšano numeričko rješavanje

Nešto brži i moćniji algoritmi za numeričko rješavanje diferencijalnih jednadžbi dostupni su putem funkcije ode_solver(). Slično kao gore problem treba organizirati kao sustav diferencijalnih jednadžbi prvog reda oblika

$$ \frac{d y_i}{dt} = f_i(y, t) \quad i=1, 2, \ldots \;, $$

te je potrebno definirati funkciju (vidi slijedeće poglavlje) koja vraća listu $[f_0(t), f_1(t), \ldots]$, gjde je dopuštena i dodatna ovisnost o nekim konstantnim parametrima sustava. Npr. za van der Polovu jednadžbu takva funkcija se može definirati ovako (uspoređujući s ranijim van der Pol sustavom imamo korespondenciju $x\to y_0$ i $y\to y_1$):

Sučelje za ode_solver() je malo drugačije: ta funkcija kreira  ode_solver objekt i sve ostalo se izvodi putem metoda tog objekta. Prilikom kreiranja možemo zadati sustav jednadžbi (putem gornje funkcije) i početne uvjete (putem argumenta y_0=$[y_{0}(0), y_{1}(0), \ldots ]$).

Metoda koja rješava jednadžbu je ode_solve(), čiji važni argumenti  su t_span=$[t_0, t_{\rm max}]$ i broj točaka integracije num_points, a ukoliko sustav ovisi o nekim konstantnim parametrima specificiramo ih opcijom params:

Rješenje je kao i gore lista točaka (u malo drugačijem formatu), a postoji i funkcija interpolate_solution() koja interpolira te točke i pretvara ih u funkciju, tako da ne moramo koristiti spline. Opcionalni argument te funkcije i=k, određuje koju od funkcija $y_k(t)$ želimo:

♦ Zadatak 2-6.6:Ponovite zadatak 2-6.5, ali s ode_solver().