Algebra Symboliczna
Wykład VII

Andrzej Odrzywolek
Instytut Fizyki, Zakład Teorii Względnosci i Astrofizyki
03-12-2008, środa, 13:15

Nowa strona wykładu: http://ribes.if.uj.edu.pl/alsymb/

Dane kontaktowe

dr Andrzej Odrzywołek
pokój 447, IV pietro
E-mail: odrzywolek@th.if.uj.edu.pl
Wykład: środy 13.15-15.00, sala 128
Cwiczenia: piątki 10.30-12.00
Konsultacje: środy 11-13, czwartki 12-13
WWW: http://ribes.if.uj.edu.pl/alsymb/

Równania różniczkowe cząstkowe

R. r. cząstkowe to równania na funkcje wielu zmiennych, w odróżnieniu od r. r. zwyczajnych które opisują funkcje jednej zmiennej. Przykład takiego prostego równania:

((∂F(x,y))/(∂x))+x ((∂F(x,y))/(∂y))==x+sin(y)

Wchodząc w zagadnienia równań różniczkowych cząstkowych należy mieć świadomość, iż jest to obszar matematyki znacznie rozleglejszy, skomplikowany i trudniejszy niż r. r. zwyczajne.

Powyższe stwierdzenie dotyczy metod zarówno symbolicznych jak i numerycznych. Zastosowanie algebry symbolicznej ogranicza się tu więc do celów edukacyjnych, wsparcia skomplikowanych rachunków i tworzenia algorytmów.

Uzyskanie rozwiązania r .r. cząstkowego porzez proste wpisanie DSolve czy NDSolve jest na obecnym etapie rozwoju algebry symbolicznej i metod numerycznych wysoce nieprawdopodobne. Zadania z Zestawu 9 są specjalnie spreparowane i służą ilustracji pewnych zagadnień.

Przykład prostego rozwiązania analitycznego

Proste równanie:

r1=D[f[x,y],x]+x D[f[x,y],y]==x +Sin[y]

I jego rozwiązanie ogólne:

DSolve[r1,f[x,y],{x,y}]

Należy koniecznie zwrócić uwagę na fakt, że ,,stała całkowania'' C[1] jest funkcją (w tym konkretnym przypadku jednej zmiennej). Dlatego przykładowe rozwiązania szczególne to:

rozw1=f[x,y] /. DSolve[r1,f[x,y],{x,y}][[1]]

rozw1a=rozw1 /. C[1]->Sin

rozw1b = rozw1 /. C[1] ->(#^2&) (* UWAGA! W 7.0 ( i 6.0.?) ,,czyste funkcje'' muszą być w nawiasie !!! *)

Rowiązanie typowego konkretnego problemu, np. znaleźć rozwiązanie szczególne które na krzywej przyjmuje jakieś wartości wymaga wyznaczenia nieznanej funkcji C[1], czyli rozwiązania równania funkcyjnego ! Tego kroku nie zrobi za nas MATHEMATICA.

Rozważmy rozwiązanie szczególne równania poniżej, które na krzywej y=x^2 przyjmuje wartości f(x,y^2)=x^4 :

r2= D[f[x,y],x]-x D[f[x,y],y]==x y

rozw2=DSolve[r2,f[x,y],{x,y}]

Wstawiamy równanie krzywej na której mamy zadane wartości:

f[x,y]/.rozw2[[1]] /. y->x^2

Żądamy aby na tej krzywej f=x^4:

rfunkcyjne=%==x^4

Mathematica nie rozwiązuje r. funkcyjnych:

Reduce[%,C[1]]

Ale łatwo zrobić to ,,ręcznie'':

rfunkcyjne /. x->Sqrt[2/3t]

Solve[%,C[1][t]]

Sprawdzamy czy funkcja x^2/6 spełnia żądane warunki:

rozwiazanie=f[x,y]/.rozw2[[1]] /. C[1]->(#^2/6&)

rozwiazanie /. y->x^2

Sprawdzamy czy rozwiązanie szczególne spełnia równanie wyjściowe:

g[x_,y_]:=(1/8) ((x^4)+4 (x^2) y+(1/3) ((((x^2)+2 y))^2))

r2 /.f->g

Simplify[%]

Równanie falowe w 1D

Bardzo dobrym przykładem ilustrującym różne zagadnienia związane z równaniami cząstkowymi jest równanie falowe:

Rozwiązanie jawne (wzór d'Alemberta)

rfalowe = D[u[x,t],t,t]==D[u[x,t],x,x]

DSolve znajduje dobrze znane rozwiązanie ogólne:

DSolve[rfalowe,u[x,t],{x,t}]

Rozwiązanie zawiera dwie funkcje dowolne, bo równanie jest drugiego rzędu

rozw=u[x,t] /. %[[1]]

Zmieniamy nazwy f. dowolnych C[1] i C[2] na f i g :

r=rozw/. C[1]->f /. C[2]->g /. t-x->x-t (* Ostatnia zamiana znaku pozwala ominąć późniejsze trudności z DSolve *)

Zadanie kształtu impulsu w chwili t=0 nie wyznacza rozwiązania jednoznacznie. Niech np. u[x,0]=Exp[-x^2].

rimpuls=(r /. t->0)==Exp[-x^2]

Oto trzy przykładowe rozwiazania:

rimpuls /. {f->(0&),g->(Exp[-#^2]&)}

rimpuls /. {g->(0&),f->(Exp[-#^2]&)}

rimpuls /. {f->(1/2 Exp[-#^2]&),g->(1/2 Exp[-#^2]&)}

impuls1=r /. {g->(0&),f->( Exp[-#^2]&)}

impuls2=r /. {f->(0&),g->( Exp[-#^2]&)}

impuls3=r /. {f->(1/2 Exp[-#^2]&),g->(1/2 Exp[-#^2]&)}

Manipulate[Plot[{impuls1,impuls2,impuls3}[[ii]] /. t->T,{x,-10,10},PlotRange->{0,1}],{T,0,10},{ii,{1,2,3}}]

Trzy przykłady powyżej mają różną pochodną cząstkową po t  dla t=0 :

D[impuls1,t] /. t->0

D[impuls2,t] /. t->0

D[impuls3,t] /. t->0

Spróbujemy znaleźć rozwiązanie zagadnienia początkowego (Cauchy'ego) dla r. falowego. Zakładamy, że dla t=0 znamy wartość funkcji i jej pochodnej na całej osi x:

Pierwszy warunek wymaga aby w chwili t=0 rozwiązanie było równe u0[x]:

eq1 = (r /. t->0) == u0[x]

Pierwszy warunek wymaga aby w chwili t=0 pochodna rozwiązanie było równe u1[x]:

eq2=(D[r,t] /. t->0)==u1[x]

Okazuje się, że DSolve znajduje funkcje f i g (jeżeli tu pojawia się komunikat o błędzie należy wcześniej zmienić znak wyrażenia t-x):

DSolve[{eq1,eq2},{f,g},x]

Wyrażenie powyżej znane jest jako wzór d'Alemberta:

dAlembert=r /. % //First //Simplify

Sprawdzamy dla wcześniej zbadanego impulsu:

dAlembert /. {u1->(0&), u0->(Exp[-#^2]&)}

Plot[UnitStep[x+1]*UnitStep[1-x],{x,-2,2}]

Inny przykład, tzw. ,,implozja'' dla -1<x<1 :

przyklad=dAlembert /. {u1->(-UnitStep[#+1]*UnitStep[1-#]&),u0->(0&)}//Simplify

Manipulate[Plot[przyklad /. {x->X,t->T},{X,-4,4},PlotRange->{-2,2}],{T,0,10}]

Manipulate[Plot[D[przyklad,x] /. {x->X,t->T},{X,-4,4},PlotRange->{-2,2}],{T,0,10}]

Manipulate[Plot[D[przyklad,x]^2+D[przyklad,t]^2 /. {x->X,t->T},{X,-4,4},PlotRange->{-2,2}],{T,0,10}]

Rozwiązanie numeryczne

Jeżeli ktoś ma wątpliwości czy powuższe rozwiązania są poprawne może rozwiązać ten sam problem numerycznie:

rozwNUM=u/. NDSolve[{
D[u[x,t],t,t]- D[u[x,t],x,x]==0 (* r. falowe *),
u[x,0]==0 (* f. w chwili t=0 *),
Derivative[0,1][u][x,0]==-UnitStep[1+x]UnitStep[1-x] (* pochodna dla t=0 *),
u[-10,t]==0 (* lewy warunek brzegowy, zakładamy że fala nie dotrze do x=-10 *) ,
u[10,t]==0 (* prawy warunek brzegowy *)},
u (* szukana funkcja *),
{t,0,8} (* dla t<8, fala nie dotrze do x=-10 lub x=10, bo c*t<10 *),
{x,-10,10} (* przedział musi być taki sam jak podany w warunkach brzegowych *),

(* reszta to opcje, nie są wymagane *)Method->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MaxStepSize"->0.1}}][[1]]

Manipulate[Plot[{rozwNUM[X,T],przyklad /. {x->X,t->T}},{X,-4,4},PlotRange->{-2,2}],{T,0,10}]

Rozwiązanie metodą transformaty Fouriera

rfalowe

FourierTransform[D[u[x,t],t,t],x,ω]

FourierTransform[D[u[x,t],x,x],x,ω]

Równanie ,,przetransformowane'':

D[Uu[t],t,t]+ω^2 Uu[t]==0

sol=DSolve[%,Uu[t],t]

C[1] i C[2] należy traktować jak funkcje ω. Dla t=0 dostajemy transformaty warunków początkowych. C[1] to transformata
u0[x], natomiast ω C[2] to transformata u1[x]. Widać to po podstawieniu:

Uu[t] /. sol[[1]] /. t->0

D[Uu[t] /. sol[[1]],t] /. t->0

Rozważmy impuls z wcześniejszych przykładów. Oto jego transformata Fouriera:

FourierTransform[Exp[-x^2],x,ω]

Transformata funkcji u1[x]=0 to oczywiście zero. Aby znaleźć rozwiązanie dokonujemy transformaty odwrotnej:

InverseFourierTransform[((E^(-((ω^2)/4)))/(Sqrt[2])) Cos[t ω],ω,x]

Wynik jest w innej postaci niż wcześniej ale ciągle  ten sam. Odejmowanie Cosh-Sinh powoduje ogromne błędy numeryczne i fałszywie wyglądający wykres. Należy zwiększyć precyzję numeryczną.

Manipulate[Plot[(1/2) (1+(E^(4 t x))) (Cosh[(((t+x))^2)]-Sinh[(((t+x))^2)]),{x,-10,10},PlotRange->{0,1},WorkingPrecision->$MachinePrecision],{t,0,10,1}]

Drugi przykład: ,,implozja'':

FourierTransform[-UnitStep[1+x]UnitStep[1-x],x,ω]

InverseFourierTransform[-(Sqrt[(2/π)]) Sin[ω]Sin[ω t]/ω^2,ω,x]

Manipulate[Plot[(1/4) (-(1+t-x) Sign[1+t-x]-(-1+t-x) Sign[1-t+x]+(-1+t+x) Sign[-1+t+x]-(1+t+x) Sign[1+t+x]),{x,-10,10},PlotRange->{-1,1},WorkingPrecision->$MachinePrecision],{t,0,10,1}]

Metoda rozdzielenia zmiennych

Metoda nadaje się szczególnie do badania rozwiązań w ograniczonych obszarach.

rozdzielenie =D[X[x]T[t],t,t]==D[X[x]T[t],x,x]

DSolve[T''[t]/T[t]==-k^2,T[t],t]

DSolve[{X''[x]/X[x]==-k^2,X[-Pi]==0,X[Pi]==0},X[x],x]

sol = X[x] /. DSolve[{X''[x]/X[x]==-k^2},X[x],x][[1]]

warbrzeg1=sol /. x->-Pi/2

warbrzeg2=sol /. x->Pi/2

Reduce[{warbrzeg1==0&&warbrzeg2==0&&(C[1]≠0||C[2]≠0)},{C[1],C[2],k}] /. C[3]->n

Jak widać poniżej wyrazy z 4n i 2n+1 można przedstawić zbiorczo jako 2 m, gdzie m∈Integers:

Manipulate[Plot[{Sin[2(1+2n )x],Sin[4 n x],Cos[-x+4 n x],Cos[x+4n x]},{x,-Pi/2,Pi/2}],{n,-2,2,1}]

Ogólne rozwiązanie jest sumą nieskończonej ilości wyrazów typu:

(∑_(n=0)^∞)(C1t(m) sin(m t)+C2t(m) cos(m t)) (C1x[m]( sin(m x)+C2x(m) cos(m x-(x/2)))

Ograniczymy się do przypadku gdzie funkcja u0[x] jest nieparzysta (znikają wyrazy z Cos[x]) a funkcja u1[x]=0 (znikają wyrazy z Sin[ n t]). Suma przyjmuje postać:

Sum[C[m]Cos[m t]Sin[2 m x],{n,0,Infinity}]

Manipulate[Plot[{Sin[2 m x]},{x,-Pi/2,Pi/2}],{m,-2,2,1}]

Niech w chwili t=0 mamy u0[x]=1-(2x/Pi)^2.Rozwijamy warunek początkowy w szereg Fouriera:

FourierSinSeries[(1-(2x/Pi)^2)x, x, 6,FourierParameters->{1,2}]

Plot[{(1-(2x/Pi)^2)x,%},{x,-Pi/2,Pi/2}]

Tworzymy sumę:

Sum[FourierSinCoefficient[(1-(2x/Pi)^2)x, x, m,FourierParameters->{1,2}]*Sin[2 m x]*Cos[2 m t],{m,0,16}]

Manipulate[Plot[((6 Cos[2 t] Sin[2 x])/(π^2))-((3 Cos[4 t] Sin[4 x])/(4 (π^2)))+((2 Cos[6 t] Sin[6 x])/(9 (π^2)))-((3 Cos[8 t] Sin[8 x])/(32 (π^2)))+((6 Cos[10 t] Sin[10 x])/(125 (π^2)))-((Cos[12 t] Sin[12 x])/(36 (π^2))),{x,-Pi/2,Pi/2},PlotRange->{-2,2}],{t,0,3}]

rozwNUM2=u/. NDSolve[{
D[u[x,t],t,t]- D[u[x,t],x,x]==0 (* r. falowe *),
u[x,0]==(1-(2x/Pi)^2)x (* f. w chwili t=0 *),
Derivative[0,1][u][x,0]==0 (* pochodna dla t=0 *),
u[-Pi/2,t]==0 (* lewy warunek brzegowy, struna zamocowana na stałe *) ,
u[Pi/2,t]==0 (* prawy warunek brzegowy *)},
u (* szukana funkcja *),
{t,0,100} (* t formalnie dowolne *),
{x,-Pi/2,Pi/2} (* przedział musi być taki sam jak podany w warunkach brzegowych *)][[1]]

Manipulate[Plot[{rozwNUM2[x,t]-(((6 Cos[2 t] Sin[2 x])/(π^2))-((3 Cos[4 t] Sin[4 x])/(4 (π^2)))+((2 Cos[6 t] Sin[6 x])/(9 (π^2)))-((3 Cos[8 t] Sin[8 x])/(32 (π^2)))+((6 Cos[10 t] Sin[10 x])/(125 (π^2)))-((Cos[12 t] Sin[12 x])/(36 (π^2)))+((6 Cos[14 t] Sin[14 x])/(343 (π^2)))-((3 Cos[16 t] Sin[16 x])/(256 (π^2)))+((2 Cos[18 t] Sin[18 x])/(243 (π^2)))-((3 Cos[20 t] Sin[20 x])/(500 (π^2)))+((6 Cos[22 t] Sin[22 x])/(1331 (π^2)))-((Cos[24 t] Sin[24 x])/(288 (π^2)))+((6 Cos[26 t] Sin[26 x])/(2197 (π^2)))-((3 Cos[28 t] Sin[28 x])/(1372 (π^2)))+((2 Cos[30 t] Sin[30 x])/(1125 (π^2)))-((3 Cos[32 t] Sin[32 x])/(2048 (π^2))))},{x,-Pi/2,Pi/2}],{t,0,10}]

Kilka uwag ogólnych

Proszę zwrócić uwagę na różnicę pojęciową pomiędzy rozwiązywaniem problemu w całej przestrzeni a rozwiązywaniem w ograniczonym obszarze. 
1. W pierwszym przypadku nie musimy podawać warunków brzegowych. Skuteczne metody to np. rozwiązania ogólne i transformaty. 
2. W drugim przypadku musimy podać obszar w którym szukamy rozwiązanie i warunki brzegowe na końcach przedziału. Skuteczne metody to np. NDSolve i rozdzielanie zmiennych

Wyspecyfikowanie warunków brzegowych jest ,,równoważne'' pewnym ,,falom'' wchodzącym do obszaru z zewnątrz i ,,kasowaniu się'' fal. Rozważmy np. rozwiązanie Sin[x] Cos[t] które spełnia równanie falowe w całej przestrzeni:

rfalowe /. u->(Sin[#1]Cos[#2]&)

a jednocześnie spełnia warunek brzegowy u[-Pi,t]==u[Pi,t]==0, czyli opisuje drgania struny:

Manipulate[Plot[If[n>0,Sin[x]Cos[t],Sin[x]Cos[t] UnitStep[-(x-Pi)(x+Pi)]],{x,-2 Pi,2 Pi},PlotRange->{-1,1}],{t,0,10},{n,{0,1}}]