Predavanje 23: Primene, Glava 5

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

5. Diferencijalne jednačine u WL

1. Obične diferencijalne jednačine prvog reda

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.

Nastavak

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 bar dve promenljive. Te promenljive mogu biti, na primer, vreme i jedna ili više prostornih koordinata.

U ovom predavanju ćemo se suočiti sa oba tipa diferencijalnih jednačina. Nećemo morati da se mučimo oko tehnike rešavanja jednačina bilo da su u pitanju obične ili parcijalne diferencijalne jednačine metodom korak po korak. Jer, WL ima funkcije koje to sa velikom lakoćom rade za nas.

Najprostija obična diferencijalna jednačina prvog reda

Najprostija ODJ je sledeća:

y'[t]==0

Koja je to funkcija čiji je prvi izvod uvek jednak nuli? Odgovor je: proizvoljna konstanta. Naime, prvi izvod konstante je jednak nuli. U WL postoji moćna funkcija DSolve koja nalazi nepoznate funkcije na osnovu diferencijalne jednačine ili sistema diferencijalnih jednačina:

DSolve[{y'[t]==0},{y[t]},{t}]

Ovako nalazimo opšte rešenje. Prvi argument je lista diferencijalnih jednačina koje treba da rešimo. Drugi argument je lista funkcija koje treba da odredimo, a treći je lista promenljivih. Kada se lista nepoznatih funkcija i argumenata sastoji samo od jednog elementa, možemo pisati i ovako:

DSolve[{y'[t]==0},y[t],t]

DSolve daje rešenje u vidu Rule. Pošto je ovo ODJ prvog reda opšte rešenje sadrži jednu proizvoljnu konstantu.  Fiksiramo rešenje pretpostavivši da je u trenutku 0 funkcija y jednaka 3. Pošto je y konstanta, ona će i u svim ostalim trenucima biti jednaka 3. Proverimo:

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

Stoga možemo da kažemo da je y[t]==3 rešenje diferencijalne jednačine y’[t]==0 uz početni uslov y[0]==3. Umesto kao listu pravila rešenje možemo dobiti i kao izraz koristeći DSolveValue:

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

DSolveValue[{y'[t]==0,y[0]==3},y[t],t]

Jednačina radioaktivnosti

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]==a n[t]

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.

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

Ukoliko se kamata obračunava neprekidno a novac ne podiže sa računa, suma koja se akumulira u banci sledi sličnu ODJ. Neka je trenutku t iznos b[t]. Ukoliko je godišnja kamata r, između t i t+dt imali bismo dodatni prihod od db =r b dt. Podelimo obe strane sa dt:

db ==r b dt

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

Opšte rešenje je opet eksponencijalno:

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

Konstantu određujemo na osnovu početnog ulaganja. Neka je to iznos b0. Posle t godina na računu bismo imali:

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

b0 (E^(r t))/.{b0-> 1000,r-> 0.03,t-> 4}

Početna suma će ekponencijalno rasti ukoliko ne podižemo novac sa računa. Brzina rasta data je godišnjom kamatnom stopom. Proverimo da je ova eksponencijalna funkcija zaista rešenje gornje jednačine, kao i da zadovoljava početni uslov.

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

Kada ubacimo u diferencijalnu jednačinu gornji izraz, vidimo da je leva strana jednačine jednaka desnoj.

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

U trenutku 0 vidimo da je:

b[0]

Stoga rešenje zadovoljava i diferencijalnu jednačinu i početni uslov.

Clear[b]

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. Pretpostavimo da smo zaključili da je marginalni trošak proizvodnje količine robe q dat izrazom:

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

Na osnovu marginalnog nаći ukupni trošak proizvodnje c[q] (fiksni trošak c[0]== 12). Ukupni trоšak nije ništa drugo nego rešenje gornje diferencijalne jednačine, odnosno funkcija c[q]. Nađimo prvo opšte rešenje:

DSolveValue[{c'[q]==3 q +((q^3)/2),c[0]==12},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

Opštiji slučaj linearne ODJ prvog reda

ODJ je linearna jednačina prvog reda ukoliko se može predstaviti na sledeći način (ovde su a, b i g funkcije koje zavise samo od parametra po kome se diferencira, t.j. od t):

a[t]y'[t]+b[t] y[t]==g[t]

Ukoliko je g[t] ==0, jednačina se zove homogena. Ukoliko to nije slučaj, onda je nehomogena. Razmotrimo homogenu linearnu ODJ prvog reda (a==1, b==3, g==0):

y'[t]+3 y[t]==0

Rešenje ove jednačine je ekponenta, slično onome što smo imali ranije.

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

Uporedimo to rešenje sa onim za nehomogenu jednačinu, na primer kada je g[t] ==5 t - 7.

D[-(26/9)+((5 t)/3),t]+3 (-(26/9)+((5 t)/3))//Simplify

-(26/9)+((5 t)/3)+(E^(-3 t)) C[1]==8/.{t-> 0}

y[0]==8

y'[t]+3 y[t]==5 t-7

DSolve[{y'[t]+3 y[t]==5 t-7},y[t],t]

Opšte rešenje je, sada, suma dva dela. Prvi deo je rešenje homogene jednačine (kad je desna strana bila jednaka nuli). To opšte rešenje zavisiti od proizvoljne konstante koju ćemo fiksirati kasnije iz dodatnog uslova. A kakav je smisao ostatka rešenja? U pitanju je zapravo konkretno rešenje nehomogene jednačine.  Naime, zamislimo da je rešenje samo:

y[t_]:=-(26/9)+((5 t)/3)

Rešenje linearne nehomogene jednačine: opšte rešenje homogene jednačine plus konkretno rešenje nehomogene

Proverimo da je ovo rešenje kompletne jednačine:

y'[t]+3 y[t]

Simplify[%]

U slučaju nehomogene linearne ODJ, rešenje će uvek biti opšte rešenje homogene jednačine plus neko rešenje nehomogene jednačine. Pretpostavimo da nismo promenili homogenu jednačinu ali jesmo nehomogeni deo. Kakvo će sad biti opšte rešenje?

D[-(8/9)+((5 t)/3)+(t^2),t]+3(-(8/9)+((5 t)/3)+(t^2))//Simplify

Clear[y]

DSolveValue[{y'[t]+3 y[t]==3 (t^2)+7t-1},y[t],t]

Primetimo da je rešenje homogenog dela ostalo nepromenjeno i samo se promenilo rešenje nehomogene jednačine pošto se sada promenila desna strana. Posle ove digresije, vratimo se na prethodni primer.

Pošto je u pitanju jednačina prvog reda, dovoljno je da zadamo jedan uslov da potpuno fiksiramo rešenje. Na primer, možemo da zadamo vrednost funkcije u trenutku 0, y[0]==8.

Clear[y]

DSolveValue[{y'[t]+3 y[t]==5 t-7,y[0]==8},y[t],t]//Expand

Malo komplikovanija linearna jednačina ODE prvog reda

Šta se desi ako neki koeficijent u ODJ zavisi od vremena, npr ako umesto a=3 imamo a= 3 t (a sve ostalo je kao malopre):

y'[t]+3 t y[t]==5 t-7

Potražimo opšte rešenje:

DSolveValue[{y'[t]+3 t y[t]==5 t-7},y[t],t]//Expand

?Erfi

Ono zavisi od funkcije Erfi. Funkcija Erfi je definisana u WL preko sledećeg integrala koji je tabulliran:

  Erf[z]=(2/(Sqrt[π]))(∫_0^z)(e^(-(t^2)))dt

yres=(5/3)+(19/3) (E^(-((3 (t^2))/2)))-7 (E^(-((3 (t^2))/2))) (Sqrt[(π/6)]) Erfi[(Sqrt[(3/2)]) t]

D[yres,t]+3 t (yres)==5 t -7//Simplify

To su primer specijalne funkcije. One su definisane preko integrala ili beskonačnih redova. Uporedimo grafički rešenja jednačina kada je a=3 i kada je a=3 t:

Plot[{-(26/9)+((98 (E^(-3 t)))/9)+((5 t)/3),(5/3)+(19/3) (E^(-((3 (t^2))/2)))-7 (E^(-((3 (t^2))/2))) (Sqrt[(π/6)]) Erfi[(Sqrt[(3/2)]) t]},{t,0,4},PlotLegends->"Expressions"]

Primetimo da za dovoljnu veliku vrednost t, prvo rešenje rasle linearno sa t, dok za drugo dobijamo

Limit[(5/3)+(19/3) (E^(-((3 (t^2))/2)))-7 (E^(-((3 (t^2))/2))) (Sqrt[(π/6)]) Erfi[(Sqrt[(3/2)]) t],t-> ∞]

Stoga će prvo rešenje će vremenom dominirati u odnosu na drugo.

Nelinearne obične diferencijalne jednačine prvog reda

Linearne ODJ često nemaju egzaktno rešenje. Riccati jednačina ga ima. 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]

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]

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

Uvod

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 (tipično, u pitanju je vektor), m je masa tela, a F je sila (takođe vektor, istih dimenzija kao ubrzanje). Za početak zamislimo da živimo u jednodimenzionom prostoru. Telo stoga ima samo jednu koordinatu, recimo x, koja zavisi od vremena. Onda je kretanje opisano kao funkcija x[t]. Drugi Newton-ov zakon je zapravo diferencijalna jednačina drugog reda, pošto je ubrzanje drugi izvod položaja po vremenu.

a = x''[t]

Sila F takođe može da zavisi od vremena. Ali za početak, razmotrimo najprostiji primer, naime slučaj kada na telo ne deluje nikakva sila.

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  x’’[t] = 0 ili, pošto je masa tela m>0, na jednačinu:

x''[t]==0

Opšte rešenje jednačine drugog reda sadrži 2 proizvoljne konstante. U našem slučaju:

DSolveValue[{x''[t]==0},x[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[{x''[t]==0,x[0]==x0,x'[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 x[t] po vremenu:

D[t v0+x0,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)

Neka se telo i dalje kreće samo u jednoj dimenziji, duž z ose. Kretanje tela je opisano funkcijom  z[t]. Na telo deluje sila f, koja je konstantna. Sila je usmerena duž z ose. Naći zakon kretanja i rešiti diferencijalnu jednačinu. Zakon kretanja je:

m z''[t]==f

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

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

Nađimo opšte rešenje:

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

U slučaju konstantne sile f, rešenje se sastoji iz 2 dela. Prvi je rešenje homogene diferencijalne jednačine, dok je drugi deo rešenje nehomogene jednačine. 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[{z''[t]==-((g m)/m),z[0]==0,z'[0]==v0},z[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[{z=-((g (t^2))/2)+t v0},Solve[D[z,t]==0,{t}]]

zmax=-((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! Ovaj iskaz važi, kako ćemo posle razumeti, samo zato što ne uzimamo u obzir otpor vazduha. 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[z]

DSolveValue[{z''[t]==-((g m)/m),z[0]==((v0^2)/(2 g)),z'[0]==0},z[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

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

brzina sa kojom će udariti u zemlju je:

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

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

Spoljna periodička sila

Telo se kreće duž y ose. Ukoliko je F =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:

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

DSolveValue[{y''[t]==((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 ω])/(m (ω^2))),{t,2}]

to se poklapa sa desnom stranom jednačine. Stoga je -((f0 Sin[t ω])/(m (ω^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[{y''[t]==((f0 Sin[ω t])/m),y[0]==0,y'[0]==v0},y[t],t]//Expand

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,m=1,ω=5},Plot[t v0+((f0 t)/(m ω))-((f0 Sin[t ω])/(m (ω^2))),{t,0,10}]]

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

Opštije linearne jednačine drugog reda

Ranije smo razmatrali jednačine tipa y’[t] + 3 y[t]==0. Kako bi izgledalo rešenje jednačine ovakvog tipa:

D[Exp[a t],{t,2}]+3 Exp[a t]==0

c1 Exp[(Sqrt[3])I t]+c2 Exp[-(Sqrt[3])I t]

Solve[{(a^2)+3==0},{a}]

D[(t^3),{t,2}]

y''[t]+3 y[t]==0

Njeno opšte rešenje bi bilo:

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

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

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

Grafički, rešenje izgleda ovako:

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

Ukoliko ima trenja, amplituda oscilacija opada sa vremenom. Uticaj trenja može se predstaviti preko člana u jednačini proporcionalnog prvom izvodu y. Ako je taj član (1/10) y'[t], dobili bismo sledeću jednačinu:

y''[t]+(1/10) y'[t]+3 y[t]==0

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

Nastavak

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

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

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

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:

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

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

yspoljna=DSolveValue[{y''[t]+(1/10) y'[t]+3 y[t]==2 Sin[t],y[0]==10,y'[0]==2},y[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[{ybeztrenja,ysatrenjem,yspoljna},{t,0,30},PlotLegends->{"Bez trenja","Sa trenjem","spoljna sila"}]

Plot[{ybeztrenja,ysatrenjem,yspoljna},{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.

3. Slobodni pad kada postoji otpor vazduha

Pad bez otpora vazduha

Pretpostavimo da smo skočili sa 1000 metara visine. Bez otpora vazduha, jednačina kretanja bi, u opštem slučaju, bila:

DSolveValue[{m y''[t]==-m g,y[0]==y0,y'[0]==0},y[t],t]//Expand

Vreme leta bismo mogli da izračunamo ovako:

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

a brzina sa kojom bismo udarili u zemlju bila bi data izrazom:

D[(1/2) (-g (t^2)+2 y0),t]/.t->(((Sqrt[2]) (Sqrt[y0]))/(Sqrt[g]))

To znači da bismo, ako padamo sa 1000 metara visine udarili o tlo brzinom od oko 140 metara u sekundi

(Sqrt[2]) (Sqrt[g]) (Sqrt[y0])/.{g-> 9.81,y0-> 1000}

Jednačina u slučaju kada postoji otpor vazduha

Šta bi se desilo ako postoji otpor vazduha? U tom slučaju jednačina kretanja bi bila:

m v'[t]==-m g+(1/2)ρ A ((v[t])^2)

Sila trenja usporava pad i proporcionalna je gustini vazduha (ρ), površini tela koje se pada (A) kao i kvadratu brzine sa kojom se telo kreće. Zbog postojanja ovog drugog člana sa desne strane diferencijalne jednačine, to jest zbog sile trenja, prilikom leta uspostaviće se (i to relativno brzo) konstantna brzina pri kojoj je desna strana jednačine kretanja postati jednaka nuli. Ta brzina se na engleskom zove Terminal Velocity:

vterminal=(Sqrt[((2 m g)/(ρ A))])

Brzina sa kojom će se telo kretati i sa kojom će udariti o tlo sada potpuno nezavisna od visine sa koje se pada! Takođe, za razliku od situacije sa standardnim slobodnim padom, ta će brzina sada biti proporcionalna našoj težini a obrnuto proporcionalna našoj površini. Upravo zbog ovog drugog padobrani rade dobro! Jednačinu kretanja rešavamo ovako:

DSolve[{m y''[t]==-m g+(1/2)ρ A ((y'[t])^2),y[0]==y0,y'[0]==0},y[t],t][[1]]//Quiet

y=((A y0 ρ-2 m Log[Cosh[(((Sqrt[A]) (Sqrt[g]) t (Sqrt[ρ]))/((Sqrt[2]) (Sqrt[m])))]])/(A ρ))

v=D[((A y0 ρ-2 m Log[Cosh[(((Sqrt[A]) (Sqrt[g]) t (Sqrt[ρ]))/((Sqrt[2]) (Sqrt[m])))]])/(A ρ)),t]

Rešenje je komplikovanije nego za slobodni pad ali se da izračunati egzaktno. Da bismo našli kako brzina zavisi od vremena, sve što treba da uradimo jeste da prodiferenciramo prethodni izraz:

Pad sa padobranom

Padamo sa padobranom. Gustina vazduha na nivou mora je ρ->1.225 kg/(m^3). Naravno, na visini od 1000 metara gustina vazduha bi bila manja, ali ćemo time radi jednostavnosti prenebregnuti. Pretpostavimo da je površina padobrana A->1.5 π (2.5^2) (m^2), kao i da je masa čoveka 80 kg. Za masu padobrana našao sam podatak da je negde oko 14 kg. U tom slučaju, brzina padobranca u trenutku t bi bila data jednačinom:

Sin[x]=((Exp[I x]-Exp[- I x])/(2 I))

Sinh[x]=((Exp[ x]-Exp[- x])/(2 ))

?Tanh

((Exp[x]-Exp[-x])/(Exp[x]+Exp[-x]))

vpadobran=v/.{A->1.5 π (2.5^2),ρ->1.225 ,g->9.81,m->( 80+14)}

Plot[vpadobran,{t,0,5},AxesLabel->{"Vreme (s)","Brzina (m/s)"},PlotLabel-> Style["Brzina pada sa padobranom",14,Bold]]

Vidimo da bi se veoma brzo brzina padobranca stabilizovala na oko 7 (m/s).

Koliko bismo dugo padali sa padobranom? To određujemo na osnovu uslova da je visina postala jednaka 0.

y==0/.{A->1.5 π (2.5^2),ρ->1.225 ,g->9.81,m->( 80+14),y0-> 1000}

Solve[%,{t}][[2]]//Quiet

Vidimo da bismo padali nešto manje od 2.5 minuta.

Šta ako se padobran ne otvori?

Razmotrimo sada neprijatnu situaciju da se padobran nije otvorio. I dalje bi postojao otpor vazduha. Ono što je sada drugačije jeste površina čoveka, koja bi bila, recimo, 1.3 (m^2). Sada bi brzina pada bilo skoro pet puta veća nega sa padobranom. Ali, priroda rešenja bi ostala ista:

vsazatvorenimpadobranom=v/.{A->1.3,ρ->1.225 ,g->9.81,m-> ( 80+14)}

Plot[vsazatvorenimpadobranom,{t,0,10},AxesLabel->{"Vreme (s)","Brzina (m/s)"},PlotLabel-> Style["Brzina pada sa zatvorenim padobranom",14,Bold]]

Koliko dugo bismo padali bez padobrana? To dobijamo opet iz uslova da visina postaje jednaka 0.

y==0/.{A->1.3,ρ->1.225 ,g->9.81,m->(80+14),y0-> 1000}

Solve[%,{t}][[2]]//Quiet

Padali bismo oko pola minuta bez padobrana. Pri tome, došli bismo do konstantne brzine negde posle 7-8 sekundi.

4. Kretanje u više dimenzija: sistemi običnih diferencijalnih jednačina

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.

Kretanje bez dejstva 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. 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}]

5. Numeričko rešavanje običnih diferencijalnih jednačin korišćenjem NDSolve

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 sledeći skup tačaka:

tacke={Range[6],{1,2,3,5,8,5}}\[Transpose]

Napravimo sada interpolaciju pomoću funkcije Interpolation

f=Interpolation[tacke]

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

f[4]

f[2.5]

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

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

A ovako možemo uporediti funkciju sa tačkama koje smo interpolirali:

Show[Plot[f[t],{t,1,6}],ListPlot[tacke,PlotStyle->{Red,PointSize[0.02]}]]

Možemo da izračunamo i njen integral od 1 do 6 (pošto je interpolirana funkcija, po definiciji, neprekidna).

Integrate[f[t],{t,1,6}]

Interpolacija omogućava da umesto da radimo sa nekim diskretnim skupom tačaka radimo sa funkcijom.

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:

Exp[1.2]

resenje[1.2]

Plot[{Exp[t],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:

Clear[resenjey]

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

6. Parcijalne diferencijalne jednačine prvog reda

Parcijalne diferencijalne jednačine (PDJ) su diferencijalne jednačine u kojima nepoznata funkcija zavisi najmanje od 2 promenljive po kojima se diferencira. Ukoliko su najviši izvodi prvog reda, takve jednačine su PDJ prvog reda. Evo ekonomskog primera.

Ako proizvodnja zavisi od dva faktora, na osnovu marginalnih troškova naći ukupni trošak

Pretpostavimo da su izmerili marginalne troškove po faktorima proizvodnje 1 i 2 te znamo da su dati kao

D[c[q1,q2],{q1}]==3 q1+2 q2;

D[c[q1,q2],{q2}]==2 q1+4 q2

Pronaći opšti izraz za c kao funkciju q1 i q2. U pitanju je sistem parcijalnih diferencijalnih jednačina. Zašto parcijalnih? Pa zato što nepoznata funkcija zavisi od dve promenljive, q1 i q2. Opšte rešenje je:

DSolveValue[{D[c[q1,q2],{q1}]==3 q1+2 q2,D[c[q1,q2],{q2}]==2 q1+4 q2},c[q1,q2],{q1,q2}]

Nije teško videti da bismo diferenciranjem funkcije ukupnog troška po q1 i po q2 dobili upravo ono što treba:

D[((3 (q1^2))/2)+2 q1 q2+2 (q2^2)+C[1],{q1}]

D[((3 (q1^2))/2)+2 q1 q2+2 (q2^2)+C[1],{q2}]

Dakle, to je rešenje ovog sistema linearnih parcijalnih diferencijalnih jednačina. Ukoliko znamo da je fiksni trošak jednak 12, što je dato uslovom c[0,0]==12, dobili bismo i konkretno rešenje za funkciju ukupnog troška.

((3 (q1^2))/2)+2 q1 q2+2 (q2^2)+12

7. Parcijalne diferencijalne jednačine drugog reda: jednačine provodjenja toplote

Uvod

PDJ drugog reda leže u osnovi većine inženjerskih modela kao i modela za vrednovanje finansijskih derivata. U finansijama se pojavljuju uglavnom jednačine u kojima postoji prvi izvod po vremenu i drugi izvod po prostornoj koordinati (takozvane jednačine provodjenja toplote). Stoga ćemo onda ćemo razmatrati kako se temperatura neke tanke niti menja sa vremenom t u zavisnosti od položaja na niti x.

Vremensko-prostorna zavisnost temperature beskonačne niti

Imamo nit koja nije ograničena ni na jednom od dva kraja. U početnom trenutku zagrevamo sistem tako da je temperatura u sredini niti najveća a da veoma brzo opada prema krajevima. Drugim rečima, zagrejali smo samo središte. Ta početna temperaturna raspodela duž x ose naziva se početni uslov (inicial condition ili ic). Osnovna jednačina prenosa temperature kaže da prvi izvod po vremenu treba da bude jednak ili proporcionalan drugom izvodu po prostornoj koordinati x (u ovom slučaju jednak).

heqn=D[u[x,t],t]==D[u[x,t],{x,2}];

ic=u[x, 0] == E^(-x^2);

ures=DSolveValue[{heqn,ic },u[x,t],{x,t}]

Nastavak primera sa beskonačnom niti

Pošto je nit beskonačno duga, ne moramo eksplicitno da zadajemo granične uslove osim da temperatura na plus i minus beskonačnosti mora biti jednaka nula (a to program podrazumeva). U početku raspodela temperature je bila gausijan. Ona je ostala gausijanska, ali se polako vremenom širi (širina gausijana (1+4 t)). To znači da postoji difuzija toplote iz centra ka krajevima.Proverimo sada da je izraz koji smo dobili zaista rešenje PDJ. Leva i desna strana jednačine su, respektivno

ures

D[ures,t]

D[ures,{x,2}]//Expand

Vidimo da se poklapaju. Takođe, početni uslov je zadovoljen zato što dobijamo da je:

ures/.{t-> 0}

Plot3D[ures,{x,-5,5},{t,0,4},PlotRange->All]

Razmotrimo sada kako se rešenje menja sa vremenom:

Manipulate[Plot[((E^(-((x^2)/(1+4 τ))))/(Sqrt[1+4 τ])),{x,-5,5},PlotRange-> 1],{τ,0.01,4}]

Vremensko-prostorna zavisnost temperature polu-beskonačne niti

Razmotrimo sličan problem samo što sada nit ima jedan fiksiran kraj (u tački 0) i jedan beskonačan kraj (u tački x-> ∞). PDJ je ista kao i pre, kao i početni uslov. Neka je granični uslov (boundary condition ili bc) u tački x=0 da je temperatura konstantna u vremenu na granici i jednaka 1. Primetimo da početni uslov mora biti konsistentan sa graničnim. Pošto je gausijan u tački x=0 jednak 1, konsistentnost je zadovoljena.

heqn=D[u[x,t],t]==D[u[x,t],{x,2}];

ic=u[x, 0] == E^(-x^2); bc=u[0,t]==1;

Ovako možemo da predstavimo rešenje (pomoću Refine definišemo da nas zanimaju samo pozitivne vrednosti x):

ures=Refine[DSolveValue[{heqn,ic,bc },u[x,t],{x,t}],x>0]//FullSimplify

Plot3D[ures,{x,0,15},{t,0,10},Exclusions->None,PlotRange->All]

Nastavak

Proverimo da je to zaista rešenje PDJ

D[ures,t]//Expand

D[ures,{x,2}]//Expand

Vidimo da se poklapaju. Takođe, početni uslov je zadovoljen zato što, kada t->0, Erf u oba dela izraza teži 1.

Razmotrimo sada evoluciju:

Manipulate[Plot[(((E^(-((x^2)/(1+4 t)))) Erf[(x/(2 (Sqrt[t (1+4 t)])))])/(Sqrt[1+4 t]))+Erfc[(x/(2 (Sqrt[t])))],{x,0.01,15},PlotRange-> 1],{t,0.01,10}]

Vremensko-prostorna zavisnost temperature niti ograničene sa obe strane

Ukoliko je nit ograničena sa obe strane, moramo zadati dva granična uslova u odnosu na x (pošto imamo izvode drugog reda po x) kao i jedan početni uslov u odnosu na t (pošto imamo izvod prvog reda po t). Opet, početni i granični uslovi moraju biti konsistentni.

Razmotrimo nit koja ide od x=0 do x=2. Na obe granice je temperatura niti jednaka 0 (uslovi bc dole). Sa druge strane, početni uslov je zadat pomoću funkcije ic=u[x, 0] ==x (2-x)^2. U tačkama 0 i 2, ova funkcija je jednaka 0, pa to znači da su granični uslovi konsistentni sa početnim. Kada imamo ograničenje sa obe strane, rešenje jednačine biti dato, u principu, preko beskonačnog reda (zašto i kako, o tome neki drugi put). U pitanju je takozvani Fourier red.

heqn=D[u[x,t],t]==D[u[x,t],{x,2}];

ic=u[x, 0] ==x(2-x)^2;bc={u[0,t]==0,u[2,t]==0};

sol=DSolve[{heqn,ic,bc },u[x,t],{x,t}]

Vidimo da je rešenje dato u vidu sinusnih funkcija odgovarajućih frekvencija čije amplitude opadaju sa vremenom (član sa eksponentom). Umesto beskonačnog broja članova u praksi je dovoljno da uzmemo nekoliko prvih (ovde 3 člana):

ures=u[x,t]/.sol[[1]]/.{∞->3}//Activate

Nastavak sa konačnom niti

Proverimo da ovo rešenje zadovoljava diferencijalnu jednačinu:

D[ures,t]//Expand

D[ures,{x,2}]//Expand

Vidimo da se leva i desna strana poklapaju.

Plot3D[ures//Evaluate,{x,0,2},{t,0,3/2},Exclusions->None,PlotRange->All]

Primetimo da se profil vremenom razmazuje i da dolazi tako do difuzije. Ukoliko želimo da dodamo još članova, samo ćemo zameniti i staviti, recimo 7 umesto 3

asol=u[x,t]/.sol[[1]]/.{∞->7}//Activate

Plot3D[asol,{t,0,3/2},{x,0,2},Exclusions->None,PlotRange->All]

Kao što vidimo rešavanje diferencijalnih jednačina, kako običnih, tako i parcijalnih, sa WL nije uopšte teško.