Algebra Symboliczna
Wykład XI

Andrzej Odrzywolek
Instytut Fizyki, Zakład Teorii Względnosci i Astrofizyki
21-01-2009, środa, 13:15

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/

Komentarz do wyników kolokwium

Zadanie 1.

Rozwiązanie:

uklad={I1+I2==I0,
I3+I4==I1+I2,
I1-I5==I3,I2+I5==I4,R I5 +R1 I1==R2 I2,R2 I3 - R I5 == R1 I4,R1 I1+R2 I3==U}

niewiadome={I0,I1,I2,I3,I4,I5}

Solve[uklad]

rozw=Solve[uklad,niewiadome]

U/I0 /. rozw[[1]]

Rz=Simplify[%] (* ogólne rozwiązanie *)

Rz/.R->0 (* R->0* )

Rz/.R->Infinity

Limit[Rz,R->Infinity] (* R->Infinity *)

Rz /. R2->R1

Simplify[%]

R1 (* R1=R2 *)

Sposób 2 (macierze od p. T. Trześniewskiego):

ℐ=Table[ToExpression["ℐ"<>ToString[i-1]],{i,1,6},{j,1,1}];
ℐ//MatrixForm

𝒜={{-1,1,1,0,0,0},{0,-1,-1,1,1,0},{0,1,0,-1,0,-1},{0,0,1,0,-1,1},{0,ℛ1,-ℛ2,0,0,ℛ0},{0,0,0,ℛ2,-ℛ1,-ℛ0},{0,ℛ1,0,ℛ2,0,0}};
𝒜//MatrixForm

Dimensions[𝒜]

ℬ={{0},{0},{0},{0},{0},{0},{𝒰}};
ℬ//MatrixForm

solution=LinearSolve[𝒜,ℬ]

𝒰/solution[[1,1]] //Simplify

ℐ//Flatten

Solve[𝒜.ℐ==ℬ,ℐ//Flatten]

Popełniane błędy:

1) pozostawienie wyniku w postaci nieuproszczonej; użycie Simplify drastycznie zmienia wartość odpowiedzi: powyżej widzimy, że Rz nie zależy od R (!) co nie jest oczywiste
2) przypisanie R1=R2 jest trwałe, i po jego wywołaniu wszystkie późniejsze obliczenia je uwzgledniają; lepiej użyć np. R1->R2
3) pochopne użycie %; użycie odpowiedniej nazwy redukuje możliwość pomyłki i konieczność przeliczania rachunku za każdym razem. Znak %   nie zawsze odnosi się do wyniku linijka wyżej.
4) błędy nie dotyczyły istoty problemu; wynika to być może ze złych przyzwyczajeń pochądzących z rachunków ręcznych: rozwiązanie układu równań na kartce stanowi 99% całej pracy jaką należy wykonać aby uzyskać wynik i nikt nie zwraca uwagi na pozostałe 1%, nawet gdy zawiera rażące błędy. W MATHEMATICE większość pracy to wpisywanie równań i manipulacja rozwiązaniami, czyli to na co zwykle nie kładzie się nacisku. Proszę pamiętać, że obliczenia wykonuje się w jakimś CELU, a wynik musi być poprawny.

Zadanie 2.

Rozwiązanie:

rownanie=r''[ϕ]+ϕ r'[ϕ]+r[ϕ]==Sin[ϕ]

rozwOgolne=DSolve[rownanie, r[ϕ],ϕ]

f1 = r[ϕ]/.rozwOgolne[[1]] //FullSimplify

Rozwiązanie, pomimo że zawiera jednostkę urojoną jest rzeczywiste dla rzeczywistego ϕ. Można to sprawdzić np. poprzez szereg potęgowy. Widać wtedy także, że rozwiązanie pomimo skomplikowanego wzoru zachowuje się w ϕ=0 jak funkcja liniowa. Z szeregu widać, że warunek r[0]==1 oznacza że C[2]=1.

Series[f1,{ϕ,0,5}]

Rozwiązujemy zagadnienie początkowe wg. treści:

rozwSzczegolne=DSolve[{rownanie,r[0]==1,r'[0]==α}, r[ϕ],ϕ]

Lub inaczej:

rozwOgolne2=rozwOgolne /. {C[1]->p, C[2]->q}

Clear[ϕ]

r1=((r[ϕ]/.rozwOgolne2[[1]])/.ϕ->0)==1

r2=(D[(r[ϕ]/.rozwOgolne2[[1]]),ϕ]/.ϕ->0)==α

staleCalkowania=Solve[{r1,r2}]

rozwSzczegolne2=rozwOgolne2 /. staleCalkowania

Ostatecznie nasze rozwiązanie szczególne zagadnienia poczatkowego to:

R=r[ϕ]/.rozwSzczegolne[[1]]

Rysujemy. Ponieważ nic nie wiadomo na temat spodziewanego zakresu α sprawdzam bardzo zakres α w postaci potęgowej od 0.001 do 10^3. Szukam takiej wartości α dla której przy ϕ=2 π (prawa krawędź wykresu) R ponownie staje się równe 1; jest to warunek konieczny aby krzywa była zamknięta.

Manipulate[Plot[{1,R /. α->znak 10^lgalpha},{ϕ,0,2 Pi},PlotRange->{-5,5},PlotLabel->"α = "<>ToString[znak 10.^lgalpha]],{lgalpha,-3,3},{znak,{-1,1}}]

Przybliżoną wartość z Manipulate wstawiamy do FindRoot:

rozw=FindRoot[(R /.ϕ->2 Pi)==1,{α,6.09}]

Rzamkniete=R /. rozw

Sprawdzenie:

wykr=PolarPlot[Rzamkniete,{ϕ,0,2 Pi},PlotStyle->Thick]

A teraz "nieszczęsne" pole figury  (http://en.wikipedia.org/wiki/Area#Area_in_calculus):

Całka : sposób 1

pole1=NIntegrate[Rzamkniete^2/2,{ϕ,0,2 Pi}] //Chop (* wzór na pole we wsp. biegunowych, łatwy do wyprowadzenia *)

pole2=NIntegrate[NIntegrate[r,{r,0,Rzamkniete}],{ϕ,0,2 Pi}] (* ogólny wzór we wsp. biegunowych *)

Całka : sposób 2

Podział na trójkąty:

trojkaty=Table[{Rzamkniete Cos[ϕ],Rzamkniete Sin[ϕ]}//Chop,{ϕ,0, 2Pi, 0.1}];

w=ListPlot[trojkaty,PlotRange->{{-4,4},{-2,6}}] (* Pole całego wyświetlonego prostokata to 64 ! Figura wewnątrz nie może mieć więcej !*)

Dimensions[trojkaty]

Wzór na pole trójkąta o zadanych 3 wierzchołkach w tym jeden w punkcie (0,0) : 1/2 Det[{a,b}] (http://pl.wikipedia.org/wiki/Trójkąt) lub identyczny wzór na pole wielokąta wypukłego (http://mathworld.wolfram.com/PolygonArea.html)

Sum[1/2Det[{trojkaty[[ii]],trojkaty[[ii+1]]}],{ii,1,Dimensions[trojkaty][[1]]-1}]+1/2Det[{trojkaty[[1]],trojkaty[[-1]]}]

Całka : sposób 3

A teraz metoda ,,Monte Carlo''. Losujemy pewną ilośc punktów na płaszczyżnie tak aby pokryły naszą ,,figurę'':

tblXY=Table[{RandomReal[{-4,4}],RandomReal[{-2,6}]},{ii,1,2000}];

Teraz musimy przeliczyć punkty do układu biegunowego. Teoretycznie robi to funkcja ArcTan , ale kąt ma zakres od -π do π, a my potrzebujemy od 0 do 2 π :

Plot[ArcTan[Cos[ϕ],Sin[ϕ]],{ϕ,0, 2 Pi}]

Dlatego trzeba nieco zmodyfikować funkcję ArcTan jak poniżej:

Plot3D[ArcTan[x,y]-Pi(Sign[y]-1),{x,-4,4},{y,-4,4},PlotPoints->50]

Teraz generujemy punkty we współrzędnych biegunowych:

tblϕR=Table[{ArcTan[tblXY[[ii,1]],tblXY[[ii,2]]]-Pi (Sign[tblXY[[ii,2]]]-1),Sqrt[tblXY[[ii,1]]^2+tblXY[[ii,2]]^2]},{ii,1,2000}];

Table[{tblϕR[[ii,2]]*Cos[tblϕR[[ii,1]]],tblϕR[[ii,2]]*Sin[tblϕR[[ii,1]]]},{ii,1,2000}]-tblXY//Tr (* sprawdzenie czy współrzędne są OK *)

Sprawdzenie graficzne:

ListPlot[tblXY] (* wsp. kartezjańskie *)

ListPolarPlot[tblϕR,PlotStyle->Red] (* wsp. biegunowe *)

Show[%,%%] (* punkty pokrywają się *)

UWAGA: naiwne losowanie bezpośrednio we wsp. biegunowych  w sposób poniżej prowadzi do poważnego błędu: rozkład nie jest równomierny, jest skoncentrowany wokół R=0:

Table[{RandomReal[{0,2 π}],RandomReal[{0,10}]},{ii,1,2000}] //ListPolarPlot (* ŹLE ! *)

Teraz potrzebujemy tylko wybrać punkty które znajdują się wewnątrz  krzywej:

Show[wykr,ListPolarPlot[tblϕR,PlotStyle->Red] ] (* w desperacji można policzyć ręcznie ...*)

Dla wygody definiujemy funkcję zwracającą promień naszej krzywej zamkniętej (metodą Ctrl+C Ctrl+V). Re na końcu usuwa numeryczne 0.0 I.

Clear[ϕ]

Rz[ϕ_]:=(1/4) (E^(-((ϕ^2)/2))) (4+(30.512172793263105`+0.` I) Erfi[(ϕ/(Sqrt[2]))]+2 (Sqrt[2 π]) Erfi[(ϕ/(Sqrt[2]))]-(Sqrt[2 E π]) Erfi[((-I+ϕ)/(Sqrt[2]))]-(Sqrt[2 E π]) Erfi[((I+ϕ)/(Sqrt[2]))])//Re

Rz[1]

Show[wykr,ListPolarPlot[Select[tblϕR,#[[2]]<Rz[#[[1]]]&],PlotStyle->Red]]

Liczymy ilość punktów wewnątrz:

in=Select[tblϕR,#[[2]]<Rz[#[[1]]]&]//Dimensions //First

all=tblϕR//Dimensions//First (* liczba wszystkich punktów *)

in/all*64.0 (* pole całego kwadratu wynosi 64 *)

Całka : sposób 4

A teraz obliczymy 4 zwykłe całki: po jednej dla każdej ćwiartki.

funkcja1=Table[{Rz[ϕ]Sin[ϕ],Rz[ϕ]Cos[ϕ]},{ϕ,0,Pi/2,0.5}]//Interpolation

Plot[funkcja1[y],{y,0,5.218}]

pole1=NIntegrate[funkcja1[y],{y,0,5.218}]

funkcja2=Table[{Rz[ϕ]Sin[ϕ],-Rz[ϕ]Cos[ϕ]},{ϕ,Pi/2,Pi,Pi/8}]//Interpolation

Plot[funkcja2[y],{y,0,5.165}]

pole2=NIntegrate[funkcja2[y],{y,0,5.165}]

funkcja3=Table[{Rz[ϕ]Cos[ϕ],-Rz[ϕ]Sin[ϕ]},{ϕ,Pi,2 Pi,Pi/16}]//Interpolation

Plot[funkcja3[x],{x,-2.92,1}] (* coś jest nie tak: oczywisty błąd interpolacji !*)

funkcja3=Interpolation[Table[{Rz[ϕ]Cos[ϕ],-Rz[ϕ]Sin[ϕ]},{ϕ,Pi,2 Pi,Pi/16}],InterpolationOrder->1] (* redukujemy interpolację do liniowej *)

Plot[funkcja3[x],{x,-2.92,1}]

pole3=NIntegrate[funkcja3[x],{x,-2.92,1}]

pole1+pole2+pole3 (* wynik *)

Całka : sposób 5

Z twierdzenia Stokesa (zamiana całki krzywoliniowej na powierzchniową) :

-NIntegrate[Rzamkniete Sin[ϕ] D[Rzamkniete Cos[ϕ],ϕ],{ϕ,0,2 Pi}] //Re

NIntegrate[Rzamkniete Cos[ϕ] D[Rzamkniete Sin[ϕ],ϕ],{ϕ,0,2 Pi}]//Re

Podręcznikowa postać tego wzoru (http://en.wikipedia.org/wiki/Area#Area_in_calculus):

1/2 NIntegrate[Rzamkniete Cos[ϕ] D[Rzamkniete Sin[ϕ],ϕ]-Rzamkniete Sin[ϕ] D[Rzamkniete Cos[ϕ],ϕ],{ϕ,0,2 Pi}] //Chop

Całka : sposób 6

Metoda dla desperatów: nakładamy siatkę i liczymy kwadraciki:

Show[wykr,GridLines->{Table[ii,{ii,-4,4,1/2}],Table[ii,{ii,-2,6,1/2}]},PlotRange->{{-4,4},{-2,6}}]

Na zewnątrz: 94, wewnątrz 103, całkowita ilość: 256, pole kratki : 0.25, wynik: 25.75<pole<40.5. Czas liczenia: 3 minuty.

Po uśrednieniu wynik jest całkiem sensowny:

Mean[{25.75,40.5}]

GeometricMean[{25.75,40.5}]

RootMeanSquare[{25.75,40.5}]

Inny prosty sposób oszacowania pola niżej: 21<pole<44. Jak widać wyniki typu 18 czy 66 są ewidentnie błedne !

Manipulate[Show[wykr,ParametricPlot[{r Cos[α]+x0,r Sin[α]+y0} ,{α,0,2 Pi}],PlotLabel->N[Pi r^2]],{{r,1},0.1,10},{x0,-2,2},{y0,0,4}]

Popełniane błędy (historia obliczania pewnej całki czyli jak uzyskać 87.5% całkowicie błędnych wyników, a każdy inny):

1) Lekcja z obliczania całki jest następująca: MATHEMATICA nie pomoże jak nie opanujemy matematyki...

2) Aby wykorzystać potencjał MATHEMATICI musimy znać bardzo wielu róznych metod zapisania rozwiązania problemu w postaci matematycznej. Albo przynajmniej wiedzieć o istnieniu metod, odpowiednie wzory i twierdzenia nie muszą tkwić w głowie, można korzystać z wikipedii, http://mathworld.wolfram.com, podręczników, notatek itp

3) Odrobina zdrowego rozsądku nie zaszkodzi: pole figury która w całości mieści się w kwadracie o wymiarach 8x8=64 nie może wynosić 66 podobnie nie moze wynosić 18 bo w figurę da się wpisać okrąg o polu 21 i czworokat o polu 26.5.

4) Pole figury musi mieć odpowiedni wymiar; jeżeli mierzymy R w np. metrach to pole musi mieć wymiar m^2. Całka
typu Integrate[R, {ϕ,0, 2 Pi}] nie może być polem ! Musi pojawić się iloczyn  dwóch wielkości o wymiarze długości, np. x dx, r^2, x y' itp.

5) Nawet zgrubne oszacowanie ma pewną wartość o ile jest poprawne, np. pole czworokąta wpisanego w figurę wynosi 26.68, a okręgu w którym się mieści 44.0. Dlatego na pewno 26.68<pole <44.0 ! Odpowiedź, że pole≈25 nie jest poprawnym ani użytecznym oszacowaniem,
bo nie wiadomo o co chodzi gdyż:  a) nie znamy definicji  ,,przybliżenia''  (sposobu obliczenia) b) dokładności tej liczby. Jeżeli przyjmiemy konwencję stosowaną w fizyce, że zapis pole=25 oznacza pole=25±1 na poziomie 1σ to oszacowanie jest fatalne z błędem 8σ.

6) Nie było wielkim wysiłkiem ustalenie, że przynajmniej część wyników jest całkowicie błędna. Wystarczyło wykonać bardzo grube oszacowanie, czyli policzyć pole kwadratu i  4 trójkatów lub 2 kół. To nie tylko pozwala wykluczyć pewne ,,metody obliczeniowe'',
ale przy okazji daje jakikolwiek poprawny wynik. Oczywiście studenci lekką ręką podają błędne wyniki, gdyż nic za to nie grozi, ale w realnych zastosowaniach  jest to katastrofalna sytuacja. Sytuację mogłyby zmienić ujemne punkty na kolokwium za podanie błędnego wyniku.

Zadanie 3.

Rozwiązanie:

rown=x^n==Exp[x^m]

Niech dla przykładu wymiar macierzy wynosi 4, czyli indeksy zmieniają się od 0 do 3. Rachunek wykonałem bez problemu dla num=13. Zaczynając od n=6 pojawiaja się przypadki z 4 rozwiązaniami rzeczywistymi.

num=6

W MATHEMATICE 7.0 zadanie jest szczególnie proste gdyż Reduce znajduje wszystkie rozwiązania korzystając z nowej postaci funkcji Root:

Table[Reduce[x^n==Exp[x^m]&&x∈Reals],{n,0,num},{m,0,num}]//MatrixForm

Ponieważ macierz nie jest bardzo duża, można policzyć ilość rozwiązań ,,ręcznie'' i wprowadzić macierz. False oznacza zero rozwiązań. Poniżej przykład automatycznego wypełnienia macierzy. Instrukcja If sprawdza czy pojawia się False; jeżeli tak to wstawia zero, w przeciwnym wypadku liczy znaki "x" których ilość musi być równa ilości rozwiązań.

M=Table[If[ToString[Reduce[x^n==Exp[x^m]&&x∈Reals]] =="False",0, StringCount[ToString[Reduce[x^n==Exp[x^m]&&x∈Reals]],"x"]],{n,0,num},{m,0,num}]

MatrixForm[M]

Pozostałe obliczenia są kwestią wywołania odpowiednich komend na macierzy M:

Det[M] (* wyznacznik *)

Tr[M] (* ślad *)

M //Flatten//Tr (* suma wszyskich elementów *)

Sum[ M[[i,j]],{i,1,num+1},{j,1,num+1}] (* inny sposób na sumę elementów *)

Sum[ M[[i,j]]^2,{i,1,num+1},{j,1,num+1}](* suma kwadratów wszyskich elementów *)

M^2 //Flatten//Tr (* suma kwadratów wszyskich elementów *)

ww=Eigenvalues[M]  (* wartości własne *)

ww //N (* numeryczne wartości wart. własnych *)

Powyżej wyraźnie widać, które wart. własne są rzeczywiste, ale możemy je wybrać także odpowiednią komendą:

Select[ww,Im[#]==0&] (* wybierz te których część urojona jest zero *)

N[%]

Użycie Solve daje inny (błędny) rezultat:

Table[Solve[x^n==Exp[x^m]],{n,0,num},{m,0,num}]//N//MatrixForm //TraditionalForm

Mzespolone=Table[Dimensions[Solve[x^n==Exp[x^m]]][[1]],{n,0,num},{m,0,num}] (* ilość rozw. zespolonych wg. Solve *)

Table[Select[Table[ x/.Solve[x^n==Exp[x^m]][[inum]],{inum,1,Mzespolone[[n+1,m+1]] }],Im[#]==0&] //Dimensions//First,{n,0,num},{m,0,num}]//MatrixForm //TraditionalForm (* ZŁY WYNIK !!! ilość rozw. rzeczywistych wg. Solve *)

M//MatrixForm//TraditionalForm

Sposób graficzny. Uzyskanie 2 rózych macierzy poprzez Solve i Reduce rodzi wątpliwości czy wyniki te są poprawne. Jedyny sposób aby się o tym przekonać to narysowanie wykresu. Użycie wykresu LogPlot pozwala na wizualizację dużego zakresu zmienności, co gwarantuje że nie przeoczymy przecięcia dla dużych wartości "x". Ustawienie n=2 (3 wiersz) i m=3 (4 kolumna) ewidentnie pokazuje, że rozwiązanie rzeczywiste istnieje. Solve daje więc błedny wynik, przed czym zresztą ostrzega sugerując użycie Reduce.

Manipulate[LogPlot[{Tooltip[x^n,"(x^n)"],Tooltip[Exp[x^m],"Exp[(x^m)]"]},{x,-10,10},PlotRange->{10^-4,10^4}],{n,0,num,1},{m,0,num,1}]

Uwagi:

1) Zadanie wymagało głównie staranności i dokładności, gdyż w gruncie rzeczy sprowadzało się do policzenia przecięć 2 krzywych dla 16 kombinacji parametrów n i m. Godne pochwały jest zapisywanie częściowych wyników na kartce, co bardzo ułatwiało zadanie, w szczególności zmniejszało prawdopodobieństwo pomyłki. 
2) Automatyzacja obliczeń, jak w rozwiązaniu wyżej  nie jest niezbędna aby rozwiązać to zadanie. Wystarczyło narysować 16 wykresów i wypełnić ręcznie macierze.
3) Praktyka rysowania wykresów w skali liniowej prowadzi czasem do przeoczenia przecięć które znajdują się poza zakresem. Najpewniejsze rezultaty daje użycie metod symbolicznych i graficznych równocześnie, co pozwala na porównanie ilości rozwiązań. W przypadku równań zawierających funkcje Exp i Power czy innych przyjmujacych duże wartości skuteczny jest LogPlot