Równania różniczkowe zwyczajne

Krótkie wprowadzenie do równań różniczkowych zwyczajnych

Przypomnienie: równaniem różniczkowym zwyczajnym nazywamy równanie (układ równań można potraktować jako równanie wektorowe) które zawiera pochodne funkcji jednej zmiennej. Np:

rownanie = y''[x]+y[x]==y'[x]+Sin[x]+2

uklad={y'[t]+x[t]==Sin[t],x'[t]+y[t]==Cos[t]}

Niewiadome są funkcje których poszukujemy ! Rozwiązanie zwykle zawiera tyle stałych dowolnych (,,stałych całkowania'') ile jest funkcji razy rząd najwyższej pochodnej w równaniu. W Mathematice rozwiązujemy równania (i tak samo układy równań) instrukcją DSolve . Posiada ona trzy argumenty: 
1. równanie lub układ równań wraz z ewentualnymi warunkami dodatkowymi (zwykle początkowymi lub brzegowymi)
2. niewiadome funkcje
3. zmienną po której różniczkujemy - z definicji r.r.zwyczajnego jest tylko jedna taka zmienna!

rozw1 = DSolve[rownanie, y[x],x] (* 2 stałe dowolne C[1] i C[2] ponieważ mamy jedną funkcję a najwyższa pochodna to druga pochodna *)

rozw2 = DSolve[uklad,{x[t],y[t]},t] (* 2 stałe dowolne C[1] i C[2] bo mamy 2 funkcje niewiadome obydwie najwyżej w pierwszej pochodnej *)

Rozwiazaniem ogólnym r. rózniczkowego jest zbiór na ogół nieskończenie wielu funkcji np. dla równań powyżej. Rozwiązaniem szczególnym nazywamy pewne konkretne rozwiązanie, np:

Manipulate[Plot[ (y[x] /. rozw1[[1]]) /. {C[1]->a, C[2]->b},{x,-10,10}],{a,-1,1},{b,-1,1}]

Manipulate[Plot[ Table[(y[x] /. rozw1[[1]]) /. { C[1]->a,C[2]->b},{b,-2,2,1}],{x,-5,5},PlotPoints->5],{a,-1,1}]

Manipulate[ParametricPlot[ ({x[t],y[t]} /. rozw2[[1]]) /. {C[1]->a, C[2]->b},{t,-2,2}],{a,-1,1},{b,-1,1}]

Znaczenie równań różniczkowych zwyczajnych

FIZYKA:
1. Równania ruchu Newtona
2  Równanie Schrodingera dla cząstki w potencjale 1D
3. Rozpad promieniotwórczy
4. Reakcje jadrowe w gwiazdach
5. Odkształcenie pretów
6. Równowaga ciekłych kryształów
INNE przykłady:
1. równania geodezyjnych w zakrzywionej czasoprzestrzeni
2. teoria chaosu
3. populacja drapiezników/roslinożerców
4. linie sił pół wektorowych

Istnieje gigantyczny zbiór REALNYCH problemów opisywanych równaniami różniczkowymi zwyczajnymi. Jest tu cały czas otwarte pole do nowych odkryć, lepszego zrozumienia problemów i rewizji  ,,starych'' metod. Co więcej, moc obliczeniowa współczesnych komputerów pozwala operować rozwiązaniami, w tym numerycznymi, równie łatwo jak fukcją Sin czy x^2.

Zagadnienie początkowe i brzegowe

O ile na rozwiązanie w postaci symbolicznej możemy nałożyć dowolne warunki dodatkowe, typowo przyjmują one postać warunków początkowych lub warunków brzegowych. Numeryczna wersja DSolve czyli NDSolve w większości wypadków także wymaga okreslenia warunków początkowych lub brzegowych. W innych wypadkach numeryczne rozwiązanie będzie znacznie utrudnione!

WARUNKI POCZĄTKOWE
Zadajemy wartości funkcji i jej pochodnych w jednym ustalonym punkcie!

DSolve[{rownanie,y[0]==1,y'[0]==2},y[x],x] (* 2 warunki zadane w punkcie x=0 *)

NDSolve[{rownanie,y[0]==1,y'[0]==2},y[x],{x,-1,1}] (* to samo numerycznie, zwróć uwagę na podany explicite zakres zmiennosci t *)

WARUNKI BRZEGOWE
Zadajemy wartości funkcji (i jej pochodnych jeżeli rząd równania (pochodnych) jest większy niż 2) na końcach przedziału:

DSolve[{rownanie,y[-1]==1,y[1]==2},y[x],x] (* warunki zadane w 2 punktach x=-1 oraz x=1 *)

NDSolve[{rownanie,y[-1]==1,y[1]==2},y[x],{x,-2,2}] (* numerycznie, przedział czasem może ,,wystawać'' poza punkty graniczne *)

WARUNKI  INNEGO TYPU

f1= y[x]/.DSolve[{rownanie,y[-1]-y'[1]==1,y[-1]+y'[1]==0},y[x],x] [[1]];

f1purefunction= y/.DSolve[{rownanie,y[-1]-y'[1]==1,y[-1]+y'[1]==0},y,x] [[1]];

Warunki wyglądają na spełnione:

{y[-1]-y'[1]==1,y[-1]+y'[1]==0} /. y->f1purefunction

rownanie /.  y->f1purefunction

Simplify[%]

Z takimi warunkami NDSolve na ogół nie radzi sobie!

f2 = y[x] /.NDSolve[{rownanie,y[-1]+y'[1]==1,y[-1]+y'[1]==0},y[x],{x,-2,2},Method->"BSD"][[1]]

Plot[{f1},{x,-2,2}]

Typowe wywołania DSolve

Sposób 1: podajemy w drugim argumencie DSolve funkcję y[s] wraz z argumentem s. Wynikiem jest wzór opisujący rozwiązanie:

DSolve[y''[s]==y[s],y[s],s]

wyrazenie=y[s] /. %[[1]]

Sposób 2: podajemy w drugim argumencie DSolve funkcję y bez argumentu s. Wynikiem jest funkcja opisująca rozwiązanie:

DSolve[y''[s]==y[s],y,s]

funkcja = y /. %[[1]]

Jeżeli chcemy sprawdzić czy otrzymalismy poprawne rozwiązanie, podstawiamy do równania funkcję a nie jej wzór:

y''[s]==y[s] /. y->funkcja (* DOBRZE *)

y''[s]==y[s] /. y[s]->(E^s) C[1]+(E^(-s)) C[2]  (* ŹLE !!! *)

y''[s]==y[s] /. y->(Exp[#]C[1]+Exp[-#]C[2] &) (* DOBRZE *)

Sposób 3: ParametricNDSolve ( od wersji 9.0 ! )

Sposób 4:  NDSolveValue ( od wersji 9.0 ! )

Powrót na poziom równań fundamentalnych (w fizyce)

Dysponując możliwościami oferowanymi przez CAS (Mathematica) operowanie gotowymi wzorami opisującymi rozmaite zjawiska fizyczne i analizowanie ,,typowych'' rozwiązań nie ma już więcej sensu. Należy wychodzić od równań fundamentalnych, w szczególnie jeżeli są to r. r. zwyczajne.

PRZYKŁAD: rzut w polu grawitacyjnym:
1. Arystoteles twierdził, ze ciało rzucone pod katem α najpierw porusza sie po linii prostej, aż wytraci „pęd”, potem po łuku okregu, a nastepnie pionowo w dół 
2. Każdy dobry uczen na to samo pytanie odpowie, że ciało porusza sie po paraboli i potrafi nawet podac na to odpowiedni wzór
PYTANIE: kto ma racje, najwiekszy filozof starozytnosci czy uczeń w
gimnazjum na początku XXI wieku?

Przykład: Rzut ukośny w polu grawitacyjnym

Podejście tradycyjne: Wyciagamy z pamięci wzór:

s = -g t^2/2 + v0 t

s /. {g->9.81, t->1.0, v0->10.0}

Podejście z użyciem CAS: Zaczynamy od II prawa Newtona

newton2=F==m a

rownanie = newton2 /.{F->-m g, a ->z''[t]}

DSolve[{rownanie,z[0]==0,z'[0]==v0},z[t],t]

z[t] /. %[[1]] /. {g->9.81, t->1.0, v0->10.0}

Pytanie dodatkowe: a jeżeli uwzględnimy opór powietrza ?

Podejście tradycyjne:

Tego nie bylo na wykladzie!

Podejście z użyciem CAS:

Ale było prawo Newtona i wzór na siłę oporu powietrza!

rownanie2 = newton2 /. {F->-m g -k z'[t], a->z''[t]}

rozw2 =DSolve[{rownanie2,z[0]==0,z'[0]==v0},z[t],t]

z[t] /. rozw2[[1]] /. {g->9.81, t->1.0, v0->10.0}

Ooops! Wynik zależy od masy ciała i współczynnika oporu! Niech będzie m=0.5 kg i k=0.1:

z[t] /. rozw2[[1]] /. {g->9.81, t->1.0, v0->10.0, m->0.5, k->0.1}

Odpowiedź została w obydwu przypadkach znaleziona równie łatwo.

Podejście tradycyjne:

To niesprawiedliwe. Opór powietrza rośnie z kwadratem prędkości i maleje z wysokością! To równie błędna odpowiedź!

Podejście z użyciem CAS:

Cx[h_,h0_]:=0.1 Exp[-h/h0]

rownanie3 = newton2 /. {F->-m g -k z'[t]-Cx[z[t],h0] z'[t]Abs[z'[t]], a->z''[t]}/. {g->9.81, m->0.5, k->0.1, h0->100}

rozw3 =NDSolve[{rownanie3,z[0]==0,z'[0]==10},z[t],{t,0.0,1}]

z[t] /. rozw3[[1]] /. t->1