Algebra Symboliczna
Wykład IX

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

Układ równań liniowych

Rozważmy układ równań:

ur={x+y+z+t==1,2x+4y-t==2,x-y+2t==0,x+y-z-t==-1}

Który możemy rozwiązać korzystając z już poznanych metod:

Solve[ur,{x,y,z,t}]

Reduce[ur,{x,y,z,t}]

FindInstance[ur,{x,y,z,t}]

FindRoot[ur,{{x,1},{y,1},{z,1},{t,1}}]

Ze względu na ogromne znaczenie teoretyczne i praktyczne do rozwiązywania takich układów stosuje się rachunek macierzowy. Ze współczynników liczbowych przy zmiennych x,y,z,t tworzymy macierz o rozmiarze 4x4. Wprowadzamy ją jako listę złozoną z 4 list:

A={{1,1,1,1},{2,4,0,-1},{1,-1,0,2},{1,1,-1,-1}}

Jak widać poniżej, kolejne 4-elementowe listy stanowią wiersze macierzy:

%//MatrixForm

Z wyrazów wolnych tworzymy wektor równoważny w MATHEMATICE liście jednowymiarowej:

B={1,2,0,-1}

%//MatrixForm

Wektor można mnożyć przez macierz lewo- i prawostronnie. Mnożenie macierzowe wykonujemy kropką:

B.A

A.B

Tworzymy wektor X  zbudowany z niewiadomych

X={x,y,z,t}

Teraz układ równań możemy zapisać jako:

A.X==B

Nadal możemy rozwiązać przez Solve (to może zależeć od wersji MATHEMATICI):

Solve[A.X==B, X]

Teraz rozwiązanie mozemy znaleźć metodami rachunku macierzowego.  Inverse oblicza macierz odwrotną:

Inverse[A].B

Identyczny rezultat otrzymamy za pomocą LinearSolve. Funkcja ta jest efektywna dla bardzo dużych macierzy.

LinearSolve[A,B]

Poniżej przykład dla większych macierzy:

Table[RandomReal[],{ii,1,500},{jj,1,500}];

A=%+Transpose[%];

B=Table[RandomReal[],{ii,1,500}];

Inverse[A].B; //Timing

LinearSolve[A,B] ;//Timing

Wektory, macierze, wektory macierzowe i listy

Przy pracy ze sformalizowanym systemem, jak MATHEMATICA musimy pamiętać o pewnych sprawach które zwykle traktuje się domyslnie w rachunkach ręcznych. Obiekty takie jak wektor, macierz 1xN czy Nx1, macierz o wymiarze 1x1 oraz liczba zespolona nie są tymi samymi obiektami!

Rozważmy dwa wektory macierzowe: poziomy "a" i pionowy "b":

a={{x,y,z}}

MatrixForm[a]

b={{X},{Y},{Z}}

MatrixForm[b]

Ich iloczym daje w wyniku macierz i wymiarze 1x1 :

a.b

Wynik powyżej nie jest równy liczbie będącej jej jedyną skladową, ponieważ są to zupełnie inne obiekty:

a.b==x X+y Y+z Z //Simplify

W rachunkach ręcznych operację poniżej wykonujemy w głowie pisząc po prostu a·b = x X + y Y + z Z. Formalnie jednak równośc zachodzi pomiędzy składową macierzy 1x1  a liczbą, a nie samą macierzą i liczbą:

(a.b)[[1,1]]==x X+y Y+z Z //Simplify

Rozważmy teraz wektory :

u={x,y,z}

v={X,Y,Z}

Ich iloczyn skalarny nie zależy od kolejności i wynosi:

u.v

v.u

v.u==u.v

Natomiast dla wektorów macierzowych, nie tylko wynik jest inny, ale rezultatem jest zupełnie inny obiekt: macierz 1x1 dla a.b i macierz 3x3 dla b.a:

a.b //MatrixForm

b.a //MatrixForm

Podstawy rachunku macierzowego

• Dodawanie i odejmowanie macierzy (i innych list o zgodnch wymiarach) odbywa się za pomocą zwykłych znaków "+" i "-":

A={{a11,a12},{a21,a22}}

B={{1,-1},{1,-1}}

A+B

A-B

Prosze zwrócić uwagę na wynik ponizszego odejmowania. Jest nią macierz zerowa o wymiarze 2x2. MATHEMATICA rozróżnia te zera, podczas gdy w rachunkach ręcznych zwykle pisze się po prostu "0".

A-A

Dlatego poniższe wyrażenia nie są prawdziwe:

A-A==0

A-A=={0}

W przypadku liczb zespolonych jest nieco inaczej:

0+0I==0

Ale ...

Infinity+I Infinity==Infinity

• Działania na pojedynczych elementach macierzy i całych macierzach

Nie ma takiej operacji jak dzielenie macierzy. Wynikiem tego typu operacji, podobnie jak wielu innych jest macierz w której każdy element macierzy A został podzielony przez odpowiadający mu element macierzy B. Pomylenie działania macierzowego z działaniem na elementach listy to najczęstszy błąd studentów!

A/B

A/A

B/B

• Transpozycja (przestawienie wierszy z kolumnami)

U={{1,x},{y,1}}

MatrixForm[U]

Transpozycja powoduje, że wiersze stają się kolumnami i na odwrót. Dla macierzy kwadratowej jest to równoważne odbiciu względem przekątnej.

Transpose[U]//MatrixForm

Wektora nie można transponować:

Transpose[{a,b,c,d,e}]

Natomiast macierz 1xn można:

{{a,b,c,d,e}}//MatrixForm

Transpose[%] //MatrixForm

Macierz symetryczną tworzymy w następujący sposób dodając do siebie macierz i macierz transponowaną:

tmp=RandomReal[1,{500,500}];

gauss=RandomReal[NormalDistribution[],{500,500}];

T=(tmp+Transpose[tmp])/2;

W=(gauss+Transpose[gauss])/2;

Macierz symetryczna posiada same rzeczywiste wartości własne:

wwT=Eigenvalues[T]//Sort //Short

Dla macierzy symetrycznej o elementach losowych z rozkładu jednorodnego na odcinku (0,1) (RandomReal[] ) największa wartość własna jest w przybliżeniu równa 1/2 wymiaru macierzy:

Max[wwT] //Round

Min[wwT]

wwW=Eigenvalues[W]//Sort;

Max[wwW]

Nie licząc największej wartości własnej dla rozkładu jednorodnego, rozkład prawdopodobieństwa pojawienia się pewnej wartości własnej jest opisany przez półokrąg Wignera:

ListPlot[{BinCounts[Drop[wwT,-1],1], BinCounts[wwW,1]} ,InterpolationOrder->0,Joined->True]

• Mnożenie macierzy

Poniższe działanie nie jest mnożeniem macierzy i działa podobnie jak dzielenie powyżej:

A*B

Prawidłowe mnożenie macierzy odbywa się za pomocą KROPKI:

A.B

A B -A. B (* mnożenie ,,element przez element'' minus mnożenie macierzy przez macierz *)

W szczególności mnożenie macierzowe nie jest przemienne:

A B-B A (* ŹLE *)

A.B -B.A (* DOBRZE *)

• Potęgowanie macierzy dokonujemy poprzez mnożenie lub komendą MatrixPower:

B^3 (* ŹLE *)

B.B.B(* DOBRZE *)

MatrixPower[B,3](* DOBRZE *)

• Obliczanie macierzy odwrotnej wykonujemy poleceniem Inverse.

Inverse[A]

Nie istnieje macierz odwrotna do macierzy osobliwej, czyli o wyznaczniku równym zero:

Inverse[B]

Det[B]

Macierz  odwrotna ma sens tylko dla macierzy kwadratowych (NxN):

Inverse[B[[1]]]

A.Inverse[A]==Inverse[A].A

A.Inverse[A]

Simplify[%]

Inverse[A]. A==IdentityMatrix[2]

Simplify[%]

Inverse[Inverse[A]] //Simplify

• Funkcje z argumentem w postaci macierzy.

Eksponenta macierzy:

MatrixExp[ {{0,ϕ},{-ϕ,0}}]

Exp[{{0,ϕ},{-ϕ,0}}] (* ŹLE *)

W postaci diagonalnej funkcja macierzowa jest równoważna zwykłej funkcji elementów na przekątnej:

ϕ={{α,0,0},{0,β,0},{0,0,γ}}

MatrixForm[ϕ]

(MatrixExp[I ϕ]+MatrixExp[-I ϕ])/2 //ComplexExpand //MatrixForm

Cos[ϕ]//MatrixForm

Solve[Sinh[x]==1,x]

I jeszcze ciekawy przykład z manuala:

Manipulate[FullSimplify[MatrixExp[ x IdentityMatrix[n]//Reverse]]//MatrixForm,{{n,8},1,30,1}]

• Obliczanie wyznacznika, śladu, wartości i wektorów własnych

Tworzymy macierz n x n:

n=2

A=Table[ToExpression["a"<>ToString[i]<>ToString[j]],{i,1,n},{j,1,n}]

Wyznacznik obliczamy poleceniem Det natomiast ślad poleceniem Tr . Warto zwrócić uwagę, że w tym wypadku stosuje się skróconą nazwę polecenia. W szczególności polecenie Trace nie ma nic wspólnego z macierzami; służy do śledzenia wywoływanych funkcji.

Det[A] (* wyznacznik *)

Tr[A] (* ślad *)

Trace[A] (* ŹLE !!!! *)

Użyteczna tożsamość:

Det[MatrixExp[A]]-Exp[Tr[A]]//FullSimplify

Wartości własne obliczamy poleceniem Eigenvalues:

Eigenvalues[A]

A wektory własne poleceniem Eigenvectors

Eigenvectors[A]//MatrixForm

Wektory własne i odpowiednie wartości własne:

Eigensystem[A]

Macierze Pauliego i przegląd ich własności:

σ={IdentityMatrix[2],{{0,1},{1,0}},{{0,-I},{I,0}},{{1,0},{0,-1}}}

Table[MatrixForm[σ[[i]]],{i,1,4}]

Tabelka mnożenia:

Table[σ[[i]].σ[[j]],{i,1,4},{j,1,4}] //TraditionalForm

Wartości własne:

Table[Eigenvalues[σ[[i]]]//MatrixForm,{i,1,4}]

Wyznaczniki:

Table[Det[σ[[i]]],{i,1,4}]

Ślady:

Table[Tr[σ[[i]]],{i,1,4}]

Macierze Pauliego są odwrotne do samych siebie:

Table[Inverse[σ[[i]]]//MatrixForm,{i,1,4}]

Table[Inverse[σ[[i]]]==σ[[i]],{i,1,4}]

Unitarność:

Table[Inverse[σ[[i]]].σ[[i]]//MatrixForm,{i,1,4}]

σ^2 to macierz jednostkowa:

Table[σ[[i]].σ[[i]]//MatrixForm,{i,1,4}]

Bezpośrednie rachunki symboliczne

Dowodzenie tożsamości macierzowych bezpośrednim rachunkiem

Niektóre z twierdzeń dotyczących algebry liniowej można udowodnić dla macierzy o określonym rozmiarze bezpośrednim rachunkiem bez zakładania czegokolwiek na temat samej macierzy.

n=3

A=Table[ToExpression["a"<>ToString[i]<>ToString[j]],{i,1,n},{j,1,n}];

MatrixForm[A]

B=Table[ToExpression["b"<>ToString[i]<>ToString[j]],{i,1,n},{j,1,n}];

Unprotect[C]

C=Table[ToExpression["c"<>ToString[i]<>ToString[j]],{i,1,n},{j,1,n}];

A.(B.C)==(A.B).C //FullSimplify

A.(B+C)==A.B+A.C //FullSimplify

(B+C).A==B.A+C.A //FullSimplify

Transpose[C.B.A]==Transpose[A].Transpose[B].Transpose[C] //Simplify

UWAGA: wyrażenie poniżej szybko prowadzi do wyczerpania pamięci

MemoryConstrained[Inverse[A.B.C]==Inverse[C].Inverse[B].Inverse[A] //Simplify,1024^3] (* 1 GB *)

Dowód twierdzenia Cayleya-Hamiltona dla dowolnej macierzy o zadanym rozmiarze

Wielomianem charakterystycznym nazywamy wyrażenie:

n=6

A=Table[ToExpression["a"<>ToString[i]<>ToString[j]],{i,1,n},{j,1,n}]

poly=Det[A-λ IdentityMatrix[n]] ;

Wspólczynnik przy (λ^(n-1)) to ślad ze znakiem (((-1))^(n+1)):

Coefficient[poly, λ,n-1]==(-1)^(n+1) Tr[A]

Natomiast wolny wyraz to wyznacznik:

Coefficient[poly,λ,0]==Det[A] //Simplify

Możemy użyc też funkcji CharacteristicPolynomial

CharacteristicPolynomial[A,λ];

Twierdzenie Cayleya-Hamiltona mówi, że wstawienie dowolnej macierzy do jej wielomianu charakterystycznego daje macierz zerową. W tym  celu musimy zamienić zwykłe potęgowanie liczby (Power)  na potęgowanie macierzy (MatrixPower) oraz wyraz wolny będący liczbą na macierz jednostkową mnożona przez tę liczbę, gdyż wynik dodawania macierzy do liczby nie ma sensu. Warto zwrócić uwagę, że macierz A jest całkowicie dowolna; może być osobliwa, nilpotentna, zespolona itd.

(Simplify[poly-Det[A]] //Expand);
%/.Power->MatrixPower;
% /. λ->A;
%+Det[A]IdentityMatrix[n];
Simplify[%] //MatrixForm