Algebra Symboliczna
Wykład X

Andrzej Odrzywolek
Instytut Fizyki, Zakład Teorii Względnosci i Astrofizyki
14-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/

Zagadnienie własne

Definicja i obliczanie wartosci i wektorów własnych: dla dowolnej macierzy kwadratowej (NxN) równaniem własnym nazywamy równanie:
𝒜·U==λ U
gdzie U jest N-elementowym wektorem a λ liczbą zespoloną. Wektor bedący rozwiązaniem powyższego równania dla pewnej wartosci λ nazywamy wektorem własnym. Na ogół rozwiazanie istnieje tylko dla pewnych wybranych wartosci λ które nazywamy wartosciami własnymi. Rozwiazanie opisanego wyzej problemu jest jednym z podstawowych problemów algebry i ma ogromną ilosc zastosowań, m. in. w fizyce oraz metodach numerycznych..

Prosty przykład:

P={{1,2},{2,1}}

V={Vx,Vy}

Do rozwiązania równania własnego możemy użyć Solve lub Reduce, ale należy uważać, gdyż rozwiązanie może być niekompletne:

Solve[P.V == λ V,V]

Solve[P.V == λ V,{λ,Vx,Vy}]

rozw=Solve[P.V==λ V]

P.V == λ V /. rozw[[1]]

P.V == λ V /. rozw[[2]]

P.V == λ V /. rozw[[3]]

Reduce[P.V==λ V]

Z powyższego widać, że istnieją 3 rozwiązania: V==0 dla dowolnego λ, oraz dwie wartości własne λ==-1 i λ==3 z odpowiadającymi im wektorami własnymi. Istnieje nieskończenie wiele wektorów własnych dla każdej z wartości własnych, różniacych się długością. Dlatego zwykłe podaje się tylko jeden konkretny wektor własny np:

λ1=-1

V1={1,-1}

P.V1==λ1 V1

λ2=3

V2={1,1}

P.V2

λ2 V2

Czasem potrzebujemy wektorów znormalizowanych, np. do jedności:

V1a=V1/Sqrt[V1.V1]

V1a.V1a

Przepisując równanie własne w postaci poniżej widać ścisły związek zagadnienia własnego z problem rozwiązania układu równań liniowych:

(P-λ IdentityMatrix[2]).V==0

MatrixForm[%[[1]]]

Taki (jednorodny) układ równań posiada niezerowe rozwiązania tylko gdy wyznacznik jest równy zero:

Det[P-λ IdentityMatrix[2]]==0

Solve[%]

Otrzymany w powyższy sposób wielomian nazywamy wielomianem charakterystycznym; można go obliczyć poleceniem CharacteristicPolynomial

CharacteristicPolynomial[P,λ]

Wartości własne jak (wektory własne) możemy policzyć dedykowaną instrukcją Eigenvalues (Eigenvectors)

Eigenvalues[P]

Eigenvectors[P]

Numeryczne vs symboliczne obliczanie wartości i wektorów własnych

Dla macierzy w postaci symbolicznej (symbole literowe, liczby całkowite, wymierne) rozwiązanie równania własnego, wielomianu charakterystycznego i wyznaczanie wartości własnych poprzez Eigenvalues niczym praktycznie się nie różni. Jeżeli macierz jest dana  w postaci numerycznej (zmiennoprzecinkowej) to Eigenvalues obliczy ww. korzystając z wysoko wydajnego algorytmu numerycznego,  do daje przyspiesznie obliczeń o kilkanaście rzędów wielkości.

macierze=Table[Table[RandomInteger[{-10,10}],{i,1,n},{j,1,n}],{n,1,199}];

symb=Table[N[Eigenvalues[macierze[[ii]]]]//AbsoluteTiming //First,{ii,1,80}]

ListLogLogPlot[symb]

Eigenvalues[N[macierze[[80]]]]//Timing

num=Table[Eigenvalues[N[macierze[[ii]]]]//AbsoluteTiming //First,{ii,1,80}]

ListPlot[num]

num=Table[Eigenvalues[N[Table[RandomInteger[{-10,10}],{i,1,n},{j,1,n}]]]//AbsoluteTiming //First,{n,1,299}]

ListLogPlot[num,Joined->True,InterpolationOrder->0,PlotRange->All]

Z punktu widzenia obliczeń numerycznych pewne rozmiary macierzy są preferowane, dlatego czasem lepiej macierz powiększyć aby szybciej dostać wynik. Poniższy przykład pokazuje jawnie że takie ,,rozszerzenie'' nie zmiwnia wartości własnych:

A=ArrayPad[{{a,b},{c,d}},{0,2}]

Eigenvalues[{{a,b},{c,d}}]

Eigenvalues[A]

Rozmiar macierzy NxN której MATHEMATICA nie poradzi sobie z obliczeniami symbolicznymi to N=1-100, w niektórych przypadkach problematyczne mogą być rachunki już dla N=3-10.  Numerycznie typowy nowoczesny (32-bitowy) komputer PC powinien bez problemu poradzić sobie z macierzami o rozmiarze do N=5000.

Eigenvalues[N[Table[RandomInteger[{-10,10}],{i,1,5000},{j,1,5000}]]]//AbsoluteTiming //First

N[Sqrt[Sqrt[8769983410567396]],1000]

Zastosowanie wartości i wektorów własnych: diagonalizacja, dekompozycja i transformacje macierzy

Własności macierzy diagonalnych

n=9

Clear[λ1,λ2,V1,V2]

Macierz diagonalna posiada różne od zera elementy wyłącznie na przekątnej:

A=DiagonalMatrix[Table[ToExpression["λ"<>ToString[ii]],{ii,1,n}]]

A//MatrixForm

B=DiagonalMatrix[Table[ToExpression["δ"<>ToString[ii]],{ii,1,n}]]

V=Table[ToExpression["V"<>ToString[ii]],{ii,1,n}]

Większość typowych operacji macierzowych jest szczególnie prosta w postaci diagonalnej:

A.B //MatrixForm (* Iloczyn macierzy *)

A.V

Inverse[A] //MatrixForm (* Macierz odwrotna *)

Det[A] (* Wyznacznik *)

Tr[A]

MatrixPower[A,k]//MatrixForm

W sposób identycznye jak poniżej dla eksponenty można obliczyć dowolną funkcję od macierzy diagonalnej, t.j. obliczając funkcję od elementów na przekątnej i pozostawiając zera. Przepis ten jest oczywisty jeżeli posłużymy się szeregiem potęgowym macierzy i popatrzymy na macierz A^k.

MatrixExp[A]//MatrixForm

Obliczmy np. MatrixSin[A]

A//Diagonal//Sin//DiagonalMatrix

Eigenvalues[A]

Diagonalizacja macierzy : przykład zastosowania

Jako przykład obliczymy e^P, MatrixSin[P] oraz MatrixLog[P] ,,ręcznie'' dokonując diagonalizacji.

P //MatrixForm

Tworzymy 2 macierze: macierz ,,zdiagonalizowaną'' Diag (zawierającą na przekątnej wartości własne)  i macierz transformacji U składającą się z wektorów własnych:

{WW,U}=Eigensystem[P]

Diag=DiagonalMatrix[WW]

Zachodzą ponizsze tożsamości, opisującą  transformację macierzy do i z postaci diagonalnej:

U.P.Inverse[U]==Diag

Inverse[U].Diag.U==P

Idea rachunków przeprowadzanych za pomocą diagonalizacji jest następująca:
1. Obliczamy wartości własne i wektory własne
2. Z wektorów własnych tworzymy macierz transformacji U, a z wartości własnych macierz diagonalną Diag
3. Wykonujemy interesującą nas operację na macierzy Diag
4. Dokonujemy transformacji do układu wyjściowego wzorem Inverse[U].wynik.U
Jak widać poniżej, procedura ta dla eksponenty daje oczekiwany wynik:

Inverse[U].(Diag//Diagonal//Exp//DiagonalMatrix).U

MatrixExp[P]

%-%%

Simplify[%]

Obliczmy sinus macierzy P:

Diag //Diagonal //Sin//DiagonalMatrix

sinP1=Inverse[U].%.U

Inny sposób: zamiana Sin na Exp

Sin[x]//TrigToExp

sinP2=(MatrixExp[I P]-MatrixExp[-I P])/2/I//Simplify

FullSimplify[sinP1-sinP2]

sinP2//N

Sin poprzez szereg potęgowy:

sinP3=Sum[(-1)^k MatrixPower[P,2k+1]/(2k+1)!,{k,0,Infinity}]

sinP2-sinP3//Simplify

Podobnie można policzyć logarytm czy pierwiastek z macierzy

logP=Inverse[U].(Diag//Diagonal//Log//DiagonalMatrix).U

sqrtP=Inverse[U].(Diag//Diagonal//Sqrt//DiagonalMatrix).U

MatrixExp[logP] //MatrixForm

sqrtP.sqrtP //Simplify//MatrixForm

Czy nawet pierwiastek 3-tego stopnia lub P^P:

sqrt3P=Inverse[U].(Diag//Diagonal//(#^(1/3)&)//DiagonalMatrix).U //ExpToTrig

sqrt3P.sqrt3P.sqrt3P//Simplify

MatrixExp[P.logP] (* P^P *)

Analiza wektorowa

Do obliczeń związanych z analizą pól wektorowych i skalarnych w krzywolniowych układach współrzędnych potrzebujemy pakietu:

Needs["VectorAnalysis`"]

Definiujemy układ współrzędnych:

SetCoordinates[Cartesian[x,y,z]]

Cartesian[x,y,z]

Definujemy pole wektorowe:, np.pole magnetyczne:

B={Bx[x,y,x],By[x,y,z],Bz[x,y,z]}

Definiujemy pole skalarne:

ϕ=Φ[x,y,z]

Oto standardowe operacje na polach wektorowych i skalarnych we współrzędnych kartezjańskich: rotacja, dywergencja i gradient

Curl[B] (* rotacja *)

Div[B] (* dywergencja *)

Grad[ϕ] (* gradient *)

Div[Curl[B]]

Div[{x,y,z}]

1/2Curl[{-y,x,0}] (* prędkość rotacji  *)

Curl[Grad[ϕ]]

Laplacian[ϕ] (* laplasjan pola skalarnego *)

Laplacian[B](* laplasjan pola wektorowego *)

V=Re[Log[Sin[(x+I y)^3]]] //ComplexExpand

Laplacian[V]

Simplify[%] (* laplasjan części rzeczywistej i urojonej f. analitycznej wynosi zero *)

Bardziej interesujące są obliczenia we współrzednych krzywolinowych. Proszę bardzo uważać na kolejnośc pojawiania się zmiennych! Może być inna niż typowo w polskich podręcznikach.

SetCoordinates[Spherical] (* wbudowane nazwy *)

Clear[ϕ]

SetCoordinates[Spherical[r, θ,ϕ]]

V=-1/r

Grad[V]

Div[Grad[V]]

Div[{r,0,0}]

Przykład: ,,czysta rotacja'' w układzie cylindrycznym

SetCoordinates[Cylindrical[r,θ,z]]

Pole wektorowe rotującego w dowolny sposób płynu:

V={0,r Ω[r,z],0}

Tzw. twierdzenie Poincare-Wavre wynikające z równań Eulera (ruchu płynu bez lepkości)

pw=Curl[Cross[V,Curl[V]]]

DSolve[pw[[2]]==0,Ω[r,z], {r,z}]

Siła odśrodkowa w poruszającym się płynie:

Cross[V,Curl[V]]-1/2Grad[V.V]//Simplify

Pole molekularne w ciekłym krysztale:

Clear[n]

n={Cos[f[r]],0,Sin[f[r]]}

Curl[n]

Div[n]

H = K1 Grad[Div[n]]-K2((n.Curl[n]) Curl[n]+Curl[(n.Curl[n]) n])+K3 (Cross[Cross[n,Curl[n]],Curl[n]]+Curl[Cross[n,Cross[n,Curl[n]]]])

Równanie równowagi ciekłego kryształu:

H-n n.H

Simplify[%]

rnematyk=% /. K3->K1

rnematyk /.f[r]->0

DSolve[{K1 Sin[2 f[r]]+2 K1 r f'[r]+2 K1 (r^2) f''[r]==0},f[r],r]

% /. C[1]->p /. C[2]->q

f[r] /. %[[1]]

sol=ϕ==%

Manipulate[Plot[-JacobiAmplitude[(Sqrt[1+2 p]) q-(Sqrt[1+2 p]) Log[r],(1/(1+2 p))],{r,0,1}],{p,-1,1},{q,0,1}]

-JacobiAmplitude[(Sqrt[1+2 p]) q-(Sqrt[1+2 p]) Log[r],(1/(1+2 p))]

N[%]