Równania cząstkowe w Mathematice

Definiowanie równań cząstkowych i warunków brzegowych

Funkcje, o których mowa była w pierwszej części zajęć (dotyczącej równań zwyczajnych) mogą posłużyć do rozwiązywania równań cząstkowych. Mathematica oferuje kilka możliwości zapisania równania:

pde1=D[u[t,x],t]+2D[u[t,x],x]==0

pde1//InputForm

pde1//OutputForm

pde1//StandardForm

pde1//TraditionalForm

Powyższe przypisanie definiuje nam równanie transportu i przy okazji pokazuje inny sposób wyświetlania pochodnych (cząstkowych). Można też jeszcze inaczej:

pde2=y (∂_x)u[x,y]-(∂_y)u[x,y]==0

gdzie do wprowadzenia pochodnych używamy skrótu Esc pd Esc. Co ciekawe, możemy też od razu wygenerować kod w TEX-u:

pde2//TeXForm

Jeśli Mathematica potrafi rozwiązać równanie cząstkowe w sposób analityczny, to używa do tego znanej funkcji:

sol1=DSolve[pde1,u[t,x],{t,x}]

Możemy tej funkcji użyć zarówno do równań liniowych i quasi-liniowych (dla których jest szczególnie użyteczna) przy szukaniu rozwiązań ogólnych, jak i szczególnych, ale także do równań nieliniowych (choć w ograniczonym zakresie).

Widzimy, że podobnie jak przy równaniach zwyczajnych pojawiają się parametry C[1] (i w przypadku równań wyższych rzędów może ich być więcej), ale tym razem są to funkcje (jednej zmiennej).

sol1A=DSolve[pde1,u[t,x],{t,x},GeneratedParameters->((Subscript[f, #])&)]

Żeby wybrać konkretną funkcję C[1], można zadać warunek brzegowy:

sol2=DSolve[{pde1,u[t,0]==Sin[t]},u[t,x],{t,x}]

i dostajemy wówczas konkretną funkcję, którą możemy wyświetlić:

Plot3D[u[t,x]/.sol2,{t,-4,4},{x,-4,4},AxesLabel->Automatic]

Warunek brzegowy nie musi być zdefiniowany na prostej równoległej do osi układu:

sol3=DSolve[{pde2,u[x,x]==Sin[x]},u[x,y],{x,y}]

Plot3D[u[x,y]/.sol3,{x,-4,4},{y,-4,4},AxesLabel->Automatic]

Zauważmy, że Mathematica rozwiązuje równanie symbolicznie, zupełnie nie zastanawiając się nad tym, gdzie rozwiązanie równania ma w ogóle sens.

Jak sprawdzić, czy znaleziona funkcja jest rozwiązaniem podanego równania?

sol3B=DSolve[{pde2,u[x,x]==Sin[x]},u,{x,y}]

pde2/.sol3B

Mathematica radzi sobie z większością równań pierwszego rzędu (choć nie ze wszystkimi), nie tylko liniowymi, np. z tymi, które były na zajęciach:

DSolve[{u[x,y] D[u[x,y],x]+D[u[x,y],y]==1,u[x,x]==(x/2)},u[x,y],{x,y}]

Możemy też rozwiązywać równania, których rozwiązanie dane jest w postaci uwikłanej:

pde3=(x-u[x,y])D[u[x,y],x]+(y-u[x,y])D[u[x,y],y]==0;

sol4=DSolve[pde3,u[x,y],{x,y}]

Jak sprawdzić, że podana funkcja jest rzeczywiście rozwiązaniem równania?

D[sol4[[1]],x]

px=Solve[D[sol4[[1]],x],D[u[x,y],x]][[1]]
py=Solve[D[sol4[[1]],y],D[u[x,y],y]][[1]]

pde3/.px/.py//Simplify

Równania drugiego rzędu

Jednym z pierwszych zadań na Równaniach cząstkowych (oraz na Analizie zespolonej) było znalezienie operatora Laplace’a w innych układach współrzędnych. Mathematica obsługuje nie tylko układ kartezjański:

(* Needs["VectorAnalysis`"] - być może potrzebny *)

Laplacian[u[x,y,z],{x,y,z},"Cartesian"]==0

Laplacian[u[x,y,z],{x,y,z},"Cylindrical"]==0 // Expand

Laplacian[u[x,y,z],{x,y,z},"Spherical"]==0 // Expand

Radzi sobie też nieźle z rozwiązywaniem równań drugiego rzędu (liniowych, o stałych współczynnikach...):

DSolve[3D[u[x,y],{x,2}]-12D[u[x,y],{y,2}]==1,u[x,y],{x,y}]

Nie tylko równania o konkretnych współczynnikach

Możemy też wprowadzać równania, które mają współczynniki będące symbolami.

sol5=DSolve[{a D[u[t,x],t]+b D[u[t,x],x]==0,u[0,x]==Sin[(x^2)]},u[t,x],{t,x}]

Można wówczas narysować wiele różnych funkcji:

F[t_,x_]:=u[t,x]/.sol5

Manipulate[
Plot3D[F[t,x]/.{a->A,b->B},{t,0,8},{x,-4,4},AxesLabel->Automatic],
{A,-2,2},{B,-2,2}
]

Rozwiązania równań nieliniowych wyższych rzędów

W pewnych szczególnych przypadkach daje się też znaleźć rozwiązania nieliniowych równań wyższych rzędów, np. równania Kortewega-de Vries:

KdV=(∂_t)u[t,x]+6u[t,x] (∂_x)u[t,x]+(∂_({x,3}))u[t,x]==0;

sol6=DSolve[KdV,u,{t,x}]

Widzimy jednak, że znaleziono tylko jedno z możliwych rozwiązań.

Rozwiązania numeryczne

Jeśli Mathematica nie radzi sobie z rozwiązaniem równania analitycznie, to warto spróbować rozwiązania numerycznego. Tutaj znowu działa standardowa funkcja, trzeba jednak uważać na warunki brzegowe - rozwiązanie będzie odbywało się w ograniczonym zbiorze!

sol7=NDSolve[
{D[u[t,x],t]==D[u[t,x],x,x],u[0,x]==0,u[t,0]==Sin[t],u[t,5]==0},
u[t,x],{t,0,10},{x,0,5}]

Plot3D[u[t,x]/.sol7,{t,0,10},{x,0,5},PlotRange->All,AxesLabel->Automatic]

Można też przedstawić wynik jako “animację”:

Manipulate[
Plot[u[t,x]/.sol7/.t->T,{x,0,5},PlotRange->{{0,5},{-1,1}},AxesLabel->{x,u[T,x]}],
{T,0,10}]

Mimo, że rozwiązania są numeryczne, to można też wprowadzić pewne zmienne symbole:

sol[k_]:=
NDSolve[{k (∂_t)u[t,x]==(x^4)(∂_(x,x))u[t,x]-(((x^3)t Cos[x Sin[t]])/((1/10)+Sin[t])),u[t,0]==0,u[t,1]==0,u[0,x]==x(1-x)},u[t,x],{x,0,1},{t,0,1}];

sol[3]

Manipulate[
Plot3D[Evaluate[u[t,x]/.sol[k]],{t,0,1},{x,0,1},PlotRange->{{0,1},{0,1},{-0.25,0.25}}],
{k,(1/100),5}]

Czemu to mimo wszystko działa? Dlaczego nie działa bez Evaluate?