Predavanje 23

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]]}]

Diferencne i diferencijalne jednačine sa jednom promenljivom

Diferencne jednačine sa jednom promenljivom

To su jednačine koje daju rekurzivan način da se nađe kako vrednost neke veličine danas zavisi od prošlih vrednosti te iste veličine kao i, potencijalno jedne ili više poznatih promenljivih. Ukoliko se, pri tome, koristi lag od jednog koraka, u pitanju su jednačine prvog reda. Razmotrimo za početak upravo njih.

Pretpostavimo da smo oročili u banci 500 evra po kamati od 3 posto godišnje. Kamata se obračunava jednom godišnje. Naći evoluciju našeg računa u banci. Koliko bismo imali novca posle 5, 10, i 20 godina?

Ukoliko je količina novca na računu u godini t jednama (Subscript[b, t]) a u godini t-1 jednaka (Subscript[b, t-1]) onda je veza između te dve veličine jednaka (ukoliko nam na račun ne seda plata i ne plaćamo ništa)

(Subscript[b, t])==1.03 (Subscript[b, t-1]), (Subscript[b, 0])==500

Vidimo da je to diferencna jednačina prvog reda (ima samo jedan lag). Nađimo rešenje na više načina. Prvi je da definišemo funkciju b[t] rekurzivno

Clear[b];b[t_]:=b[t]=1.03b[t-1];b[0]=500;

Vrednosti posle 5, 10 i 20 godina bi bile jednake

{b[5],b[10],b[20]}

a u opštem slučaju

Clear[b];b[t_]:=b[t]=(1+r)b[t-1];b[0]=b0;

{b[5],b[10],b[20]}

Drugi način izračunavanja

Istu diferencnu jednačinu možemo da rešimo i koristeći Nest. Naime, ukoliko nas zanima da nađemo b[5] u numeričkom primeru, možemo to uraditi i ovako

Nest[1.03 #&, 500,5]

Na sličan način možemo dobiti sve tri tražene vrednosti

Table[Nest[1.03 #&, 500,n],{n,{5,10,20}}]

a u opštim brojevima

Table[Nest[(1+r) #&, b0,n],{n,{5,10,20}}]

Primetimo da dok god je r>0, odnosno bruto prinos 1+r veći od 1, iznos na našem računu će rasti bez ograničenja.  Naravno, ukoliko bi se naplaćivala, kao u Švajcarskoj, negativna kamata na štednju, to ne bi bio slučaj. Na primer, ukoliko je kamata negativna i jednaka -1 posto, posle 5, 10, ili 20 godina imali bismo na računu

Table[Nest[0.99 #&, 500,n],{n,{5,10,20}}]

Vremenom bi naš račun konvergirao nuli.

Evolucija cene jednog dobra

Pretpostavimo da su tražnja i ponuda neke robe linearne funkcije cene. Pri tome, tražnja raste od trenutne cene (u trenutku t) dok ponuda zavisi od cene u prethodnom periodu. Konkretno

(Subscript[Q, dt])==1-(Subscript[p, t]), (Subscript[Q, st])==0.1+0.65(Subscript[p, t-1])

Pretpostavimo, dalje, da je početna cena (u trenutku 0) jednaka 1. Naći ravnotežnu putanju cene u ovom slučaju.

U svakom trenutku ponuda je u ravnoteži jednaka tražnji, pa važi da je

1-(Subscript[p, t])==0.1+0.65(Subscript[p, t-1])

ili, posle male transformacije

(Subscript[p, t])==0.9-0.65 (Subscript[p, t-1]), (Subscript[p, 0])==1

Putanja stoga zadovoljava sledeću jednačinu

p[t_]:=p[t]=0.9 -0.65 p[t-1];p[0]=1;

Evolucija cene jednog dobra

Izračunajmo čemu bi bila jednaka cena posle 20 perioda

Table[p[t],{t,0,20}]

Ovo možemo predstaviti i grafički

ListLinePlot[Table[p[t],{t,0,20}],PlotRange-> All]

Evoluciju možemo izračunati i koristeći NestList. U tom slučaju:

NestList[0.9 -0.65 #&,1,20]

ListLinePlot[NestList[0.9 -0.65 #&,1,20],PlotRange-> All,LabelStyle->{16,Bold},ImageSize->500]

Primetimo da sistem konvergira ka fisnoj tački. U tom slučaju kažemo da je sistem dinamički stabilan, odnosno da postoji stabilna dugoročna dinamička ravnotežna cena.

Kako naći fiksnu tačku za ovu diferencnu jednačin

Sa grafika vidimo da se cena, posle oscilacija, stabilizuje oko stanja dugoročne ravnoteže. Ovo je fiksna tačna procesa cena. Fiksne tačke su od primarnog interesa u klasičnoj ekonomskoj teoriji. U WL postoji i posebna funkcija koja pomaže da se one nađu:

FixedPoint[0.9 -0.65 #&,1]

We can figure out the fixed point from the first principles as follows. Let us come back to the original difference equation:

p[t]=0.9 -0.65 p[t-1]

Subtracting p[t-1] from both sides of this equation, on the left hand side we obtain the change in price Δp = p[t] - p[t-1]. The right hand side becomes:

0.9-0.65p[t-1]-p[t-1]==0.9 - 1.65 p[t-1]

Thus, difference equation can be rewritten as: Δp = 0.9 - 1.65 p[t-1]. Price stops changing (i.e. the left hand side is equal to zero) when the right hand side is zero as well. This is the fixed point.

Clear[p];Solve[0.9-1.65 p[t-1]==0,p[t-1]]

Uvod u diferencijalne jednačine

Pojam diferencijalnih jednačina uveli su Newton i Leibnitz. Čuvena je Newton-ova izreka: “Ja znam diferencijalne jednačine koje opisuju kretanje tela. Dajte mi početne uslove za svako telo i moći ću da izračunam u svakom trenutku gde se koje nalazi.”

Diferencijalne jednačine opisuju dinamiku kako prirodnih tako i mnogih društvenih procesa. Uglavno se fokusiramo na takozvane obične diferencijalne jednačine (ODJ). One zavise od jednog parametra po kome se diferencira.

(((d^n)y)/(dt^n)) = f[y[t], (dy/dt),...,(((d^(n-1))y)/(dt^(n-1))),t]
{y,dy/dt, ...,(((d^(n-1))y)/(dt^(n-1)))}[t0]==const

Prvi red definiše klasu mogućih rešenja. Drugi red specifikuje o kome se rešenju tačno rešenju radi. Taj drugi red su početni uslovi i-ili granični uslovi. Funkcija y[t] je rešenje gornje diferencijalne jednačine ukoliko, kada je ubacimo u tu jednačinu, dobijamo algebarsku identičnost i ukoliko su zadovoljeni početni odnosno granični uslovi. Cilj rešavanja diferencijalnih jednačina jeste da na osnovu znanja izvoda neke funkcije nađemo funkcionalnu zavisnost za tu funkciju.

Red obične diferencijalne jednačine je najviši izvod koji se pojavljuje u diferencijalnoj jednačini. Gornja ODJ je n-tog reda.

Ukoliko razmatramo ODJ reda n, da bismo mogli da preciziramo njeno konkretno rešenje potrebno je da zadamo n nezavisnih uslova. Na primer, često su zadate početna vrednost funkcije, kao i početne vrednosti izvoda sve do reda n-1. Umesto jedne ODJ, dinamika sistema može biti opisana i sistemom ODE.

Ukoliko je nepoznata funkcija zadata preko izvoda više promenljivih, onda je to parcijalna diferencijalna jednačina (PDJ). U tom slučaju, rešenje je funkcija više promenljivih. Te promenljive mogu biti, na primer, vreme i jedna ili više prostornih koordinata. Na ovom kursu nećemo proučavati parcijalne diferencijalne jednačine. Stoga ćemo dalje izbaciti reč obična i govoriti samo o diferencijalnim jednačinama.

Diferencijalne jednačine prvog reda

Diferencijalne jednačine smo zapravo već i videli i rešavali, makar da to nismo tako zvali. Naime, kada smo pokušavali da rekonstruišemo kako pređeni put nekog tela zavisi od vremena na osnovu vremenske zavisnosti brzine tog tela, mi smo zapravimo rešavali jednačinu tipa

s'[t]==v[t]

gde je s[t] pređeni put kao funkcija vremena, a v[t] je poznata funkcija brzine kao funkcije vremena. U tom slučaju, nepoznata funkcija y je s[t]. Funkcija f u ovom slučaju zavisi samo od vremena ekplicitno, a ne preko nepoznate funkcije i njenih izvoda kao u opštem slučaju. To jest, f = v[t] je poznata funkcija. Pošto je

s'[t]==v[t]

pomnožimo sa obe strane sa dt i dobijamo da je

s'[t]\[DifferentialD]t==\[DifferentialD]s==v[t]\[DifferentialD]t

Rešenje je, stoga, u opštem slučaju, neodređeni integral

s[t]==∫v[t]\[DifferentialD]t

Takvo opšte rešenje zavisi od konstante integracije koja može biti bilo šta. Ukoliko, međutim, znamo početni položaj (na primer da je s[0] ==s0, rešenje bi bilo

s[t]==s0 +(∫_0^t)v[u]\[DifferentialD]u

Da vidimo kako to radi na primeru

Pretpostavimo da se vozimo konstantnom brzinom od 90 km na sat a da počnemo od tačke nula. Pošto počinjemo iz nule, konstanta integracije je jednaka 0. U tom slučaju, pređeni put bi bio

s[t]==∫90\[DifferentialD]t

što se dalo pogoditi unapred. Pretpostavimo sada malo kompleksniji problem istog tipa, a naime brzina nije konstantna nego nego funkcija, npr

s'[t]== 90 -0.7 t+0.5 (t^2)

Rešenje bi ponovo bilo integral

s[t]==Integrate[90 -0.7 t+0.5 (t^2),t]

Kako izgleda ta putanja?

DSolve i DSolveValue su funkcije za rešavanje diferencijalnih jednačina u egzaktnoj formi

Pretpostavimo da želimo da nađemo rešenje prethodno diferencijalne jednačine koristeći DSolve. Evo kako bismo to uradili

DSolve[{s'[t]==90 -0.7 t+ 0.5(t^2),s[0]==0},s[t],{t}]

Primetimo da smo na početku zadali listu sa diferencijalnom jednačinom i početnim uslovom. Sledeća stavka je nepoznata funkcija. Najzad, poslednja lista je lista varijabli. Na taj način dobijamo konkretno egzaktno rešenje diferencijalne jednačine. Kao što vidimo, dobijamo rešenje u obliku pravila. Ako želimo umesto toga da dobijemo rešenje u obliku tražene funkcije direktno, koristimo

DSolveValue[{s'[t]==90 -0.7 t+ 0.5(t^2),s[0]==0},s[t],{t}]

Tu smo pokazali kako da nađemo konkretno rešenje diferencijalne jednačine. Šta ukoliko želimo da nađemo opšte rešenje? U tom slučaju, sve što treba da uradimo jeste da izostavimo početni uslov

DSolveValue[{s'[t]==90 -0.7 t+ 0.5(t^2)},s[t],{t}]

Pretpostavimo da smo počeli da merimo pređeni put ne od nule, nego, recimo, od 150. U dom slučaju, sve što bi se promenilo jeste konstanta integracije, i mi bismo dobili

DSolveValue[{s'[t]==90 -0.7 t+ 0.5(t^2),s[0]==150},s[t],{t}]

Kako iz marginalnog naći ukupan trošak proizvodnje?

Evo i primera iz ekonomije gde parametar diferencijalne jednačine nije vreme nego količina proizvedene robe q. Fiksni trošak je c[0]== 12. Pretpostavimo da smo zaključili da je marginalni trošak proizvodnje količine robe q dat izrazom:

c'[q]==3 q +((q^3)/2)

Primetimo da je i u ovom slučaju, f funkcija samo parametra q (koji sada figurira umesto vremena). Da  bismo rešili jednačinu peške, možemo da pomnožimo obe strane sa dq i sagledamo da je

c'[q]dq==d c==(3 q +((q^3)/2))dq

ili, integrirajući obe strane

c[q]==Integrate[3 q +((q^3)/2),q]+12

Dakle, ukupni trоšak nije ništa drugo nego rešenje gornje diferencijalne jednačine, odnosno funkcija c[q].

Kako dobiti to rešenje preko DSolveValue

Nađimo prvo opšte rešenje:

DSolveValue[{c'[q]==3 q +((q^3)/2)},c[q],q]

Kao i obično, opšte rešenje zavisi od proizvoljne konstante. Nju fiksiramo na osnovu uslova za fiksni trošak. Naime:

DSolveValue[{c'[q]==3 q +((q^3)/2),c[0]==12},c[q],q]//Expand

Vidimo da je konstanta zapravo fiksni trošak, a ostatak je varijabilni troškak! Proverimo da su zadovoljena i ODE i početni uslov:

D[12+((3 (q^2))/2)+((q^4)/8),q]==3 q +((q^3)/2)

12+((3 (q^2))/2)+((q^4)/8)/.q->0

Novac uložen u banku ukoliko se kamata obračunava neprekidno

Obe prethodne jednačine su bile jednačine prvog reda najprostijeg tipa gde se rešenje moglo naći prostim integraljenjem obe strane jednačine. To je bilo stoga što brzina odnosno marginalni trošak nije zavisio od pređenog puta odnosno ukupnog troška. Razmotrimo sada situaciju u kojoj je f linearna funkcija nepoznate funkcije, naime jednačinu ukamaćenja novca na računu u neprekidnom vremenu.

Ukoliko se kamata obračunava neprekidno (godišnja neprekidna kamata je r) a novac ne podiže sa računa, akumulirani iznos na računu u banci u nekom vrlo malom intervalu vremene dt može da se predstavi kao

Clear[b]

b[t+dt]-b[t]==db ==r b[t] dt

Ukoliko podelimo obe strane sa b[t] (pretpostavimo da iznos nije jednak nuli u trenutku t) dobijamo da su varijale razdvojene

(db/b)==r dt

Ukoliko uspešno prepišemo diferencijalnu jednačinu tako da sa jedne strane imamo samo nepoznatu funkciju a sa druge strane parameter t, obe strane možemo, bar u principu, da prointegralimo i da dobijemo rešenje. Leva strana daje Log[b] dok desna strane daje r t. Stoga je

Log[b[t]]==r t +c1

Ovo opšte rešenje možemo da prepišemo ekponencirajući obe strane jednačine

b[t]==Exp[c1]Exp[r t]==c2 Exp[r t]

Ukoliko je u trenutku t=0 na račun stavnjeno b[0] = b0, vidimo da konstanta c2 nije proizvoljna nego jednaka b0, pa je konkretno rešenje onda

b[t]==b0 Exp[r t]

Kako to isto uraditi koristeći DSolveValue

Ukoliko jednačinu

db ==r b[t] dt

podelimo sa obe strane sa dt, dobijamo da je

b'[t]==r b[t]

Opšte i konkretno rešenje dobijamo sasvim jednostavno

DSolveValue[{b'[t]==r b[t]},b[t],t]

DSolveValue[{b'[t]==r b[t],b[0]==b0},b[t],t]

Početna suma će ekponencijalno rasti ukoliko ne podižemo novac sa računa pod uslovom da je kamatna stopa r pozitivna, i opadati ukoliko je kamatna stopa negativna. Proverimo da konkretno rešenje zadovoljava diferencijalnu jednačinu. Ukoliko ga ubacimo u nju umesto proizvoljne funkcije b[t] treba da dobijemo algebarski identičan izraz

D[b0 (E^(r t)),t]==r b0 (E^(r t))

Slično proveravamo da to rešenje zadovoljava početni uslov

b0 (E^(r t))==b0/.{t->0}

b[t_]:=b0 (E^(r t))

Clear[b]

Jednačina radioaktivnosti

Evo još jedne jednačine istog tipa kao prethodna, samo što oznako malo drugačije.

Na Generalštab u Nemanjinoj ulici 1999. je bačen 1 kg obogaćenog urana. Ekipe za radiološku dekontaminaciju očistile su 90 posto urana pa je ostalo 100 gr. Ukoliko je poznato da je period poluraspada 4.468×(10^9) godina, naći koliko će radioaktivnog urana na tom mestu biti kroz 100 godina.

Neka je n[t] količina urana u trenutku t, a T je period poluraspada, t.j. vreme potrebno da se količina radioaktivne materije smanji na pola. Brzina smanjenja količine radioaktine materije proporcionalna je preostaloj količini urana pa n[t] zadovoljava sledeću ODJ prvog reda:

n'[t]==((-Log[2])/T)n[t]

Opšte rešenje je ekponencijalnog oblika:

DSolveValue[{n'[t]==((-Log[2])/T)n[t]},n[t],t]

Pošto je ODJ prvog reda, da bi našli jedinstveno rešenje treba da fiksiramo jedan početni uslov. Neka je količina urana u trenutku 0 je jednaka n0:

DSolveValue[{n'[t]==((-Log[2])/T)n[t],n[0]==n0},n[t],t]

Vremenom, količina radioaktivnog urana će se smanjivati. Ali, pitanje je koliko brzo će se to desiti. Ubacimo podatke:

(2^(-(t/T))) n0/.{n0-> 100,T-> 4.468×(10^9),t-> 100}

Na žalost, količina radioaktivnog urana za 100 godina ostaće praktično ista kao danas.

Opštiji slučaj linearne ODJ prvog reda: naći pređeni put ukoliko brzina zavisi i od pređenog puta

Pretpostavimo da je brzina vozila data izrazom (opet, radi jednostavnosti, počinjemo put da merimo od trenutka 0)

s'[t]==90-0.1 s[t]

To bi značilo da brzina ima konstantnu komponentu od 90 (recimo km na sat) ali i da opada kako se udaljavamo od početka puta (recimo, nailazimo na sve veću gužvu). To je primer diferencijalne jednačine u kome je f sada linearna funkcija nepoznate funkcije. Takve jednačine se zovu linearne diferencijalne jednačine. One se, kao pravilo, manje više uvek mogu egzaktno rešiti. Pošto u ovom slučaju desna strana jednačine ne zavisi ekplicitno od vremena (nego samo od nepoznate funkcije), možemo ponovo razvoditi promenljive i dobiti

ds ==(90-0.1 s[t])dt

(ds/(90-0.1 s))==dt

ili, ako pomnožimo sa -1 obe strane

(ds/(0.1 s-90))==-dt

Naći pređeni put ukoliko brzina zavisi i od pređenog puta

Obe strane se sada mogu prointegraliti (probajte da vidite šta ćete dobiti ako to uradite!) i onda srediti tako da dobijemo s kao funkciju vremena. Umesto toga, možemo da rešimo problem koristeći DSolveValue. Ako nas zanima odmah konkretno rešenje, imaćemo da

DSolveValue[{s'[t]==90-0.1 s[t],s[0]==0},s[t],{t}]//Expand

Pretstavimo rešenje grafički

Plot[900.` -900.` (E^(-0.1` t)),{t,0,50}]

Primetimo da rešenje konvergira ta tački 900. Drugim rečima, nećemo se, ukoliko je ovo jednačina koja reguliše brzinu udaljiti više od 900 kilometara od početne tačke! To možemo da sagledamo odmah iz oblika same jednačine. Naime, 900 liči da je fiksna tačka funkcije s[t]. Fiksna tačka diferencijalne jednačine je ona u kojoj se s ne menja sa vremenom, odnosno s’[t] je jednaka nuli.

s'[t]==0==90-0.1 s[t]

Iz čega sledi da je fiksna tačka upravo 90/0.1 = 900.

Naći pređeni put ukoliko brzina zavisi i od pređenog puta

U opštijem slučaju, ne možemo da razdvojimo promenljive za vreme i put, odnosu t i y[t]. Predpostavimo, na primer, da umesto konstante 90, izraz za brzinu izgleda ovako (opet je radi jednostavnosti početni pređeni put jednak 0)

s'[t]== 90 -0.7 t+0.5 (t^2)-0.1 s[t]

Prilično je jasno da prethodni trik ne bi radio. Međutim, pošto je u pitanju linearna jednačina postoji način da se i ovako nešto reši. Ako postoji interesovanje pokazaću vam na tabli kako se tako nešto radi. Međutim, nama će biti dovoljno da to rešimo koristeći DSolveValue.

DSolveValue[{s'[t]== 90 -0.7 t+0.5 (t^2)-0.1 s[t],s[0]==0},s[t],{t}]//Expand

Evo kako izgleda ovo rešenje. Predstavićemo ga preko Manipulate

Manipulate[Plot[1969.9999999999998` -1969.9999999999998` (E^(-0.1` t))-107.` t+5.` (t^2),{t,0,τ}],{τ,0.1,50}]

Leva strana jednačine

D[1969.9999999999998` -1969.9999999999998` (E^(-0.1` t))-107.` t+5.` (t^2),t]

Desna strana jednačine

90 -0.7 t+0.5 (t^2)-0.1 (1969.9999999999998` -1969.9999999999998` (E^(-0.1` t))-107.` t+5.` (t^2))//Expand

Nelinearne obične diferencijalne jednačine prvog reda

Nelinearne ODJ često nemaju egzaktno rešenje. Riccati jednačina ga ima (kako bismo tu rastavili promenljive?). Ona služi za modeliranje dinamike rasta stanovništva, širenja epidemije ili glasina:

dy/dt == r y(1-(y/k))//TraditionalForm

Ovde je r udeo obolelih a 1/k udeo obolelih koji se izleči u jedinici vremena. Egzaktno rešenje je:

epidemija=DSolve[{y'[t]==r y[t]-(r/k)(((y[t]))^2),y[0]==1},y[t],t]//Quiet

Vremenom broj obolelih teži ka fiksnoj tački. U njoj, prvi izvod mora postati jednak nuli:

Solve[r y[t]-(r/k)(((y[t]))^2)==0,{y[t]}][[2]]

Limit[((k Exp[r t])/(k+Exp[r t]-1)),t-> ∞, Assumptions->r>0]

Nivo udela obolelih na dugi rok je određen samo time koliki procenat ljudi se izleči u jedinici vremena! Ali zanima nas i dinamika:

Clear[T]

Manipulate[Plot[((k (E^(r t)))/(k+(E^(r t))-1)),{t,0,τ}, AxesOrigin-> {0,0},PlotRange-> {0,5},PlotLabel->  Style[Framed["Model epidemije"],14,Blue,Background->Lighter[Yellow]], Epilog-> Inset[Grid[{{"Fiksna tacka",k}},Frame-> All, Alignment-> Center, Spacings-> {2,1},ItemStyle->Directive[Bold,15]],{Center,Center}]],{{r,0.01,"Brzina infekcije r"},0.01,0.4,Appearance-> "Labeled"},{{k,2,"Vrednost k"},2,5,Appearance-> "Labeled"},{{τ,0.1,"Vremenski period τ"},0.1,1000, Appearance-> "Labeled"} ,SaveDefinitions-> True]

Obične diferencijalne jednačine drugog reda: Newton-ove jednačine kretanja

Kao što sam rekao na početku, diferencijalne jednačine su bile izmišljene da bi opisivale kretanje tela. Prisetimo se drugog Newton-ovog zakona:

m a = F.

Ovde je a ubrzanje (to jest promena brzine u jedinici vremena), m je masa tela, a F je sila koja deluje na telo. Na ovom predavanju zamišljamo da živimo u jednodimenzionom sistemu, pa su sve ove veličine skalarne. Zanima nas pređeni put u svakom trenutku t, ukoliko znamo početni položaj i početnu brzinu tela.

Primetimo da je

a = v’[t] == s''[t]

Dakle, ubrzanje jeste drugi izvod pređenog puta po vremenu. Na osnovu prethodnih jednačina, imamo da pređeni put, saglasno Njutnovom drugom zakonu, zadovoljava ODJ drugog reda

s''[t]==((F[t,s[t],s'[t]])/m)

Ovde smo u opštem slučaju pretpostavili da sila može da zavisi kako ekplicitno od vremena, tako i od brzine. Pošto treba da rešimo jednačinu drugog reda, postojaće dve konstante integracije. Stoga, da bismo potpuno definisali rešenje moramo da imamo dva dodatna uslova. Tipično to su početni uslovi na položaj i brzinu s[0]=s0, i v[0]=v0. Intuitivno, jasno je da ako ne znamo početnu brzinu nego samo položaj, ne znamo ni da li ćemo sve vreme stajati u mestu ili otpočeti kretanje. Tako da samo početni položaj nije dovoljan.

Prvi Newton-ov zakon kao posledica drugog Newton-ovog zakona

Ukoliko je F =0, drugi Newton-ov zakon se svodi na diferencijalnu jednačinu m  s’’[t] = 0 ili, pošto je masa tela m>0, na jednačinu:

s''[t]==0

Primetimo da ukoliko je ubrzanje jednako 0, brzina mora biti konstantna. To se dobija prostom integracijom obe strane gornje jednačine

s''[t]dt==0 dt

iz čega sledi

s'[t]==c1

Ukoliko je brzina konstantna, onda to znači i da je pređeni put afina (linearna plus možda konstanta) funkcija vremena

s'[t]dt==c1 dt

ili, ponovo integraleći

s[t]==c1 t + c2

Prvi Newton-ov zakon kao posledica drugog Newton-ovog zakona

Opšte rešenje jednačine drugog reda sadrži 2 proizvoljne konstante.:

DSolveValue[{s''[t]==0},s[t],t]

Vidimo da je rešenje ove jednačine linearna funkcija po vremenu. Da bismo fiksirali putanju tela potrebno su dva uslova, npr, početni položaj x0 i početna brzina v0:

DSolveValue[{s''[t]==0,s[0]==s0,s'[0]==v0},x[t],t]

Primetimo da položaj tela linearno zavisi od vremena, kao i od početne brzine. Brzina tela na koje ne deluje nikakva sila će, stoga, uvek biti konstantna i jednaka početnoj brzini jer je, po definiciji, brzina prvi izvod vremenske zavisnosti pomeraja s[t] po vremenu:

D[t v0+s0,t]

Kada na telo ne deluju nikakve sile, ukoliko u početku telo miruje, t.j. ako je v0=0, ostaće zauvek u stanju mirovanja. Ako ima neku početnu brzinu v0 zauvek će se kretati tom brzinom. To je Prvi Newton-ov zakon, ili zakon inercije. Vidimo da on sledi iz jednačine kretanja, t.j. iz drugog Newton-ovog zakona, kada je spoljna sila koja deluje na telo jednaka nuli.

Slučaj konstantne sile (npr gravitacija)

Na telo deluje sila f, koja je konstantna. Naći zakon kretanja i rešiti diferencijalnu jednačinu. Zakon kretanja je:

m s''[t]==f

Ovu jednačinu možemo prepisati i kao jednačinu za ubrzanje:

s''[t]==(f/m)

Nađimo opšte rešenje:

DSolveValue[{s''[t]==(f/m)},s[t],t]//Expand

U slučaju konstantne sile f, rešenje se sastoji iz 2 dela. Prvi je rešenje homogene diferencijalne jednačine (kada je sila jednaka nuli), dok je drugi deo rešenje nehomogene jednačine.

Slučaj konstantne sile (npr gravitacija)

Pretpostavimo da nas zanima kretanje tela bačenog sa Zemlje vertikalno uvis brzinom v0. Njegov početni položaj je 0. Pošto je sila Zemljine teže f = - m g, jednačina kretanja postaje:

DSolveValue[{s''[t]==-((g m)/m),s[0]==0,s'[0]==v0},s[t],t]//Expand

Primetimo da položaj tela sada više ne opisuje linearna nego kvadratna funkcija vremena. Maksimalna visinu je ona u kojoj brzina postaje jednaka nuli:

Module[{s=-((g (t^2))/2)+t v0},Solve[D[s,t]==0,{t}]]

smax=-((g (t^2))/2)+t v0/.{t->(v0/g)}

Nastavak

Maksimalna visina ne zavisi od mase bačenog tela, ali zavisi od brzine sa kojom je telo bačeno! Sa kojom brzinom će se telo vratiti na zemlju? Diferencijalna jednačina je ista kao i pre, ali su početni uslovi sada drugačiji:

Clear[s]

DSolveValue[{s''[t]==-((g m)/m),s[0]==((v0^2)/(2 g)),s'[0]==0},s[t],t]//Expand

Da bismo našli vrednost brzine, ili prvog izvoda, u trenutku povratka na zemlju, prvo treba da nađemo vreme povratka. A to je trenutak kada visina postaje jednaka 0:

Solve[-((g (t^2))/2)+((v0^2)/(2 g))==0,t]

Pošto prvo rešenje nema smisla, dobijamo da je vreme potrebno da telo padne sa te visine jednako (v0/g). Brzina

s=-((g (t^2))/2)+((v0^2)/(2 g));

brzina sa kojom će udariti u zemlju je:

D[s,t]/.{t->(v0/g)}

Dakle, telo će pasti sa istom onom brzinom sa kojom je inicijalno bačeno u vis (samo suprotnog smera!).

Clear[s]

Spoljna periodička sila

Kako se kreće telo pod dejstvom periodičke F =m f0 Sin[ω t], šta bi onda bila jednačina kretanja? Pretpostavljamo da je inicijalno y[0]=0, kao i da je početna brzina v0. Prvo nadjimo opšte rešenje:

DSolveValue[{y''[t]==((m f0 Sin[ω t])/m)},y[t],t]//Expand

Vidimo da je rešenje opet rešenje homogene jednačine (prva dva člana) plus rešenje nehomogene jednačine (treći član). Zaista, ukoliko bismo našli drugi izvod po vremenu funkcije

D[-((f0 Sin[t ω])/( (ω^2))),{t,2}]

to se poklapa sa desnom stranom jednačine. Stoga je -((f0 Sin[t ω])/( (ω^2))) rešenje nehomogene jednačine. Opšti iskaz glasi: opšte rešenje nehomogene linearne ODJ bilo kog reda uvek je dato kao opšte rešenje homogene jednačine plus konkretno rešenje nehomogene jednačine. Sada možemo da odredimo konkretno rešenje gornje jednačine. Naime:

DSolveValue[{s''[t]==f0 Sin[ω t],s[0]==0,s'[0]==v0},s[t],t]//Expand

Kako izgleda pređeni put kao funkcija vremena? Definišimo neke parametre i pogledajmo kako izgleda dinamika u ovom slučaju. Onda možemo da menjamo parametre.

Module[{v0=0.1,f0=1,ω=1},Plot[t v0+((f0 t)/ω)-((f0 Sin[t ω])/( (ω^2))),{t,0,10}]]

Vidimo, dakle, da je rešenje kombinacija pravolinijskog kretanja i oscilacija indukovanih spoljnom silom.

Slučaj kada je sila proporcionalna rastojanju: periodičko kretanje

Razmotirmo situaciju u kojoj je sila data izrazom - 3 s[t]. U tom slučaju, što je veće rastojanje za koje se telo pomerilo, veća je sila koja se trudi da ga vrati ka početnom položaju (to vidimo po znaku za silu). Na primer, neka je

s''[t]==-3 s[t]

Njeno opšte rešenje bi bilo:

DSolveValue[{s''[t]+3 s[t]==0},s[t],t]

Dakle, ovo je jednačina za periodičko kretanje. Ukoliko je početna pozicija 10 a brzina 2:

sbeztrenja=DSolveValue[{s''[t]+3 s[t]==0,s[0]==10,s'[0]==2},s[t],t]//Expand

Grafički, rešenje izgleda ovako:

Plot[sbeztrenja,{t,0,10}]

Nastavak

Pretpostavimo sada da pored periodičke sile proporcionalne s[t] deluje i sila trenja koja smanjuje brzinu kretanja odnosno amplitudu oscilacija. Uticaj trenja može se predstaviti preko člana u jednačini proporcionalnog prvom izvodu y. Ako je taj član - (1/10) s'[t], dobili bismo sledeću jednačinu:

s''[t]== -3 s[t]-(1/10) s'[t]

Bez problema nalazimo konkretno rešenje za iste početne uslove kao ranije.

ssatrenjem=DSolveValue[{s''[t]+(1/10) s'[t]+3 s[t]==0,s[0]==10,s'[0]==2},s[t],t]//Expand

Na grafikonu upoređujemo rešenje sa i bez trenja, za iste početne uslove.

Plot[{sbeztrenja,ssatrenjem},{t,0,30},PlotLegends->{"Bez trenja","Sa trenjem"}]

Nastavak

Kao što smo i pretpostavljali, amplituda rešenja sa trenjem se vremenom smanjuje, dok u slučaju kada nema trenja ona, razume se, ostaje ista. Dodajmo sada na to, takođe i periodičku spoljnu silu F=2 Sin[t]. U tom slučaju jednačina kretanja postaje:

s''[t]==-3 s[t]-(1/10) s'[t]+2 Sin[t]

čije je rešenje, za iste početne uslove, jednako:

sspoljna=DSolveValue[{s''[t]+(1/10) s'[t]+3 s[t]==2 Sin[t],s[0]==10,s'[0]==2},s[t],t]//Expand

Vidimo da se rešenja postepeno komplikuju, ali ih DSolve uspešno pronalazi! Sačuvajmo ovo poslednje rešenje takođe. Najzad grafički predstavimo sva tri rešenja za dva različita intervala vremena.

Plot[{sbeztrenja,ssatrenjem,sspoljna},{t,0,30},PlotLegends->{"Bez trenja","Sa trenjem","spoljna sila"}]

Plot[{sbeztrenja,ssatrenjem,sspoljna},{t,30,60},PlotLegends->{"Bez trenja","Sa trenjem","spoljna sila"}]

Primetimo da treće rešenje vremenom počinje da dominira nad drugim rešenjem. To je stoga što spoljna sila dodaje energiju u oscilacije koje se inače polako gube.

Malo razbibrige posle ovolikih jednačina

https://vimeo.com/67907543