2-6-Diferencijalne_jednadzbe-S
system:sage


<h3>Vektorska polja</h3>
<p>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.</p>
<p>Neka se čestica u ravnini nalazi u potencijalu $U(x,y)= -4 x^2 y^3$.&nbsp; Odredimo odgovarajuću silu ${\bf F} = -{\bf \nabla} U$</p>
<p><span style="font-size: medium;"><br /></span></p>

{{{id=178|
var('x y z')
U = -4*x^2*y^3
F = -U.gradient(); F
///
}}}

<p><span style="font-size: medium;">Dakle, polje sile je ${\bf F} = 8 x y^3\, {\bf i} + 12 x^2 y^2\, {\bf j}$&nbsp; i možemo ga skicirati funkcijom plot_vector_field():</span></p>

{{{id=181|
plot_vector_field(F, (x,-2,2), (y,-2,2))
///
}}}

<p>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 &scaron;to ne ovisi o $z$, onda imamo:</p>

{{{id=297|
F3 = -U.gradient([x,y,z]); F3
///
}}}

<p>Sage nema specijalizirane funkcije za divergenciju i rotaciju, ali nije ih te&scaron;ko napisati. Jedan primjer nalazi se u datoteci dostupnoj ovako:</p>

{{{id=201|
load('http://www.phy.hr/~kkumer/sp/divcurl.py')
///
}}}

{{{id=176|
print div(F3, [x,y,z])
curl(F3, [x,y,z])
///
}}}

<p>Provjerimo da je divergencija rotacije proizvoljne funkcije nula:</p>

{{{id=186|
div(curl([x^2 + 2*y, x^3 + sin(z), y*z + 2], [x,y,z]), [x,y,z])
///
}}}

<p>Crtanje 3D vektorskog polja:</p>

{{{id=270|
plot_vector_field3d(F3, (x,-2,2), (y,-2,2), (z,-2,2))
///
}}}

<h2>Diferencijalne jednadžbe</h2>

<h3>Simboličko rje&scaron;avanje</h3>

<p>Osnovna Sage funkcija koja traži analitička rje&scaron;enja običnih diferencijalnih jednadžbi je desolve(). Rije&scaron;imo pomoću nje jednadžbu gu&scaron;enog harmoničkog oscilatora</p>
<p>$$ \frac{d^2 y}{d t^2} + 2 \frac{dy}{dt} + 4y = 0 \,.$$</p>
<p>Prilikom definiranja diferencijalne jednadžbe treba eksplicitno deklarirati funkciju funkcijom function()</p>

{{{id=59|
t = var('t')
y = function('y',t)
desolve(diff(y, t, 2) + 2*diff(y, t) + 4*y == 0, y)
///
}}}

<p>Napomene:</p>
<ol>
<li>Prisjetite se da je znak jednakosti u jednadžbama '=='</li>
<li>Kao &scaron;to je kod običnih jednadžbi trebalo eksplicirati po kojoj varijabli se traži rje&scaron;enje, tako se kod diferencijalnih specificira po kojoj funkciji se traži rje&scaron;enje</li>
<li>Rje&scaron;enja uključuju nepoznate konstante ($k_1$ i $k_2$) koje treba odrediti iz početnih uvjeta</li>
</ol>
<p>Ukoliko želimo partikularno rje&scaron;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)$]:</p>

{{{id=60|
des = desolve(diff(y, t, 2) + 2*diff(y, t) + 4*y == 0, y, ics=[0, 1, 1]); des
///
}}}

{{{id=69|
plot(des,(0,5)) + text('Guseni h.o.', (2,0.5))
///
}}}

<p><span id="cell_outer_18"><span id="cell_outer_10"><span id="cell_outer_5">&diams; </span></span><strong>Zadatak 2-6.1</strong><strong>: </strong></span>Rije&scaron;ite simbolički diferencijalnu jednadžbu</p>
<p>$$ \frac{dy}{dt} - y^2 - 1 = 0 \;, $$</p>
<p>uz početni uvjet $y(0)=1/2$, provjerite rje&scaron;enje uvr&scaron;tavanjem te ga nacrtajte za t u rasponu $0&lt;t&lt;1$.</p>


<h3>Numeričko rje&scaron;avanje</h3>
<p>Neke jednadžbe se ne mogu rije&scaron;iti analitički, pa moramo pribjeći numeričkom integriranju, npr. funkcijom desolve_system_rk4(). (rk4 = Runge-Kutta metoda 4. reda, &scaron;to je vrlo popularan algoritam). Npr. na kolegiju Diferencijalne jednadžbe i dinamički sustavi proučava se tzv. van der Polova jednadžba.</p>
<p>$$ \frac{d^2 x}{d t^2} - (1-x^2) \frac{dx}{dt} + x = 0 \,.$$</p>
<p>Jednadžbe vi&scaron;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</p>
<p>$$ \frac{dx}{dt} = y$$</p>
<p>$$ \frac{dy}{dt} = (1-x^2)y -x $$</p>
<p>Prvi argument funkcije desolve_system_rk4&nbsp; 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:</p>

{{{id=50|
var('y')  # gore je y bio function()
res = desolve_system_rk4([y,  (1 - x^2)*y  - x], [x, y], ics=[0, 1, 1], ivar=t, end_points=15)
///
}}}

<p>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.)</p>
<p>Prvih nekoliko točaka možemo očitati koristeći tzv slice sintaksu (vidi kasnije odjeljak 3.1):</p>

{{{id=214|
print 'Broj točaka = ' + str(len(res))
res[:7]
///
}}}

<p>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:</p>

{{{id=197|
P1 = list_plot([[t,x] for t,x,y in res], plotjoined=True)
P2 = list_plot([[t,y] for t,x,y in res], plotjoined=True, color='red')
P1+P2
///
}}}

<p><span id="cell_outer_426"><span id="cell_outer_18"><span id="cell_outer_10"><span id="cell_outer_5">&diams; </span></span><strong>Zadatak 2-6.2</strong><strong>: </strong></span></span>Nacrtajte graf brzine gornjeg rje&scaron;enja Van der Polove jednadžbe y(t)=dx(t)/dt, ali ne u ovisnosti o vremenu t, već&nbsp; o položaju x(t)&nbsp; (tzv. fazni dijagram).</p>


<p><span id="cell_outer_426"><span id="cell_outer_18"><span id="cell_outer_10"><span id="cell_outer_5">&diams; </span></span><strong>Zadatak 2-6.3</strong><strong>: </strong></span></span>Proučite funkciju spline() koja listu točaka pretvara u funkciju koja interpolira te točke, te nacrtajte gornja dva dijagrama rje&scaron;enja Van der Polove jednadžbe (obični i fazni) pomoću plot() i parametric_plot().</p>

<p><span id="cell_outer_426"><span id="cell_outer_18"><span id="cell_outer_10"><span id="cell_outer_5">&diams; </span></span><strong>Zadatak 2-6.4</strong><strong>: </strong></span></span>Rije&scaron;ite diferencijalnu jednadžbu iz zadatka 2-6.1 numerički, nacrtajte rje&scaron;enje, te usporedite sa simboličkim rje&scaron;enjem. Smanjite korak integracije tako da krivulja bude glatka.</p>


<p><span id="cell_outer_18"><span id="cell_outer_10"><span id="cell_outer_5">&diams; </span></span><strong>Zadatak 2-6.5</strong><strong>: </strong></span>Neka na česticu u ravnini djeluje&nbsp; 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&scaron;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&scaron;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.</p>


<h3>(*) Pobolj&scaron;ano numeričko rje&scaron;avanje</h3>

<p>Ne&scaron;to brži i moćniji algoritmi za numeričko rje&scaron;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</p>
<p>$$ \frac{d y_i}{dt} = f_i(y, t) \quad i=1, 2, \ldots \;, $$</p>
<p>te je potrebno definirati funkciju (vidi slijedeće poglavlje) koja vraća listu $[f_0(t), f_1(t), \ldots]$, gjde je dopu&scaron;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$):</p>

{{{id=240|
def syst(t, y, params):
    return [y[1], params[0]*(1 - y[0]^2)*y[1] - y[0]]
///
}}}

<p>Sučelje za ode_solver() je malo drugačije: ta funkcija kreira&nbsp; 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 ]$).</p>

{{{id=241|
S = ode_solver(syst, y_0=[1,0])
///
}}}

<p>Metoda koja rje&scaron;ava jednadžbu je ode_solve(), čiji važni argumenti&nbsp; 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:</p>

{{{id=242|
S.ode_solve(t_span=[0,80], params=[1.7], num_points=100)
///
}}}

<p>Rje&scaron;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:</p>

{{{id=279|
print 'Broj točaka = ' + str(len(S.solution))
S.solution[:3]
///
}}}

{{{id=235|
plot(S.interpolate_solution(i=0), (0,15)) + plot(S.interpolate_solution(i=1), (0,15), color='red')
///
}}}

{{{id=295|
xy = (S.interpolate_solution(i=0), S.interpolate_solution(i=1))
///
}}}

{{{id=294|
@interact
def phase_diag(tmax=(10..80)):
    parametric_plot(xy, (0, tmax)).show()
///
}}}

<p><span id="cell_outer_426"><span id="cell_outer_18"><span id="cell_outer_10"><span id="cell_outer_5">&diams; </span></span><strong>Zadatak 2-6.6</strong><strong>:</strong></span></span>Ponovite zadatak 2-6.5, ali s ode_solver().</p>

{{{id=301|

///
}}}