Predavanje 24

FlipView[{Text[Style[" Branko Urošević, imqfbrankourosevic@gmail.com ",20,"Text",Italic]],Text[Style[" Faculty of Economics, University of Belgrade, NBS and CESIfo Munich",20,"Text", Italic]]}]

Sistemi diferencnih i diferencijalnih jednačina. Numeričko rešavanje diferencijalnih jednačina

Sistemi diferencnih jednačina prvog reda

U ekonomiji, evolucija sistema tipično podrazumeva da više varijabli imaju neku međuzavisnost. Pri tome, vrednosti varijabli u trenutku t mogu zavisi od vrednosti varijabli u prethodnom trenutku t-1. U tom slučaju kažemo da je dinamika zadata pomoću sistema diferencnih jednačina prvog reda. Razmotrimo sledeći dinamički sistem:

x[t] = 0.4 x[t-1] +0.2, y[t]= - 0.3 x[t-1] - 0.2 y[t-1] +0.3

x[0]=1,  y[0] =1.

Primetimo da evolucija x zavisi samo od x dok evolucija y zavisi i od x i od y. To bi značilo da možemo da prvo nađemo zavisnost x od vremena, onda zamenimo taj izraz u jednačinu za y, a onda opet rešimo i po y. Kao rezultat dobićemo rešenje sistema kao funkciju x i y od vremena.

Međutim, taj metod ne radi u opštem slučaju kada i x zavisi od x and y, a isto tako i y. Stoga, razmotrimo dva načina rešavanja sistema koji mogu da se primene u opštem slučaju. Oni će biti analogni situaciji kada smo imali diferencne jednačine sa samo jednom promenljivom. Ako uvedemo vektor vec[t] = {x[t], y[t]}, onda sistem možemo da predstavimo kao

A={{0.4,0},{-.3,-.2}};b={.2,.3};

{x[t],y[t]}==A.{x[t-1],y[t-1]}+b

Ovde je A matrica sistema, a b vektor sistema. Dati su u našem slučaju izrazima

Rešenje korišćenjem NestList

Naćimo prvih 20 vrednosti

resenje1=NestList[Dot[A,#]+b&,{1,1},20]

Pretstavimo prvo vremensku zavisnost x i y od vremena

ListLinePlot[resenje1\[Transpose],PlotRange->All,PlotLegends-> {"x[t]","y[t]"},LabelStyle->{16,Bold},ImageSize->500,AxesLabel->{"t","x,y"}]

A sada predstavimo i trajektoriju, odnosno kako y zavisi od x. To ćemo dobiti tako što prikažemo parove tačaka za x i y koje smo dobili rešavanjem sistema diferencnih jednačina.

ListLinePlot[resenje1,PlotRange->All,PlotLegends-> "y[x]",AxesLabel->{"x","y"}]

Rešenje korišćenjem dinamičkog programiranja

Definišimo sada dinamički program

x[t_]:=x[t]= 0.4 x[t-1] +0.2; y[t_]:= y[t]=- 0.3 x[t-1] - 0.2 y[t-1] +0.3;

x[0]=1;y[0]=1;

Sada možemo da generišemo rešenje pomoću Table, na primer

resenje2=Table[{x[t],y[t]},{t,0,20,1}]

Pretstavimo prvo vremensku zavisnost x i y od vremena

ListLinePlot[resenje2\[Transpose],PlotRange->All,PlotLegends-> {"x[t]","y[t]"},LabelStyle->{16,Bold},ImageSize->500,AxesLabel->{"t","x,y"}]

A sada predstavimo i trajektoriju, odnosno kako y zavisi od x. To ćemo dobiti tako što prikažemo parove tačaka za x i y koje smo dobili rešavanjem sistema diferencnih jednačina.

ListLinePlot[resenje2,PlotRange->All,PlotLegends-> "y[x]",AxesLabel->{"x","y"}]

Stanje dugoročne ravnoteže (fiksna tačka)

Postoji više načina kako možemo da nađemo fiksnu tačku dinamičkog sistema (ukoliko ona postoji, a sa grafika vidimo da postoji). Naime, jedan način je sledeći. Kada je fiksna tačka dostignusta, sistem se više ne menja, to jest vec[t-1] = vec]t]=vec. Stoga dobijamo da u takvoj tački

vec==A.vec+b

Ovaj sistem jednačina možemo da prepišemo i ovako: (I -A).vec =b, ili, u kodu, pošto je sistem dimenzije 2:

(IdentityMatrix[2]- A).vec ==b

Posle invertovanja dobijamo rešenje za fiksnu tačku

vec=Inverse[IdentityMatrix[2]-A].b

To isto možemo dobiti i korišćenjem FixedPoint funkcije u WL:

FixedPoint[A.#+b&,{1,1}]

Transformacija diferencne jednačine drugog reda u sistem diferencnih jednačina prvog reda

Ispostavlja se da se diferencne jednačine drugog i viših redova uvek mogu prepisati u ekvivalentan sistem diferencnih jednačina prvog reda. To se postići tako što ćemo uvoditi nove varijable umesto iste varijable sa lagovima. Hajde da vidimo to na primeru.

Cena nekog proizvoda zadovoljava sledeću diferencnu jednačinu drugog reda

(Subscript[p, t])==1.1-0.65(Subscript[p, t-1])-0.3 (Subscript[p, t-2]), (Subscript[p, 0])==1,(Subscript[p, 1])==0.7

Možemo je rešiti, na primer, korišćenjem dinamičkog programiranja za prvih 30 koraka

Clear[p];p[t_]:=p[t]=1.1-0.65 p[t-1]-0.3 p[t-2];p[0]=1;p[1]=0.7;

resenje3=Table[p[t],{t,0,30}]

Nastavak

Uvedimo sada varijablu (Subscript[x, t]) = (Subscript[p, t-1]) kao drugu jednačinu. A u prvoj jednačini, (Subscript[p, t-1]) ostavimo onako kako je, a umesto (Subscript[p, t-2]) koristimo (Subscript[x, t-1]). U tom slučaju dobijamo sledeći sistem jednačina.

(Subscript[p, t])==1.1-0.65(Subscript[p, t-1])-0.3 (Subscript[p, t-2])==1.1-0.65(Subscript[p, t-1])-0.3 (Subscript[x, t-1])

(Subscript[x, t])==(Subscript[p, t-1])

Da bismo sistem napravili u matričnom obliku, možemo da ga prepišemo malo detaljnije

(Subscript[p, t])==-0.65(Subscript[p, t-1])-0.3 (Subscript[x, t-1])+1.1

(Subscript[x, t])==1 (Subscript[p, t-1])+0 (Subscript[x, t-1])

Matrica sistema i vektor sistema su, respektivno

A={{-0.65,-0.3},{1,0}};b={1.1,0};

a rešenje sistema bi bilo

NestList[A.#+b&,{0.7,1},30]

resenje4=Last/@NestList[A.#+b&,{0.7,1},30]

Fiksna tačka

Primetimo da su u dugoročnoj ravnoteži, (Subscript[p, t]) i (Subscript[x, t]) isti. To ima smisla, s obzirom na definiciju varijable.

FixedPoint[A.#+b&,{0.7,1}]

Naravno, to bismo mogli dobiti i direktno preko izraza

Inverse[(IdentityMatrix[2]-A)].b

Međuzavisnost inflacije, inflacionih očekivanja i stope nezaposlenosti

Na osnovu ekonometrijskih studija poznato je da postoji međuzavisnosti između gornje tri veličine. Neka je stopa inflacije u trenutku t data sa (Subscript[p, t]), stopa očekivane inflacije u trenutku t data sa (Subscript[ep, t]), a stopa nezaposlenosti sa (Subscript[u, t]). Pretpostavimo, dalje, da važe sledeći odnosi

(Subscript[p, t])==0.13-3 (Subscript[u, t])+(Subscript[ep, t])

(Subscript[p, t+1])-(Subscript[p, t])==0.75 ((Subscript[p, t])-(Subscript[ep, t]))

(Subscript[u, t+1])-(Subscript[u, t])==-0.5(0.03-(Subscript[p, t+1]))

U prvoj jednačini vidimo da je trenutna inflacija obrnuto proporcionalna stopi nezaposlenosti. Drugi rečima, ukoliko stopa nezaposlenosti raste, proizvođači neće moći da dižu cene u nebesa zato što onda njihove proizvode neće kupovati. Sa druge strane, trenutna inflacije je direktno proporcionalna inflacionim očekivanjima.

Druge dve jednačine prikazuju dinamiku sistema. Vidimo da je promena inflacije proporcionalna razlici između realizovane i očekivane inflacije. Sa druge strane, promena stope nezaposlenosti zavisi.

Kako transformisati sistem

Pretpostavimo da nas zanima evolucija inflacije p, i nezaposlenosti u. U tom slučaju, iz prve jednačine možemo da zamenimo izraz za ep u drugu jednačinu. Druga i treća jednačina poslaju

(Subscript[p, t+1])-(Subscript[p, t])==0.75 ((Subscript[p, t])-(Subscript[ep, t]))==0.75 (0.13-3 (Subscript[u, t]))

(Subscript[u, t+1])-(Subscript[u, t])==-0.5(0.03-(Subscript[p, t+1]))

Sada smo skoro tu gde treba. Naime, ukoliko iz gornje jednačine zamenimo izraz za (Subscript[p, t+1]) u donju, dobijamo

(Subscript[p, t+1])-(Subscript[p, t])==0.75 (0.13-3 (Subscript[u, t]))

(Subscript[u, t+1])-(Subscript[u, t])==-0.5(0.03-((Subscript[p, t])+0.75 (0.13-3 (Subscript[u, t]))))

Nastavak

Sada možemo da predstavimo ove izraze kao standardni sistem dve diferencne jednačine gde lag varijable idu sa desne strane jednačine. Naime

(Subscript[p, t+1])==(Subscript[p, t])+0.75 (0.13-3 (Subscript[u, t]))//Expand

(Subscript[u, t+1])==(Subscript[u, t])-0.5(0.03-((Subscript[p, t])+0.75 (0.13-3 (Subscript[u, t]))))//Expand

Ukoliko uvedemo vektor vec čija je prva komponenta p a druga u, sistem se može predstaviti u standardnoj matričnoj formi gde je

A={{1,-2.25},{0.5,-0.125}};b={0.0975,0.03375};

Rešenje

Sada možemo da rešimo ovaj problem. Prvo ćemo naći p i u kao funkciju vremena. Pretpostavićemo da je početna inflacija 5 posto a stopa nezaposlenosti 11 posto.

resenje5=NestList[A.#+b&,{0.05,0.11},30]

ListLinePlot[resenje5\[Transpose],PlotRange->All,PlotLegends-> {"p[t]","u[t]"},LabelStyle->{16,Bold},ImageSize->500,AxesLabel->{"t","p,u"}]

A sada predstavimo i trajektoriju, odnosno kako u zavisi od p. To ćemo dobiti tako što prikažemo parove tačaka za x i y koje smo dobili rešavanjem sistema diferencnih jednačina.

ListLinePlot[resenje5,PlotRange->All,PlotLegends-> "u[p]",AxesLabel->{"p","u"}]

Sistemi običnih diferencijalnih jednačina

Kao što se jedna diferencna jednačine drugog ili višeg reda može predstaviti preko sistema diferencnih jednačina prvog reda, tako se i skalarna ODJ drugog ili višeg reda može predstaviti preko sistema diferencijalnih jednačina prvog reda.

Pretpostavimo da vozimo automobil brzinom od 90 km/h i da pred nama odjednom stane automobil na nekom rastojanju. Koliko to rastojanje mora biti da bismo uspeli da ukočimo na vreme. Poznato je da je usporednje a=13 metara u sekundi na kvadrat.

Jedan način da to rešimo jeste da postavimo diferencijalnu jednačinu

s''[t]=-13

Početni uslovi kažu da je početna brzina jednaka 90 kmh ili

v0=90/3.6

metara u sekunti

Dakle, da bismo našli pređeni put kao funkciju vremena, treba da rešimo jednačinu

DSolveValue[{s''[t]==-13 ,s[0]==0,s'[0]==25},s[t],{t}]//Expand

Nastavak

Da bismo saznali koje je minimalno rastojanje na kome bismo morali biti da ne udarimo u kola, potrebno je da nađemo trenutak t za koji je brzina postala jednaka nuli. To dobijamo tako što izjednačimo

Solve[D[25 t-((13 (t^2))/2),t]==0,t]

To znači da je najmanje rastojanje na kome bismo morali biti jednako

25 t-((13 (t^2))/2)/.{t->(25/13)}//N

to je oko 24 metara.

Rešimo sada ovaj problem malo drugačije. Naime, ukoliko uvedemo, pored s[t], kao nepoznatu funkciju i brzinu brz[t] koja je prvi izvod po vremenu od s, problem možemo predstaviti i kao sistem dve diferencijalne jednačine prvog reda sa dve nepoznate funkcije, s i brz.

DSolve[{brz'[t]==-13,s'[t]==brz[t],s[0]==0,brz[0]==25},{brz[t],s[t]},t]

Rešenje opet možemo dobiti preko DSolve. Znači, pored zavisnost s[t] odmah dobijamo i zavisnost brzine od vremena. Odatle odmah dobijamo traženi odgovor.

Kretanje u više dimenzija

Do sada smo razmatrali samo kretanje tela koje je moguće opisati preko pređenog puta, odnosno preko skalarne funkcije vremena. U realnosti tela se kreću u 3 dimenzije. Ali, da ne bismo komplikovali, princip rešavanja sistema diferencijalnih jednačina možemo da prikažemo i u dve dimenzije. Razmotrimo kretanje tela na koja ne deluju nikakve sile

Clear[x,y,z,g,v]

DSolve[{x''[t]==0,y''[t]==0,x[0]==x0,y[0]==y0,x'[0]==vx0,y'[0]==vy0},{x[t],y[t]}, t]

Napomenimo da je važno da prilikom ispisivanja sistema, prvo pišemo sve diferencijalne jednačine, onda početne uslova bez izvoda, pa početne uslove sa prvim izvodima itd. To je zgodno i radi toga da nešto ne zaboravimo. Jasno je da u ovom slučaju rešavamo svaku jednačinu pojedinačno pošto ne postoji veza između jednačine po x i po y. Time bi se, da radimo peške, problem sveo na rešavanje jedne diferencijalne jednačine, ali dva puta.

Kako izgleda putanja koja odgovara ovom kretanju? Pretpostavimo da imamo zadat vektor početne brzine {3,4} i da smo na početku u koordinatnom početku. U tom slučaju dobijamo sledeće rešenje:

DSolve[{x''[t]==0,y''[t]==0,x[0]==0,y[0]==0,x'[0]==3,y'[0]==4},{x[t],y[t]}, t]

Trajektorija ili putanja se može predstaviti pomoću ParametricPlot:

ParametricPlot[{3 t,4 t},{t,0,10}]

Dobijamo opet Prvi Newton-ov zakon. Naime, telo će se kretati pravolinijski ukoliko na njega ne deluje nikakva sila. Smer kretanja biće određen vektorom početne brzine.

Kosi hitac

Neka ubrzanje Zemljine teže deluje naniže duž y ose. Bacili smo loptu sa tačke (0, y0), t.j. sa visine y0. Rešenje se može razdvojiti na rešenje za x i na rešenje za y (pošto su same jednačine nezavisne jedna od druge, jer u jednačini za x nemamo y i obrnuto).

DSolve[{x''[t]==0,y''[t]==-g,x[0]==0,y[0]==y0,x'[0]==vx0,y'[0]==vy0},{x[t],y[t]}, t]//Expand

Vremenski trenutak u kome se hitac ponovo spusti na zemlju dat jednačinom:

t1=Solve[y0+vy0 t - g ((t^2)/2)==0,t][[2,1,2]]

Vreme povratka na tlo t1 zavisi samo od početne visine i od y komponente početne brzine (jer nema otpora vazduha). Sa druge stane, ono ne zavisi od početne brzine u horizontalnom pravcu. Pošto se maksimalna visina ostvaruje u trenutku:

Module[{y=y0+vy0 t - g ((t^2)/2)},Solve[D[y,t]==0,{t}]]

maksimalna visina je data izrazom:

ymax=y0+vy0 t - g ((t^2)/2) /. t-> vy0/g

Najzad, rastojanje koje će hitac preći po x osi je dato izrazom i zavisi od početne brzine po x osi:

xmax=vx0 t1;

Manipulate[g=9.81;ParametricPlot[{vx0 t, y0+vy0 t - g ((t^2)/2)},{t,0,T},AxesOrigin ->{0,0},PlotRange-> {{0,180},{0,100}},PlotStyle->Directive[Blue,Thick], PlotLabel->  Style[Framed["Kosi hitac"],16,Blue,Background->Lighter[Yellow]],ImageSize-> 500,
Epilog-> Inset[Grid[{{"Vreme leta",((Sqrt[(vy0^2)+2 g y0])+vy0)/g},{"Maks rastojanje",vx0* ((Sqrt[(vy0^2)+2 g y0])+vy0)/g},{"Maks visina",((vy0^2)/(2*g))+y0}},Frame-> All, Alignment-> Center, Spacings-> {2,1},ItemStyle->Directive[Bold,15]],{135,80}] ],{{vx0,10,"početna brzina vx0"},10,20,Appearance-> "Labeled"},{{vy0,10,"početna brzina vy0"},10,30,Appearance-> "Labeled"},{{y0,20,"početna visina y0"},20,50,Appearance-> "Labeled"},{{T,0.1,"vreme leta T"},0.1,7.5},SaveDefinitions->True]

Dinamički sistemi

Svi sistemi ODJ mogu biti, uvođenjem novih promenljivih ako je potrebno, predstavljeni kao sistemi diferencijalnih jednačina prvog reda Stoga možemo bez gubitka opštosti razmatrati samo sisteme ODJ prvog reda.

U ovom primeru jednačina za x zavisi od y kao i obrnutno. Stoga se sistem se ne može rešavati separatno za svaku promenljivu, nego samo zajedno. Srećom, zahvaljujući WL nama nije neophodno da znamo detalje kako se to radi peške (nije teško).

Clear[x,y,t]

x'[t]==2y[t]+10

y'[t]==- x[t]-y[t]-3

x[0]==1,y[0]==-1

Evo odmah konkretnog rešenja:

{resenjex,resenjey}=DSolveValue[{x'[t]==2y[t]+10,y'[t]==-x[t]-y[t]-3,x[0]==1,y[0]==1},{x[t], y[t]},t]//Expand

Primetimo da se ovo rešenje može grafički predstaviti ovako

ParametricPlot[{resenjex,resenjey},{t,0,30}]

Numeričko rešavanje običnih diferencijalnih jednačina

Mada DSolveValue odnosno DSolve može da reši veliki broj diferencijalnih jednačina egzaktno, najveći broj jednačina, posebno onih nelinearnih, moraju se rešavati numerički. Ovde ćemo prvo prikazati dva standardna metoda za numeričko rešavanje jednačina, a onda razmotriti i NDSolve, moćnu funkciju za numeričko rešavanje diferencijalnih jednačina.

Ispostavlja se da se numeričko rešavanje diferencijalnih jednačina defakto svodi na rešavanje diferencnih jednačina. Kada, pri tome, interpoliramo rešenje (to jest od nekih tačaka napravimo diferencijabilnu funkciju), takvo rešenje može biti i precizno i podjednako lako za korišćenje kao i egzaktno rešenje.

Objasnićemo prvo najstariji metod za numeričko rešavanja ODJ, takozvani Euler metod, a onda i Runge-Kutta koji se u praksi najčešće koristi. Najzad ćemo poraditi sa NDSolve

Euler-ov metod

Pretpostavimo da želimo da rešimo diferencijalnu jednačinu u intervalu između a i b (a<b) tako da:

Clear[a]

y'[t]==f[t,y[t]]
y[a]==(Subscript[y, 0])

Najjednostavniji način da to rešimo numerički jeste metod koji je smislio čuveni Leonard Euler (on je inače najveći deo života bio potpuno slep!). Prva observacija jeste da se ODJ prvog reda može ekvivalentno predstaviti u integralnoj formi. Počnimo od:

y'[t]==f[t,y[t]]

Ukoliko pomnožimo obe strane jednačine sa dt dobijamo identitet

dy==y'[t]dt==f[t,y[t]]dt

Sada prointegralimo od a do t obe strane:

y[t]-y[a]==(∫_a^t)f[u,y[u]]\[DifferentialD]u

Nastavak

Početnu diferencijalnu jednačinu smo, stoga, pretvorili u integralnu jednačinu  za nepoznatu funkciju y[t]:

y[t]==(Subscript[y, 0])+(∫_a^t)f[u,y[u]]\[DifferentialD]u

Konstruišimo sada približno rešenje za t ∈[a,b] koristeći tu integralnu formu. Naime, podelimo interval u kome traži rešenje na jednake intervale dužine Δt = (b-a)/n. Onda predstavimo približno integral preko sume

y[b]==(Subscript[y, 0])+(∫_a^b)f[u,y[u]]\[DifferentialD]u==(Subscript[y, 0])+(∑_(i=1)^n)(∫_(Subscript[t, i-1])^(Subscript[t, i]))f[u,y[u]]\[DifferentialD]u==(Subscript[y, 0])+(∑_(i=1)^n)f[(Subscript[t, i-1]),(Subscript[y, t-1])]Δt

Saglasno Euler-ovom metodu, integral po malo intervalu se aproksimira preko vrednosti funkcije u početnoj tački intervala puta dužina intervala Δt. Rešenje onda dobijamo iterativno, počevši od {(Subscript[t, 0]), (Subscript[y, 0])}

{(Subscript[t, i]),(Subscript[y, i])}=={(Subscript[t, i-1])+Δt,(Subscript[y, i-1])+f[(Subscript[t, i-1]),(Subscript[y, i-1])]Δt}

Rešenje je, znači, lista tačaka (parova)

{{(Subscript[t, 0]),(Subscript[y, 0])},{(Subscript[t, 1]),(Subscript[y, 1])},...,{(Subscript[t, n]),(Subscript[y, n])}}

Implementacija Euler-ovog metoda

Pretpostavimo da želimo da rešimo numerički diferencijalnu jednačinu (znači, f je y[t]+Exp[t])

y'[t]==y[t]+Exp[t]
y[0]==1

u intervalu između a =0 i b=0.2. Radi jednostavnosti, uzmimo samo n =4 iteracije. Evo parametara.

t0=a=0;y0=1;a=0.0;b=0.2;n=4;

Sada izračunamo korak integracije Δt

Δt=(b-a)/n

i definišemo funkciju f

f[t_,y_]:=y +Exp[t]

Ovako implementiramo metod Eulera-a:

NestList[{#[[1]]+Δt,#[[2]]+Δt f[#[[1]],#[[2]]]}&,{a,y0},n]

Metod Runge-Kutta

Ukoliko želimo profesionalnu preciznost, koristimo metod Runge-Kutta koji aproksimira integral u svakom malom intervalu preko metoda trapeza na sledeći način:

y1 = y0 + (k1 + 2 k2 + 2 k3 +k4)/6  gde je

k1 = Δt f[t0,y0]
 k2 = Δt f[t0+Δt/2, y0+k1/2]
 k3 = Δt f[t0+Δt/2, y0+k2/2]
 k4 = Δt f[t0+Δt, y0+k3]

Da bismo to implementirali, koristimo Module koji će računati k1, k2,k3, i k4.

NestList[{#[[1]]+Δt,Module[{k1,k2,k3,k4},k1=Δt f[#[[1]],#[[2]]];k2=Δt f[#[[1]]+Δt/2,#[[2]]+k1/2];k3=Δt f[#[[1]]+Δt/2,#[[2]]+k2/2];k4=Δt f[#[[1]]+Δt,#[[2]]+k3];#[[2]]+(1/6)*(k1+2 k2+2 k3+k4)]}&,{a,y0},n]

Egzaktno rešenje na istom gridu

Koje od ovih vrednosti su bolja aproksimacija tačnog rešenja? Da bismo na to odgovorili, treba da nađe egzatne vrednosti. Prvo rešimo diferencijalnu jednačinu preko DSolveValue

Clear[y]

DSolveValue[{y'[t]==y[t]+Exp[t],y[0]==1},y[t],t]

a onda izračunamo i rešenje za te vremenske tačke:

Table[{t,(E^t) (1+t)},{t,0,0.2,0.05}]

Drugim rečima, ne možemo da razlikujemo rešenje Runge-Kutta od tačnog rešenja ali je Euler-ov metod manje precizan.

Uporedimo sada rezultate kada je b veće i imamo veći broj koraka n

Pretpostavimo da smo sada povećali b pa je sada 2 (umesto 0.2). Početna tačka ostaje ista, kao i Δt.  U tom slučaju:

t0=a=0;y0=1;b=2.;n=40;Δt=(b-a)/n

euler=NestList[{#[[1]]+Δt,#[[2]]+Δt f[#[[1]],#[[2]]]}&,{a,y0},n]

rungekutta=NestList[{#[[1]]+Δt,Module[{k1,k2,k3,k4},k1=Δt f[#[[1]],#[[2]]];k2=Δt f[#[[1]]+Δt/2,#[[2]]+k1/2];k3=Δt f[#[[1]]+Δt/2,#[[2]]+k2/2];k4=Δt f[#[[1]]+Δt,#[[2]]+k3];#[[2]]+(1/6)*(k1+2 k2+2 k3+k4)]}&,{a,y0},n]

exact=Table[{t,(E^t) (1+t)},{t,0,2.0,0.05}]

Primetimo da Runge-Kutta metod daje izvrsne rezultate u ovom slučaju dok  Euler-ov metod greši oko 10 percent! Drugim rečima taj metod jeste prostiji, ali je i neprecizniji

ListLinePlot[{euler,rungekutta,exact},PlotStyle-> {{Red,Dashed},{Black,Thick,Dotted},LightGray},ImageSize->500,LabelStyle->{16,Bold},PlotLegends->{"euler","rungekutta","exact"}]

Interpolacija i interpolaciona funkcija

Svako numeričko rešavanje diferencijalih jednačina kao rezultat daje skup tačaka (vrednosti parametra kao što je, na primer, vreme, i procenjene funkcije-rešenja). Interpolacija pretvara skup tačaka u glatku funkciju sa kojom možemo da radimo isto kao i sa svakom čistom funkcijom. Razmotrimo skup tačaka koje daje Runge-Kutta metod

rungekutta

Napravimo sada interpolaciju pomoću funkcije Interpolation

frunge=Interpolation[rungekutta]

Funkcija frunge se koristiti kao bilo koja čista funkcija. Možemo da izračunamo vrednost funkcije u nekoj tački koja nije deo uzorka:

frunge[0.73]

Možemo je predstaviti funkciju grafički, kao i njen prvi izvod, na primer:

Plot[{frunge[t],frunge'[t]},{t,1,6},PlotLegends->"Expressions"]

Kako naći numeričko rešenje neke obične diferencijalne jednačine preko NDSolve?

Pođimo od veoma prostog primera za koji znamo egzaktno rešenje. Naime, hoćemo da koristeći NDSolve rešimo jednačinu y'(x)==y(x) za 0<x<2. Pretpostavljamo da je početna vrednost jednaka 1. Primetimo da je sintaksa skoro ista kao za DSolve, sa dve razlike:

NDSolve[{y'[t]==y[t],y[0]==1},y,{t,0,2}]

Naime, ne stavljamo argument t u y (pa dobijemo čistu funkciju) a kao drugo moramo na kraju umesto promenljive t da stavimo ceo interval u kome želimo da nađemo rešenje (obe sintakse takve vrste moguće su, doduše, i u DSolve).

Na ovaj način definišemo čistu funkciju interpoliranog rešenja:

resenje=NDSolveValue[{y'[t]==y[t],y[0]==1},y,{t,0,2}]

Sa njom radimo šta hoćemo. Recimo, možemo predstaviti rešenje grafički ovako:

Plot[resenje[t],{t,0,2}]

Možemo i da izračunamo vrednost prvog izvoda funkcije u tački 0.4, na primer

resenje'[0.4]

Numeričko rešenje jednostavnog sistema diferencijalnih jednačina

Ovako možemo da rešimo jednostavan sistem običnih diferencijalnih jednačina numerički

NDSolve[{y'[t]==z[t],z'[t]==-y[t],y[0]==0,z[0]==1},{y,z},{t,0,Pi}]

Na ovaj način definišemo čiste funkcije koje daju par rešenja:

{resenjey,resenjez}=NDSolveValue[{y'[t]==z[t],z'[t]==-y[t],y[0]==0,z[0]==1},{y,z},{t,0,Pi}]

Gledano pojedinačno, kao funkcije vremena, imamo:

Plot[{resenjey[t],resenjez[t]},{t,0,Pi},PlotLegends->{"y[t]","z[t]"}]

Prikažimo putanju (u y i z osama) grafički pomoću ParametricPlot:

ParametricPlot[{resenjey[t],resenjez[t]},{t,0,Pi}]

Prikažimo sada i prve izvode kao funkcije vremena:

Plot[{resenjey'[t],resenjez'[t]},{t,0,Pi},PlotLegends->{"y[t]","z[t]"}]

A ovo je takozvani vazni dijagram gde prikazujemo z’ kao funkciju y’:

ParametricPlot[{resenjey'[t],resenjez'[t]},{t,0,Pi}]

Malo složeniji primer

Razmotrimo sada numeričko rešenje za nelinearni sistem diferencijalnih jednačina:

NDSolve[{x'[t]==-y[t]-x[t]^2,y'[t]==2 x[t]-y[t],x[0]==y[0]==1},{x,y},{t,10}]

Kao i ranije možemo izvući definicije interpolacionih funkcija za x i y. Naime

Clear[resenjex,resenjey]

{resenjex,resenjey}=NDSolveValue[{x'[t]==-y[t]-x[t]^2,y'[t]==2 x[t]-y[t],x[0]==y[0]==1},{x,y},{t,10}]

Predstavimo rešenje za obe funkcije kako preko Plot tako i preko ParametricPlot

Plot[{resenjex[t],resenjey[t]},{t,0,10},PlotLegends->{"x[t]","y[t]"}]

A evo i putanje i faznog dijagrama

ParametricPlot[{resenjex[t],resenjey[t]},{t,0,10},AxesLabel->{"x","y"}]

ParametricPlot[{resenjex'[t],resenjey'[t]},{t,0,10},AxesLabel->{"x'","y'"}]

Interesantno bi bilo videti evoluciju, to jest kako se rešenje menja sa vremenom.

Manipulate[ParametricPlot[{resenjex[t],resenjey[t]},{t,0,τ},AxesLabel->{"x","y"},PlotRange-> 1.3],{τ,0.01,10}]

Lotka-Volterra jednačine za evoluciju broja  grabljivaca i plena

Ovaj sistem ODJ opisuje brojnost predatora i plena (recimo lisice i zečevi). Za pozitivne konstante A (brzina rasta broja zečeva), B (brzina sa kojom lisice uništavaju zečeve), C (stopa smrtnosti lisica), i D (stopa po kojoj broj lisica raste zato što jedu zečeve), u intervalu dt važi:

1. Populacija zečeva x raste proporcionalno broju zečeva za dx=Axdt ali i smanjuje propocionalno proizvodu broja zečeva i lisica za dx=-Bxydt.

2. Populacija lisica y opada proporcionalno broju lisica po stopi dy=-Cydt ali i raste po stopi dy=Dxydt. Stoga:

(dx)/(dt) =  Ax-Bxy	

(dy)/(dt) = -Cy+Dxy,

Neka je zečeva 3 puta više nego lisica na početku a odgovarajući parametri su dati na sledeći način. Ovako možemo da definišemo broj zečeva i lisica kao čiste funkcije:

{zecevi,lisice}=NDSolveValue[{y'[t]==y[t](x[t]-1),x'[t]==x[t](2-y[t]),x[0]==3,y[0]==1},{x,y},{t,0,10}]

Predstavimo rešenja prvo kao funkciju vremena, a posle i kao zavisnost broj lisica od broja zečeva:

Plot[{zecevi[t],lisice[t]},{t,0,10},PlotLegends->{"zecevi","lisice"}]

A evo i direktne zavisnosti brojeva zečeva i lisica.

ParametricPlot[{zecevi[t],lisice[t]},{t,0,10},AxesLabel->{"zecevi","lisice"}]

Proverimo kako se formira putanja

Manipulate[ParametricPlot[{zecevi[t],lisice[t]},{t,0,τ},AxesLabel->{"zecevi","lisice"},PlotRange-> 5.5],{τ,0.01,10}]

Upoređenje Euler and Runge-Kutta solutions; NDSolveValue

Uporedimo sada interpolisane Euler i Runge-Kutta rešenja koje smo ranije našli sa onim što daje NDSolve. Prethodno smo već našli interpolaciju za Runge-Kutta. Sada uradimo to isto za Euler

feuler=Interpolation[euler]

Evo i rešenja koje daje NDSolveValue

ndSolvesolution=NDSolveValue[{y'[t]==y[t]+Exp[t],y[0]==1},y,{t,0,2}];

Izračunajmo sada vrednost funkcije i njenog prvog izvoda u tački 1.024 za različite metod

{ndSolvesolution[1.024],ndSolvesolution'[1.024]}

{frunge[1.024],frunge'[1.024]}

{feuler[1.024],feuler'[1.024]}

{Exp[t](1+t),D[Exp[t](1+t),t]}/.{t->1.024}

Vidimo da NDSolveValue daje dobar rezultat za vrednost funkcije i za nijansu slabije za prvi izvod.

Malo razbibrige posle ovolikih jednačina

https://www.youtube.com/watch?v=cWGE9Gi0bB0