Równania różniczkowe

Przygotował: Maciej Maliborski
Uzupełnił: Jacek Golak

Wstęp

Bardzo często, chcąc zrozumieć zachowanie układu fizycznego, napotykamy równania, lub układy równań w których występują nie konkretne wartości liczbowe, ale FUNKCJE. Jeśli warunki na te funkcje zadane są z użyciem pochodnych, to mamy do czynienia z równaniami różniczkowymi. 
Jeśli w tych równaniach występują funkcje jednej zmiennej, to są to równania różniczkowe zwyczajne. Jeśli w równaniach występują funkcje wielu zmiennych, to są to równania różniczkowe cząstkowe.

Nie ma fizyki bez równań różniczkowych !


gdyż w ogólnym przypadku siła może zależeć w sposób jawny od czasu, ale także od prędkości i położenia. (Zakładamy, że nie zależy od przyspieszenia, ani wyższych pochodnych wektora położenia.)

Równania różniczkowe napotykamy w zasadzie w każdej dziedzinie fizyki: w mechanice punktów materialnych i brył sztywnych, w dynamice ośrodków ciągłych, w elektromagnetyźmie (słynne równania Maxwella), w nauce o cieple (równanie przewodnictwa cieplnego),  w nauce o zjawiskach falowych, w mechanice kwantowej, w OTW, ...

Główne narzędzia: DSolve[], NDSolve[]

W NOF_2.nb pokazaliśmy, że  A*Cos[Sqrt[k/m]*t]+B*Sin[Sqrt[k/m]*t] jest rozwiązaniem równania m*x''[t]+k*x[t]=0, gdzie A i B są zupełnie dowolnymi stałymi. To wymagało jedynie podstawienia funkcji do równania z pochodnymi:

x[t_]:=A*Cos[Sqrt[k/m]*t]+B*Sin[Sqrt[k/m]*t]

Simplify[m*x''[t]+k*x[t]==0]

ClearAll[x,A,B,k,m];

Teraz stawiamy sobie cel bardziej ambitny, bo chcemy znaleźć rozwiązanie dla zadanego równania różniczkowego

Zaczynamy od prostego równania: szukamy funkcji, która jest równa swojej pierwszej pochodnej. Równanie różniczkowe jest pierwszego rzędu, bo mamy w nim tylko pierwszą pochodną.

sol=DSolve[x'[t]==x[t],x[t],t]

Rozwiązanie (tzw. całkę ogólną równania różniczkowego) dostajemy w postaci reguły, jako listę, ponieważ w ogólności może istnieć więcej niż jedno rozwiązanie.  (C[1]) oznacza stałą całkowania. Symbol (C[i]) posiada atrybut “Protected”.

(* Tak ...  *)

(C[1])=2

(* ... ani tak nie można nadać stałej całkowania wartości *)

C[1]=2

(* Można jednak za stałą coś podstawić i wtedy przypisanie wartości jest dozwolone *)

x[t]/.First[sol]/.C[1]->K

K=2;(E^t) K

Kolejny przykład pokazuje, jak wydobyć rozwiązanie z DSolve

(* Można w dwóch krokach, analogicznie jak to robiliśmy dla Solve ... *)

s2=DSolve[{x'[t]==2x[t]+Cos[t]},x[t],t]

x[t]/.s2[[1]]

(* ... ale można po prostu tak: *)

x[t]/.First@DSolve[{x'[t]==2x[t]+Cos[t]},x[t],t]

Dodanie warunków początkowych prowadzi do jednoznacznego rozwiązania

(* Dla pierwszego prostego równania I rzędu wystarczy wartość funkcji x[t] dla t=0 *)

x[t]/.First@DSolve[{x'[t]==x[t],x[0]==2},x[t],t]

(* Dla równania oscylatora harmonicznego, które jest równaniem II rzędu, trzeba podać x[0] oraz x'[0], czyli wartość funkcji i pochodnej dla t=0 *)

Simplify[x[t]/.First@DSolve[{x''[t]==-(ω^2)x[t],x[0]==x0,x'[0]==v0},x[t],t],ω>0]

ClearAll[x];

W wielu przypadkach zakładamy, że równanie różniczkowe i warunki początkowe (w danym przypadku położenie początkowe i prędkość początkowa) wyznaczają jednoznaczne rozwiązanie, ale trzeba mieć świadomość, że twierdzenia matematyczne, które to zagadnienie opisują są bardzo skomplikowane.

Kolejny znany i łatwy do rozwiązania przykład: ruch punktu materialnego w jednorodnym ziemskim polu grawitacyjnym. Ruch wystarczy rozpatrywać w płaszczyźnie xy, przy czym zakładamy, że przyspieszenie ziemskie jest skierowane przeciwnie do osi y. 
II zasada dynamiki daje tylko              
                                                                       m (OverVector[a]) = m (OverVector[g]), czyli  (OverVector[a]) = (OverVector[g]) .
Dopiero warunki początkowe prowadzą do jednoznacznego rozwiązania (rzut pionowy w górę, rzut pionowy w dół, spadek swobodny z zadanej wysokości, rzut poziomy, rzut ukośny). Poniżej pokazuję rozwiązanie dla najbardziej ogólnego przypadku, rozważając rzut ukośny z wysokości H.

r[t_]={x[t],y[t]}; (* wektor położenia *)

v[t_]=r'[t] ;(* wektor prędkości *)

a[t_]=r''[t];  (* wektor przyspieszenia *)

r0={0,H}; (* wektor położenia w chwili t=0 *)

v0={(Subscript[v, 0])*Cos[alfa],(Subscript[v, 0])*Sin[alfa]};  (* prędkość w chwili t=0 *)

vecg={0,-g}; (* wektor przyspieszenia ziemskiego, gdzie g > 0 *)

rzut=DSolve[{m*a[t]==m*vecg, r[0]==r0,r'[0]==v0},{x[t],y[t]},t]

Zamiast warunków początkowych do określenia stałych całkowania można (próbować) użyć warunków brzegowych, to znaczy zadać wartości funkcji dla dwóch różnych argumentów.  
Uwaga: Nie zawsze rozwiązanie istnieje !
Przykład

DSolve[{y''[x]==x,y[0]==1,y[1]==-1},y[x],x]

(* Tu się udało *)

DSolve[{y''[x]==-y[x],y[0]==0,y[2*Pi]==1},y[x],x]

(* W tym przypadku się nie udało. Dlaczego ? *)

Równania różniczkowe często mają rozwiązania w postaci funkcji specjalnych. W wielu przypadkach ich rozwiązania wręcz definiują funkcje specjalne.

Przykład

DSolve[ (x^2)y''[x]+x y'[x]+(x-n)(x+n)y[x]==0,y[x],x]

(* W tym przypadku ogólne rozwiązanie jest dowolną kombinacją liniową funkcji Bessela pierwszego (BesselJ[n,x]⟷ (Subscript[J, n])(z)) i drugiego rodzaju (BesselY[n,x]⟷ (Subscript[Y, n])(z)). *)

Plot[{BesselJ[2,x],BesselY[2,x]},{x,0,2},PlotLegends->"Expressions"]

Inny przykład

DSolve[y''[x]+Cos[x] y[x]==0,y,x]

Jeszcze inny przykład

DSolve[y''[x]-x y[x]==0,y[x],x]

Jak pokazano dla rzutu ukośnego, przy pomocy funkcji DSolve możemy rozwiązywać układy równań różniczkowych zwyczajnych. Proszę przyjrzeć się składni tej komendy używającej “@” !

{x[t],y[t]}/.First@
DSolve[{x'[t]==2 y[t]-x[t],y'[t]==-x[t]+y[t]},{x[t],y[t]},t]

(* Dla mnie prościej wygląda taki zapis (polecam): *)

s13=DSolve[{x'[t]==2 y[t]-x[t],y'[t]==-x[t]+y[t]},{x[t],y[t]},t]

{x[t],y[t]}/.s13[[1]]

(* Gdybym chciał mieć do dyspozycji te rozwiązania w postaci funkcji, wtedy: *)

x[t_]=x[t]/.s13[[1]]

y[t_]=y[t]/.s13[[1]]

ClearAll["Global`*"];

Równania w układzie równań mogą mieć charakter różniczkowo-algebraiczny, co oznacza, że nie wszystkie równanie muszą być różniczkowe. (W drugim równaniu nie ma pochodnej !)

DSolve[{y'[x]+3z'[x]==4 y[x]+1/x,y[x]+z[x]==1},{y[x],z[x]},x]

(* Tu także pojawiają się funkcje specjalne *)

Poprawność rozwiązania analitycznego możemy sprawdzić poprzez różniczkowanie wyniku. Jednak nie zawsze równanie daje się scałkować. Rozwiązania w postaci algebraicznej (inaczej analityczne) istnieją tylko dla wybranej klasy równań. Co więcej, rozwiązanie nietrywialnego równania może być bardzo skomplikowane, a jego analiza może przyprawiać o zawrót głowy. Dlatego czasem trzeba się zadowolić rozwiązaniem w postaci numerycznej.

Wynik procedur numerycznych podawany jest jako funkcja interpolacyjna (interpolująca ?). Zobaczmy najpierw, co to oznacza dla znanego przypadku. Funkcja interpFnc ma imitować sinus i jest zbudowana z wartości funkcji sinus w punktach x= 0, Pi/2, Pi, 3*Pi/2, 2*Pi, 5*Pi/2, 3*Pi, 7*Pi/2,4*Pi.

Table[{x,Sin[x]},{x,0,4 Pi,Pi/2}]

interpFnc=Interpolation[ Table[{x,Sin[x]},{x,0,4 Pi,Pi/2}]];

Plot[ {interpFnc[t],Sin[t]},{t,0,4Pi} ,PlotLegends->"Expressions"]

Interpolującą funkcją można (i należy) posługiwać się w ten sam sposób jak zwykłymi funkcjami.

Plot[ {Cos[t],interpFnc'[t]},{t,0,4 Pi}]

dwiecałki=Integrate[{Sin[t],interpFnc[t]},{t,0,Pi}]

N[dwiecałki]

N[dwiecałki[[2]]/dwiecałki[[1]]]

(* Przy tej całce, która jest polem powierzchni pod wykresem Sin[t] i interpFnc[t] dla 0<= x <= Pi, mamy około 5 % różnicy  *)

Funkcja interpolująca nie daje absolutnie dokładnych wyników. Aby rachunki z taką funkcją były precyzyjne, trzeba znać interpolowaną funkcję w odpowiednio wielu "punktach podparcia".

W jaki sposób uzyskać dostęp do punktów, w których została spróbkowana funkcja podczas interpolacji ?

interpFnc=Interpolation[ Table[{x,Sin[x]},{x,0,4 Pi,Pi/2}] ]

(* Tak dostaniemy pełną informację o interpFnc *)

interpFnc//N//FullForm

(* Tak wydobywamy punkty x użyte do zbudowania funkcji, ... *)

xi=interpFnc[[3,1]]

(* ... a tak wartości funkcji w tych punktach *)

fi=N[interpFnc][[4,3]]

Sin[xi]

Sin[xi]==fi

Opcje funkcji Interpolation[].
Metoda “Spline” generuje funkcję o ciągłej pochodnej, w przeciwieństwie do metody “Hermite”

Options[Interpolation]

ClearAll["Global`*"];

Do rozwiązywania numerycznego równań różniczkowych służy NDSolve. Składnia funkcji NDSolve[] różni się od składni DSolve tym, że wymaga podania warunków początkowych (lub brzegowych) oraz zakresu zmiennej niezależnej.
Wszystkie wartości numeryczne parametrów równania muszą być określone.

Przykład 1

(* Tu mamy proste równanie I rzędu z jednym warunkiem początkowym. Dla tego równania mamy rozwiązanie dokładne i możemy sprawdzić rozwiązanie numeryczne *)

ana1=DSolve[{y'[x]==3*y[x],y[0]==1},y[x],x]

num1=NDSolve[{y'[x]==3*y[x],y[0]==1},y[x],{x,0,2}]

Plot[{(y[x]/.ana1[[1]]),(y[x]/.num1[[1]])},{x,0,2},PlotLegends->"Expressions"]

(* Poniżej rysunek względnej różnicy między rozwiązaniem analitycznym (dokładnym) i numerycznym (przybliżonym)*)

Plot[{((y[x]/.ana1[[1]])-(y[x]/.num1[[1]]))/(y[x]/.ana1[[1]])},{x,0,2},PlotLegends->"Expressions"]

Przykład 2 (S. Wolfram)

(* Tu mamy równanie II rzędu z dwoma warunkami brzegowymi *)

num2=NDSolve[{y''[x]+x*y[x]==0,y[0]==1,y[1]==-1},y[x],{x,0,1}]

Plot[y[x]/.num2[[1]],{x,0,1}]

Dalsze przykłady

(* Teraz pokazuję po kolei składniki problemu *)

(* Samo równanie: *)

eq=x'[t]+t ((Cos[Pi t])^2) x[t]==t Cos[t];

(* Warunki początkowe *)

wp=x[0]==0;

(* Przedział, w którym szukamy rozwiązania *)

{tmin,tmax}={0,10};

sol=x/.First@NDSolve[{eq,wp},x,{t,tmin,tmax}]

Plot[Evaluate[sol[t]],{t,0,10}]

Używamy InterpolatingFunction[] jako zwyczajnej funkcji. Tym razem użyjemy jej jako argumentu funkcji FindRoot[].

FindRoot[sol[t]==-3,{t,10}]

Możemy sprawdzić poprawność rozwiązania poprzez wstawienie wyniku funkcji NDSolve[] do równania

Plot[eq/.x->sol,{t,tmin,tmax},Evaluated->True,PlotRange->All]

Domyślną wartością opcji AccuracyGoal (PrecisionGoal) funkcji NDSolve jest (dla precyzji maszynowej) 7 lub 8. Dlatego można się było spodziewać, że uzyskany wynik będzie rozwiązaniem z dokładnością (10^(-7)) czy (10^(-8)). Powyższy wykres wskazuje jednak coś zupełnie innego, bo widzimy znacznie większe wartości.

Okazuje się, że jeśli sprawdzić rozwiązanie w punktach, dla których algorytmy zaimplementowane w NDSolve wyliczyły rozwiązanie, wynik będzie zgodny z oczekiwaniem.

First@sol["Coordinates"]//Short (* punkty siatki w których jest znane rozwiązanie oraz pochodne *)

Table[{t,x'[t]+t ((Cos[Pi t])^2) x[t]-t Cos[t]/.x->sol},{t,First[sol["Coordinates"]]}];
ListPlot[%,PlotRange->All,Joined->True,Mesh->All]

(* Teraz widzimy znacznie mniejsze wartości, tak jak tego oczekiwaliśmy *)

ClearAll["Global`*"];

Pułapki obliczeń numerycznych

Typowe błędy w rachunkach numerycznych spowodowane są przez zaokrąglenia, które są efektem ubocznym skończonej precyzji obliczeń. Co więcej, rachunki związane z równaniami różniczkowymi posiadają dodatkową pułapkę, mianowicie małe błędy mogą propagować się w rozwiązaniu. Jeśli problem jest bardzo wrażliwy na tego typu niewielkie zaburzenia, otrzymane wyniki mogą być dalekie od prawdziwych. 
Obrazują to poniższe przykłady.

Przykład 1

ode={x'[t]==2 x[t]+Cos[t],x[0]==-2/5};

xs=x/.First@NDSolve[ode,x,{t,0,3 π}];

(* Wykres dla rozwiązania numerycznego *)

Plot[xs[t],{t,0,3 π},PlotRange->All]

Dla tego problemu, możemy znaleźć rozwiązanie analityczne

X=x/.First@@DSolve[ode,x,t];

(* Możemy więc porównać wynik numeryczny (niebieski) z analitycznym (pomarańczowym) *)

Plot[{xs[t],X[t]},{t,0,3 π},PlotRange->All]

Olbrzymia rozbieżność wyników jest związana z wartością początkową -(2/5) (do przeanalizowania).

To niezwykle szczęśliwy przypadek dla którego znamy rozwiązanie analityczne. W ogólnym przypadku, kiedy takiego rozwiązania nie mamy, należy stosować i polegać na innych metodach, aby zdobyć pewność poprawności przybliżonego numerycznego rozwiązania.

Przykład 2: równanie Duffinga

Czteroparametrowe równanie pojawiające się przy analizie ruchu oscylatora wymuszonego
 
x''(t)-a x(t)+b ((x(t))^3)+c x'(t)=d cos(t) .

Dla pewnych wartości początkowych zachowanie rozwiązania jest bardzo skomplikowane.

Clear[ode,x,xs]

odeDuff[a_,b_,c_,d_]:=x''[t]-a x[t]+b ((x[t])^3)+c x'[t]==d Cos[t]

ode={odeDuff[1,1,15/100,3/10],x[0]==1,x'[0]==-1};
xs=x/.First@NDSolve[ode,x,{t,0,100}];

(* Teraz wprowadzamy bardzo niewielką zmianę (eps) w warunkach początkowych i sprawdzamy, jak zmieni się rozwiązanie *)

eps=(10^(-8));
ode={odeDuff[1,1,15/100,3/10],x[0]==1+eps,x'[0]==-1+eps};
xseps=x/.First@NDSolve[ode,x,{t,0,100}];

(* Porównujemy wyniki dla przypadków, które bardzo nieznacznie różnią się warunkami początkowymi *)

Plot[ {xs[t],xseps[t]},{t,0,100} ]

Począwszy od około t=75 zaburzone rozwiązanie różni się znacząco od oryginalnego rozwiązania. To mocny dowód na to, że oba rozwiązania są nieprawidłowe dla dużych t. Co więcej, zgodność obu rozwiązań do tego punku może być potraktowana jako dowód na to, że oba rozwiązania są poprawne.

(* Mamy sposób, poprzez WorkingPrecision, by uczynić rozwiązanie bardziej dokładnym *)

ode={odeDuff[1,1,15/100,3/10],x[0]==1,x'[0]==-1};
xs=x/.First@NDSolve[ode,x,{t,0,100},WorkingPrecision->20];
eps=(10^(-8));
ode={odeDuff[1,1,15/100,3/10],x[0]==1+eps,x'[0]==-1+eps};
xseps=x/.First@NDSolve[ode,x,{t,0,100},WorkingPrecision->20];

Plot[ {xs[t],xseps[t]},{t,0,100} ]

Zwiększając precyzję obliczeń (poprzez ustawienie parametru WorkingPrecision->20) możemy wydłużyć czas, po jakim oba rozwiązania się rozbiegają. Jest to pesymistyczny przypadek, dla którego całkowanie dla dużych czasów wymaga coraz większej precyzji.

Przykład 3

ode={x''[t]==((x'[t])^2)-x[t],x[0]==14,x'[0]==0};

xs=x/.First@NDSolve[ode,x,{t,0,16.5}];

Plot[xs[t],{t,0,16.5},PlotRange->All]

Częstym sposobem sprawdzenia poprawności otrzymanych wyników numerycznych jest całkowanie równania wstecz.

ode={x''[t]==((x'[t])^2)-x[t],x[16.5]==xs[16.5],x'[16.5]==xs'[16.5]};
xb=x/.First@NDSolve[ode,x,{t,16.5,0}];

Plot[{xs[t],xb[t]},{t,0,16.5},PlotRange->All]

{xs[0],xb[0]}

(* Widać, że te wyniki nie zgadzają się idealnie i dają informację o dokładności obliczeń. *)

ClearAll["Global`*"];

Informacje o problemach związanych z numerycznym rozwiązywaniem równań różniczkowych nie powinny Państwa zniechęcać. Jako fizycy mamy własne sposoby, by sprawdzić jakość rozwiązań - zasady zachowania !
Pokazuję to dla nieco ambitniejszego problemu, w którym szukamy rozwiązania układu równań różniczkowych II rzędu z czterema funkcjami. Chodzi o tzw. problem Keplera.

Dwie masy, (Subscript[m, 1]) i (Subscript[m, 2]) , poruszają się w ustalonej płaszczyźnie pod wpływem wzajemnej siły grawitacji. Opisujemy ruch tego układu, rozwiązując numerycznie równanie Newtona dla wektorów położenia (Subscript[(OverVector[r]), 1]) i  (Subscript[(OverVector[r]), 2]) .

 (Subscript[m, 1])= 4 MS,  (Subscript[m, 2])= MS

r1={x1[t],y1[t]};

r2={x2[t],y2[t]};

r=r2-r1;

(* II zasada dynamiki dla masy m1 *)

eq1=m1*D[r1,{t,2}]==G*m1*m2/(Sqrt[r.r])^3*r

(* II zasada dynamiki dla masy m2 *)

eq2=m2*D[r2,{t,2}]==-G*m1*m2/(Sqrt[r.r])^3*r

(* MS=2*10^(30) *) (* kg , masa naszego Słońca *)

MS=1;

m1=4*MS;

m2=MS;

AU=1; (* jednostka astronomiczna, czyli 150 milionów kilometrów jest teraz jednostką długości ! *)

rmin=24*AU;

rmax=36*AU;

a=(rmin+rmax)/2;

(* ziemski rok jest teraz jednostką czasu *)

rok=1;

T=80*rok;

(* Po wyborze okresu T oraz parametru elipsy w ruchu względnym, wartość stałej G jest już określona i możemy ją znaleźć z III prawa Keplera *)

solG=Solve[2*Pi*a^(3/2)/Sqrt[G*(m1+m2)]==T,G];

G=G/.solG[[1]];

tmax=1500;

(* Warunki początkowe w odpowiednich jednostkach. Są one tak dobrane, by środek masy układu pozostawał w spoczynku *)

WP={x1[0]==-912/100,y1[0]==0,x2[0]==3648/100,y2[0]==0,x1'[0]==0,y1'[0]==-48/100,x2'[0]==0,y2'[0]==192/100};

(* Używam funkcji NDSolve, aby numerycznie rozwiązać układ równań różniczkowych *)

numerycznie=NDSolve[{eq1,eq2,WP},{x1[t],y1[t],x2[t],y2[t]},{t,0,tmax}];

(* Wydobywam rozwiązania w postaci funkcji *)

x1[t_]:=Evaluate[x1[t]/.Flatten[numerycznie]]

y1[t_]:=Evaluate[y1[t]/.Flatten[numerycznie]]

x2[t_]:=Evaluate[x2[t]/.Flatten[numerycznie]]

y2[t_]:=Evaluate[y2[t]/.Flatten[numerycznie]]

ClearAll[r1,r2,r];

(* Na podstawie numerycznych rozwiązań tworzę wektory położenia obu mas oraz wektor położenia względnego *)

r1[t_]:={x1[t],y1[t]};

r2[t_]:={x2[t],y2[t]};

r[t_]:=r2[t]-r1[t];

(* całkowita (dla obu mas) energia kinetyczna *)

Ekin[t_]=(1/2)*m1*D[r1[t],t].D[r1[t],t]+(1/2)*m2*D[r2[t],t].D[r2[t],t];

(* całkowita (dla obu mas) energia potencjalna *)

Epot[t_]=-G*m1*m2/Sqrt[r[t].r[t]];

(* Całkowita energia mechaniczna układu dwóch mas *)

Etot[t_]=Ekin[t]+Epot[t];

(*
px1=Plot[x1[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [rok]","(Subscript[x, 1])(t) [AU]"}];
py1=Plot[y1[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [rok]","(Subscript[y, 1])(t) [AU]"}];
px2=Plot[x2[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [rok]","(Subscript[x, 2])(t) [AU]"}];
py2=Plot[y2[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [rok]","(Subscript[y, 2])(t) [AU]"}];
*)

(* Tory obu mas, niebieski dla większej masy m1 *)

p1=ParametricPlot[{{x1[t],y1[t]},{x2[t],y2[t]}},{t,0,tmax},AxesLabel->{"(Subscript[x, i])(t) [AU]","(Subscript[y, i])(t) [AU]"}]

(* Sprawdzam jakość rozwiązania: całkowita energia mechaniczna nie może się zmieniać w czasie *)

Plot[Etot[t],{t,0,tmax},PlotRange->All,AxesLabel->{"t [rok]","Etot(t) "}]

(* Numeryczne zmiany całkowitej energii są bardzo niewielkie, ale warto je dodatkowo porównać z wartościami energii kinetycznej, potencjalnej i całkowitej *)

Plot[{Ekin[t],Epot[t],Etot[t]},{t,0,tmax},PlotRange->All,AxesLabel->{"t [rok]","Ekin(t), Epot(t), Etot(t) "},PlotLegends->"Expressions"]

(* Na tym rysunku widać, że spełniona jest zasada zachowania całkowitej energii mechanicznej *)

Modelowanie ruchu układu przy pomocy wbudowanych obiektów graficznych

masa1[t_]=Disk[{x1[t],y1[t]},m1];

masa2[t_]=Disk[{x2[t],y2[t]},(1/2)*m1];

gmasa1[t_]:=Graphics[masa1[t]]

gmasa2[t_]:=Graphics[masa2[t]]

animacja=Animate[Show[p1,gmasa1[t],gmasa2[t],PlotRange->{{-140,50},{-70,70}},AspectRatio->Automatic],{t,0,tmax}]

(* SetDirectory[NotebookDirectory[]]; Export["2masy_2.avi",animacja]; *)

ClearAll["Global`*"];

Równania różniczkowe cząstkowe

DSolve służy nie tylko do rozwiązywania równań i układów równań różniczkowych zwyczajnych dla funkcji jednej zmiennej, ale także równań różniczkowych cząstkowych, gdzie występują funkcje wielu zmiennych.

W tym przypadku uzyskanie jednoznacznego rozwiązania wymaga zwykle podania dodatkowych warunków (tzw. warunków początkowych i brzegowych). Rozwiązanie ogólne może zależeć nie od dowolnych parametrów, ale wręcz od dowolnych funkcji !

Przykład 1 : równanie I rzędu na funkcję dwóch zmiennych 
(S. Wolfram)

DSolve[D[y[x1,x2],x1]+D[y[x1,x2],x2]==1/(x1 x2),y[x1,x2],{x1,x2}]

(* Rozwiązanie ogólne zależy od dowolnej funkcji różnicy x2-x1. Mathematica nazywa ją (C[1]). *)

Przykład 2 : problem dotyczący pewnego rodzaju sił, dla których ma istnieć energia potencjalna.
Mamy sprawdzić, czy siły F={2*x*z*z-2*y,-2*x-6*y*z,2*x*x*z-3*y*y} i G={x*x*z,-x*y,5} należą do tej kategorii i jeśli tak, znaleźć odpowiadające im energie potencjalne.

F:={2*x*z*z-2*y,-2*x-6*y*z,2*x*x*z-3*y*y};

rotF=Curl[F,{x,y,z}]

G:={x*x*z,-x*y,5};

rotG=Curl[G,{x,y,z}]

(* Gradient funkcji skalarnej to wektor, którego składowymi są pochodne odpowiednio po x, y i z. W ogólnym zapisie mamy: *)

Grad[f[x,y,z],{x,y,z}]

(* Powyższy zapis może nie być zbyt czytelny, dlatego sprawdzam bezpośrednio, czym jest gradient *)

Grad[f[x,y,z],{x,y,z}]=={D[f[x,y,z],x],D[f[x,y,z],y],D[f[x,y,z],z]}

(* Szukamy energii potencjalnej dla siły F *)

DSolve[Grad[VF[x,y,z],{x,y,z}]==-F,VF[x,y,z],{x,y,z}]

(* Rozwiązanie na energię potencjalną jest określone z dokładnością do stałej *)

DSolve[Grad[VG[x,y,z],{x,y,z}]==G,VG[x,y,z],{x,y,z}]

(* Zgodnie z oczekiwaniem, dla siły G Mathematica nie znalazła rozwiązania *)

Przykład 3: jednowymiarowe (jedna współrzędna przestrzenna x) równanie falowe.
Mamy do czynienia ze struną rozciągającą się wzdłuż osi x, której wychylenie z położenia równowagi zależy od położenia wzdłuż osi x i od czasu.

ClearAll["Global`*"];

$Assumptions=c>0;

eq3={D[u[x,t],{x,2}]-(1/c^2)*D[u[x,t],{t,2}]==0}

s3=Simplify[DSolve[eq3,u[x,t],{x,t}]]

W równaniu powyższym nie występowały ani warunki początkowe, ani warunki brzegowe. Dlatego rozwiązaniem jest suma dwóch dowolnych funkcji jednej zmiennej, (C[1]) i (C[2]), z argumentami równymi odpowiednio t-x/c oraz t+x/c. Te dwie części rozwiązania odpowiadają fali biegnącej zgodnie ze zwrotem osi x (argument t-x/c) oraz przeciwnie do zwrotu osi x (t+x/c). 

Teraz wprowadzamy warunki początkowe i brzegowe.
Warunki początkowe dla równania drugiego rzędu w czasie oznaczają, że musimy podać postać funkcji u i jej pierwszej pochodnej w chwili t=0.

Dla fali na strunie o skończonej długości, x ∈ [-L,L], odpowiada to podaniu kształtu struny w chwili t=0 oraz rozkładu prędkości punktów struny w chwili t=0.

wp3={u[x,0]==Exp[-(x^2)],Derivative[0,1][u][x,0]==0}

Warunki brzegowe: zakładamy, że fala występuje na strunie, która jest sztywno umocowana na obu końcach. Konkretne wartości wybieram tak, by warunki początkowe były zgodne z warunkami brzegowymi.

L=1;

wb3={u[-L,t]==E^(-1),u[L,t]==E^(-1)}

c=1;tmax=5;

s31=NDSolve[{eq3,wp3,wb3},u,{x,-L,L},{t,0,tmax}]

Plot3D[Evaluate[u[x,t]/.s31[[1]]],{t,0,tmax},{x,-L,L},AxesLabel->{"t","x","u(x,t)"}]

Animate[Plot[Evaluate[u[x,t]/.s31[[1]]],{x,-L,L},PlotRange->{-1/2,1.1},AxesLabel->{"x","u(x,t)"}],{t,0,4}]

Dalsze przykłady są dostępne na stronie przedmiotu
http://users.uj.edu.pl/~golak/zestawyNOF.html
oraz na stronie
http://users.uj.edu.pl/~golak/zestawyF.html