Macierze i wektory

Podstawowe operacje na wektorach

Wprowadzanie wektorów i operacje arytmetyczne

Najprostszym i naturalnym sposobem wprowadzenia wektora jest lista (List), na ogół wprowadzana z użyciem nawiasów wąsatych:

u={1,-3,4}

v={2,3,-1}

u=List[1,-3,4]

u//FullForm

Mnożenie wektora przez liczbę:

2*u

Dodawanie wektorów:

u+v

Odejmowanie:

u-v

Kombinacja liniowa:

2u-3v

Składowe wektorów mogą być symbolami:

r={x,y,z}

r+r

3r

-r

r-r

Większość funkcji programu Mathematica ma atrybut Listable, co oznacza, że zastosowane do listy (wektora), podziałają na każdy element (składową) osobno! Jest to źródłem bardzo częstych nieporozumień:

Attributes[Sin]

Sin[u]

u

u^2 (* Każda składowa została podniesiona do kwadratu ! *)

x^2//FullForm

Attributes[Power]

u u

u*v  (* To NIE JEST iloczyn wektorów, składowe zostały pomnożone przez siebie *)

u.v (* Tak liczymy iloczyn skalarny *)

u.u (* A tak kwadrat długości wektora ! *)

Atrybut Listable ma charakter ,,przeciążenia’’ operacji matematycznych takich jak  np: Plus, co oznacza, że natrafiając na obiekt inny niż liczba lub pojedynczy symbol, program stara się wykonać sensowną operację, analogiczną do tej dla liczb. Czasem jest to oczywiste, jak przy dodawaniu (dodajemy składowe osobno), czasem zaskakuje, jak przy mnożeniu.

Dlatego Mathematica wykonuje poniższą operację, pomimo że nie ma ona na ogół sensu matematycznego (nie można dodawać liczby do wektora):

1+u (* istnieje sensowne ,,przeciążenie'' dodawania w tym wypadku *)

Dopiero próba dodania wektorów z przestrzenie o różnej liczbie wymiarów produkuje błąd:

{1,1}+u

{1}+u

Iloczyn skalarny

u

v

Iloczyn skalarny, a także inne operacje analogiczne do mnożenia na złożonych obiektach jak macierze i tensory, wykonujemy operatorem Dot, skrótowo zapisywanym jako         .         (kropka dziesiętna!).

Dot[u,v]

u.v

v.u

Total[u*v]

u.u

Clear[a]

a={ax,ay,az}

b={bx,by,bz}

Szkolny wzór na iloczyn skalarny:

a.b

Iloczyn skalarny jest przemienny:

a.b-b.a

a={a1,a2,a3,a4}

b={b1,b2,b3,b4}

a.b

Kwadrat długości wektora:

a.a

Kropka ma dwojakie znaczenie: iloczyn skalarny lub ,,przecinek’’ w ułamku dziesiętnym:

2.a (* Tu kropka jest interpretowana jako część liczby zmiennoprzecinkowej 2.0 *)

2.u

FullForm[%]

λ.u (* λ jest pewnym niesprecyzowanym symbolem; jeżeli stanie się wektorem o 3 składowych to wyrażenie będzie miało sens, na razie nic się nie dzieje *)

Iloczyn skalarny ma sens tylko dla wektorów o tej samej liczbie wymiarów:

{2,3}.{1,2,3}

Długości i kąty pomiędzy wektorami

Kąt pomiędzy wektorami można obliczyć następująco:

u

v

VectorAngle[u,v]

Wartość numeryczna w radianach:

N[%]

... i w stopniach:

%/Degree

Długość wektora, a także ,,normę’’ innych obiektów (macierze itp.), obliczamy poleceniem Norm:

Norm[u]

Sqrt[u.u]

Szkolny wzór na kosinus kąta pomiędzy wektorami:

u.v/Norm[u]/Norm[v]

Ogólne symbole literkowe, jak zawsze, traktowane są jako liczby zespolone:

a

b

VectorAngle[a,b]

ComplexExpand[%,TargetFunctions->{Re,Im}]

Powyższy przykład demonstruje konieczność rozumienia co kryje się za poleceniem VectorAngle (ang: kąt między wektorami ):

cosAlpha=a.b/Sqrt[a.a]/Sqrt[b.b]

Norm[a]

ComplexExpand[%] (* frustrujące... *)

Sqrt[a.a]  (* upragniony wzór na długość wektora ! *)

Podsumowanie: operator Dot jest niesłychanie użyteczny, a jego stosowanie uczy pisania wzorów w formie zwektoryzowanej, co jest czytelne i w przyszłości zaprocentuje szybciej działającymi programami, jeżeli tylko software będzie umiał wykorzystać wektorowe możliwości procesorów (SSE, AVX, karty graficzne). Funkcje typu Norm ( długość wektora)  oraz VectorAngle (kąt pomiędzy wektorami ) są w praktyce uciążliwe w stosowaniu do wzorów zawierających symbole literkowe, a powodu założenia, że reprezentują liczby zespolone (przydają się w obliczeniach na konkretnych liczbach).

Iloczyn wektorowy

Iloczyn wektorowy obliczamy poleceniem Cross:

Cross[u,v]

Iloczyn wektorowy wektorów równoległych to zero:

Cross[2u,-u]

Iloczyn wektorowy jest prostopadły do każdego ze składników:

Cross[u,v].u

Cross[u,v].v

Iloczyn wektorowy jest wektorem, więc można go składać wielokrotnie:

Cross[Cross[u,v],v]

Nest[Cross[#,v]&,u,5]

NestList[Cross[#,u]&,v,9]

NestList[Cross[#,{1,0,0}]&,{0,1,0},9]

Iloczyn wektorowy w klasycznej szkolnej postaci, ma sens tylko w 3 wymiarach:

Cross[{1,2},{3,-1}]

W ogólności, iloczyn wektorowy daje wektor prostopadły do podprzestrzeni n-1 wymiarowej, wyznaczonej przez n-1 wektorów:

Cross[{1,2}] (* wektor prostopadły do prostej równoległej do wektora {1,2} *)

Cross[{1,0,0},{0,1,0}] (* wektor prostopadły do 2 wektorów w przestrzeni 3 wymiarowej *)

Tak obliczylibyśmy wektor prostopadły do każdego z trzech wektorów w przestrzeni 4-wymiarowej:

Cross[{2,3,4,5},{-1,2,3,0}, {1,1,1,1}]

{2,3,4,5}.%

Przykład: udowodnij wzór a x (b x c) = b a.c - c a.b

Wprowadzamy 3 dowolne wektory:

a={a1,a2,a3}

b={b1,b2,b3}

c={c1,c2,c3}

Lewa strona tożsamości:

lewa=Cross[a,Cross[b,c]]

Prawa strona:

prawa=b a.c - c a.b

lewa==prawa

Simplify[%]

Lub:

lewa-prawa

Simplify[%]

Operowanie macierzami kwadratowymi

Wprowadzanie macierzy

Macierz jest po prostu listą list. Każda ,,mała’’ lista stanowi jeden wiersz macierzy:

m={{1,2},{3,4}}

m={
{1,2}, (* pierwszy wiersz *)
{3,4}(* drugi wiersz *)
}

Polecenia MatrixForm i pokrewne służą tylko i wyłącznie do prezentacji graficznej! Nie są potrzebne do wykonywania operacji matematycznych; co gorsza ich niewłaściwe zastosowanie uniemożliwia rachunki.

MatrixForm[m]

TableForm[%]

m.m

m=MatrixForm[{{1,2},{3,4}}] (* ŹLE! Polecenie MatrixForm nie słuzy do tworzenia macierzy, tylko ładnego wyświetlania ich zawartości na ekranie! *)

m.m (* Na skutek błędu powyżej mnożenie jest niemożliwe *)

FullForm[%]

Sposoby tworzenia macierzy o elementach będących symbolami literowymi:

a={{a11,a12},{a21,a22}} (* ręczne wpisania, optymalne dla małych macierzy *)

Array[q,{3,3}] (* Użycie Array *)

%[[1,1]]

FullForm[%]

Poniżej preferowany przeze mnie sposób na tworzenie dużych macierzy, polegający na składaniu łańcuchów tekstowych:

Table[ToExpression["a"<>ToString[i]<>ToString[j]],{i,1,10},{j,1,10}]

MatrixForm[%]

Wyznaczniki

Wyznacznik obliczamu poleceniem Det, podobnie jak to czasem zapisuje się na tablicy:

m={{1,2},{3,4}}

Det[m]

Macierz 5x5 zawierająca losowe liczby całkowite z przedziału 0 do 4 :

RandomInteger[4,{5,5}]

MatrixForm[%]

Det[%%]

Clear[a,b,c,d]

Znany wzór na wyznacznik macierzy 2x2:

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

Mniej znany, ale jednak znany, wzór:

Det[{{a,b,c},{d,e,f},{g,h,i}}]

Mathematica radzi sobie także z większymi macierzami, ale w przypadku macierzy wypełnionych symbolami nie wykonamy obliczeń ze znacznie większymi od tej poniżej. U mnie pamięci brakło już przy macierzy 12x12.

Table[ToExpression["a"<>ToString[i]<>ToString[j]],{i,1,6},{j,1,6}]

MatrixForm[%]

MemoryConstrained[Det[%%],1024^2]

?MemoryConstrained

Jeżeli macierz zawiera większość zer, to prawdopodobieństwo, że wyznacznik jest równy zero jest duże:

RandomInteger[1,{8,8}]

Det[%]

MaxMemoryUsed[]

Mnożenie i dodawanie macierzy

Wprowadźmy trzy macierze A, B i F (symbole C, D, E są zarezerwowane: C - numeracja stałych całkowania, D - pochodna/różniczkowanie, E - liczbe e):

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

B={{1,1},{1,0}}

F={{1,-1},{1,-1}}

Dodawanie i odejmowanie macierzy:

A+B

A-B

Mnożenie przez liczbę:

2A

Mnożenie macierzowe:

A.B

UWAGA: poniższe operacje NIE SĄ mnożeniem macierzy !

A B (* ŹLE! Mnożenie elementów macierzy przez siebie *)

A*B (*ŹLE! Jak wyżej *)

A.B (* DOBRZE! *)

Mnożenie macierzy nie jest przemienne:

A.B==B.A

Ale jest łączne, co pozwala na sensowne wykonywanie poniższej operacji bez nawiasów:

A.B.F

(A.B).F==A.(B.F)

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

B={{b11,b12},{b21,b22}}

F={{f11,f12},{f21,f22}}

A.B.F//MatrixForm

(A.B).F-A.(B.F)

Simplify[%]

A.B-B.A

Simplify[%]

PossibleZeroQ[Det[%]]

Jeszcze raz o różnicy pomiędzy mnożeniem macierzy (Dot) a przeciążonym operatorem mnożenia (Times):

A.B//MatrixForm (*  Mnożenie macierzy *)

A*B //MatrixForm (* Mnożenie ELEMENTÓW MACIERZY ! *)

Dot[A,B]//MatrixForm (*  Mnożenie macierzy, instrukcja Dot *)

Times[A,B]//MatrixForm (*  Mnożenie liczb, instrukcja Times (ang. razy ) *)

A+B//MatrixForm (* W przypadku dodawania jest jeden operator: Plus *)

Poniżej widać wyraźnie różnicę pomiedzy podniesieniem macierzy do kwadratu, a podniesieniem każdego elementu macierzy do kwadratu:

A^2//MatrixForm

A*A//MatrixForm

MatrixPower[A,2]//MatrixForm

Power[A,2]//MatrixForm

A.A//MatrixForm

Macierz jednostkowa i zerowa

Macierz jednostkowa, zawiera same jedynki na przekątnej. Generuje ją polecenie IdentityMatrix:

IdentityMatrix[4]

MatrixForm[%]

Table[If[i==j,1,0],{i,1,3},{j,1,3}] (* Inny sposób *)

Array[KroneckerDelta,{3,3}]//MatrixForm

f=#1^2+#2^2&

f[x,y]

Sort[{12,1,243,34,345,35,53},#1<#2&]

Array[If[#1==#2,1,0]&,{4,4}]

Array[KroneckerDelta,{3,3,3}]

id=IdentityMatrix[2]

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

Macierz jednostkowa nie zmienia mnożonej przez nią macierzy a:

id.a

id*a (* ŹLE !!! To nie jest mnożenie macierzowe *)

a.id

Macierz zerowa jest wypełniona zamymi zerami:

zero=id-id

zero.a

a+zero

Inne sposoby utworzenia macierzy zerowej:

Array[0&,{3,3}]

Table[0,{i,1,3},{j,1,3}]

Często w praktyce mamy do czynienia z macierzami prawie wszędzie wypełnionymi zerami. Takie macierze nazywamy rzadkimi, i istnieje specjalny mechanizm zapisywania ich w pamięci (SparseArray), tak aby nie marnować jej na zera.

SparseArray[{},{3,3}]

MatrixForm[%]

Np: niech tylko środkowy element (2 rząd, 2 wiersz)  będzie niezerowy:

mS=SparseArray[{{2,2}->1},{3,3}]

MatrixForm[mS]

mN=Normal[mS] (* Konwersja na ,,normalną'' macierz *)

MatrixForm[mN]

Różnicę w ilości zajętej pamięci widać dopiero przy macierzach o rozmiarze 100x100, teraz pogorszyliśmy sytuację:

ByteCount[mS]

ByteCount[mN]

Macierz odwrotna

Kluczowym pojęciem w zagadnieniach teoretycznych, a także w obliczeniach, jest pojęcie macierzy odwrotnej  (ang. inverse matrix ). Obliczamy ją poleceniem Inverse.

m={{1,2},{3,4}}

MatrixForm[%]

mINV=Inverse[m]

MatrixForm[%]

Macierz pomnożona przez odwrotną do niej naje macierz jednostkową:

m.mINV//MatrixForm

mINV.m//MatrixForm

Nie każda macierz może być odwrócona:

m={{1,1},{1,1}}

Inverse[m]

Warunkiem koniecznym jest niezerowy wyznacznik:

Det[m]

Możemy oczywiście operować na macierzach zawierających symbole literowe:

Clear[a,b,c,d]

mac={{a,b},{c,d}};
MatrixForm[mac]

Inverse[mac]//MatrixForm

Elementy macierzy odwrotnej są minorami (mniejszymi wyznacznikami) dzielonymi przez wyznacznik główny:

Det[mac]Inverse[mac]//MatrixForm

Macierzy odwrotnej można użyć do rozwiązywania liniowych układów równań, aczkolwiek zwykle szybsze jest użycie Solve lub LinearSolve. Zapis z użyciem macierzy odwrotnej jest natomiast bardziej czytelny i elegancki.

Benchmarking

Obliczanie macierzy odwrotnej i pokrewne rachunki są standardowym testem szybkości współczesnych komputerów. W programie Mathematica mamy dodatkowo możliwość porównania rachunku symbolicznego (np: na liczbach wymiernych) i numerycznego (korzystającego z koprocesora numerycznego CPU).

timing[n_]:=timing[n]=AbsoluteTiming[Inverse[RandomInteger[{-4,4},{n,n}]]][[1]]

dat=Table[{n,timing[n]},{n,2,64,1}]

dat2=Join[dat,{{128,timing[128]}}];

ListPlot[dat2,PlotRange->All]

ListLogLogPlot[dat2,Joined->True,InterpolationOrder->0]

FindFit[dat2,a X^p,{a,{p,3}},X]

MaxMemoryUsed[]

timingNUM[n_]:=timingNUM[n]=AbsoluteTiming[Inverse[N@RandomInteger[{-4,4},{n,n}]]][[1]]

datNUM=Table[{n,timingNUM[n]},{n,2,256,1}];

ListPlot[datNUM,PlotRange->All,Joined->True,InterpolationOrder->0]

ListLogLogPlot[datNUM,PlotRange->All,Joined->True,InterpolationOrder->0]

FindFit[datNUM,a X^p,{a,{p,3}},X]

Wyniki można porównać z teoretyczną wydajnością: https://en.wikipedia.org/wiki/Computational_complexity_of_mathematical_operations#Matrix_algebra

Macierze o wymiarze 1 x n, n x 1 a wektory

W matematyce ,,szkolnej’’ rozróżnianie pomiędzy wektorem i macierzą 1 na N, czy też liczbą a macierzą 1x1 nie jest zwykle potrzebne, gdyż obiekty te są izomorficzne, a odpowiednią konwersję dokonuje się domyślnie. To może prowadzić do nieporozumień podczas pracy z programami komputerowymi. Rozważmy 3 izomorficzne obiekty, trójki liczb (a konkretnie symboli interpretowanych jako liczby):

a={a1,a2,a3} (* ,,wektor''*)

b={{a1,a2,a3}} (* macierz 1x3 ,,wierszowa'' *)

c={{a1},{a2},{a3}} (* macierz 3x1 ,,kolumnowa'' *)

Graficznie Mathematica przedstawia te obiekty nastepująco:

MatrixForm[a]

MatrixForm[b]

MatrixForm[c]

Iloczyn skalarny wektorów obliczamy bez problemu:

a.a

Ale wektory ,,kolumnowe’’ i ,,wierszowe’’ nie są wektorami, ale macierzami. Obowiązuję więc reguły mnożenia macierzy: liczba kolumn pierwszego składnika musi być równa liczbie wierszy drugiego, co w poniższych przykładach nie zachodzi:

b.b

c.c

a.b

Mnożenie wektora przez macierz jest dopuszczalne, z odpowiedniej strony:

a.c

b.a

c.a

Jeszcze inne wyniki uzyskamy mnożąć macierze 1xN i Nx1 przez siebie. Obie operacje są dopuszczalne z punktu widzenia operatora Dot. W pierwszym przypadku wynikiem jest macierz 1x1:

b.c

Nie jest więc prawdziwe wyrażenie:

(b.c)[[1,1]] == a.a

%/.{a1->1,a2->2,a3->3}

Simplify[%]

W przypadku mnożenia w odwrotnej kolejności, wynikiem jest macierz 3x3, tzw. iloczyn diadyczny (https://pl.wikipedia.org/wiki/Iloczyn_diadyczny).

b.c

c.b//MatrixForm

b[[1]]

Table[Sum[Signature[{i,j,k}]a[[i]] b[[1,j]],{i,1,3},{j,1,3}],{k,1,3}]

Transpozycja macierzy

Transpozycja macierzy polega na zamianie kolumn i wierszy:

Clear[a,b,c,e,f,g,h,i]

m={{a,b,c},{d,e,f},{g,h,i}}

MatrixForm[m]

MatrixForm[Transpose[m]]

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

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

Wybieranie bloków macierzy

m//MatrixForm

m[[2,2]]

m[[-1,-1]]

m[[1;;2]]//MatrixForm

m[[2;;,2;;]]//MatrixForm

m[[1;;3,2]]

?Part

Diagonalizacja macierzy i zaawansowane operacje

Wartości i wektory własne na przykładzie macierzy (({{0, -1, 0}, {-1, -2, 0}, {0, 0, 2}}))

Sztandarowym zadaniem na ćwiczeniach i egzaminie z algebry jest znalezienie wartości własnych i wektorów własnych. Polecenia z powodów historycznych mają nazwy pochodzące z języka niemieckiego: Eigenvalues (wartości własne) Eigenvectors (wektory własne).

(* RandomInteger[{-1,1},{3,3}]*)

m={{0,-1,0},{-1,-2,0},{0,0,2}};
MatrixForm[m]

Det[m] (* Wyznacznik *)

ww=Eigenvalues[m] (* Wartości własne *)

vw=Eigenvectors[m] (* Wektory własne *)

Zgodnie z definicją, macierz m pomnożona przez wektor własny, powinna być równa wartości własnej pomnożnej przez wektor własny (m.v==λ v):

v=vw[[1]] (* pierwszy wektor własny *)

λ=ww[[1]] (* pierwsza wartość własna *)

m.v

λ v

Simplify[%]

m.v==λ v

Podobnie mozna sprawdzić pozostałe wektory i wartości własne. Zrobimy to automatycznie:

Table[m.vw[[i]]==ww[[i]] vw[[i]],{i,1,3}]

Wartości własne nie są posortowane:

ww//N

Wartości i wektory własne możemy obliczyć jednym poleceniem Eigensystem:

Eigensystem[m]

Proces diagonalizacji polega na mnożeniu macierzy przez macierz utworzoną z wektorów własnych oraz macierz do niej odwrotną.
Poniżej widzimy, że iloczyn 3 macierzy:
1) odwrotna do macierzy utworzonej z wektorów własnych Inverse[vw]
2) macierz, która na przekątnej ma wartości własne: DiagonalMatrix[ww]
3) macierz utworzona z wektorów własnych vw

daje wyjściową macierz m:

Inverse[vw].DiagonalMatrix[ww].vw//Simplify//MatrixForm

m//MatrixForm

Inverse[vw].DiagonalMatrix[ww].vw==m//Simplify

Analogicznie, poniższe mnożenie daje w wyniku macierz zdiagonalizowaną, czyli posiadającą na przekątnej wartości własne:

vw.m.Inverse[vw]//Simplify//MatrixForm

Wielomian charakterystyczny, ślad macierzy

W szkolnej algebrze, najpierw wyznaczamy wielomian charakterystyczny, a nastepnie znajdujemy jego pierwiastki, czyli wartości własne. Z definicji wielomian ten wynosi:

Clear[λ]

Det[m-λ IdentityMatrix[3]]

Solve[%==0,λ]

W. charakterystyczny można utworzyć poleceniem CharacteristicPolynomial:

poly=CharacteristicPolynomial[m,λ]

Współczynnik przy najwyższej potędze to 1, przy kolejnej jest równy śladowi macierzy (Tr), a ostatni to wyznacznik:

Det[m]

Tr[m]

m//MatrixForm

CoefficientList[poly,λ]

UWAGA: polecenie Trace NIE MA NIC wspólnego ze śladem macierzy! Służy do śledzenia obliczeń i jest narzędziem programisty.

Trace[m]

?Trace

Funkcje argumentu macierzowego: potęgowanie . Przykład:oblicz  n-tą potęgę macierzy(({{1, 1}, {1, 0}}))

Funkcji znane dla liczb rzeczywistych i zespolonych mogą być zastosowane do macierzy. Najprostszy sposób to użycie szeregu potęgowego, w którym potęgowanie zastepujemy potęgowaniem macierzy.

Rozważmy przykład:

m={{1,1},{1,0}};MatrixForm[m]

m.m //MatrixForm (* m ,,do kwadratu'' *)

m.m.m //MatrixForm (* m ,,do sześcianu'' *)

m.m.m.m//MatrixForm

m.m.m.m.m//MatrixForm

Zwracam uwagę, że poniższe wyrażenie NIE JEST potęgą macierzy:

m^5

Należy użyć polecenia MatrixPower. Ogólnie, funkcja która ma działać na macierz jako całość, a nie na jej elementy osobno, ma przedrostek Matrix...

MatrixPower[m,5]

Manipulate[MatrixPower[m,n]//MatrixForm,{n,0,100,1}]

MatrixPower[m,n]

FullSimplify[%]

Table[MatrixPower[m,n][[1,1]],{n,0,10,1}]

FindSequenceFunction[%]

Table[MatrixPower[m,n]-{{Fibonacci[n+1],Fibonacci[n]},{Fibonacci[n],Fibonacci[n-1]}},{n,0,12}]

%/.n->12

N[%]

FullSimplify[ MatrixPower[m,n]-{{Fibonacci[n+1],Fibonacci[n]},{Fibonacci[n],Fibonacci[n-1]}},n>0&&n∈Integers]

Przykład: oblicz eksponentę macierzy (({{0, 1}, {-(ω^2), 0}})) (rozwiązanie r. oscylatora harmonicznego)

Clear[k,m]

Równanie oscylatora harmonicznego można zapisać w postaci macierzowej {x’[t],v’[t]}==a.{x[t],v[t]}, gdzie a jest macierzą:

a={{0,1},{-ω^2,0}};MatrixForm[a]

Rozwiązanie to po prostu (e^(a t)):

m=MatrixExp[{{0,1},{-ω^2,0}}*t]

Rozwiązanie zagadnienia początkowego daje mnożenie macierzy m przez wektor zawierający położenie i prędkość początkową:

m.{x0,v0}

Clear[m,v]

Na tablicy zwykle nie stosuje się takich metod rozwiązania, gdyż samo obliczenie eksponenty macierzy jest bardzo czasochłonne. W sytuacji gdy jedyne co musimy zrobić, to wpisać MatrixExp, metoda okazuje się, krótka, elegancka i wygodna.

Przykład: twierdzenie Cayleya-Hamiltona

Twierdzenie mówi, że jeżeli wstawimy macierz m do jej wielomianu charakterystycznego, traktując mnożenie (potęgowanie) jako mnożenie macierzowe, to otrzymamy macierz zerową. Dowód jest prosty, gdy macierz jest diagonalizowalna. Twierdzenie jest jednak prawdziwe
dla każdej macierzy kwadratowej. Poniżej sprawdzamy to bezpośrednim rachunkiem, na moim komputerze aż do macierzy 6x6:

num=2

m=Table["m"<>ToString[i]<>ToString[j],{i,1,num},{j,1,num}]

poly=CharacteristicPolynomial[m,λ]

a=CoefficientList[poly,λ]

Sum[a[[i]]MatrixPower[m,i],{i,1,num+1}]

Simplify[%]//AbsoluteTiming

MaxMemoryUsed[]

Przykład zastosowania eksponenty macierzy: obrót dookoła dowolnej osi a zadany kąt

RotationMatrix[α,{1,1,1}]

n={1,1,1}/Norm[{1,1,1}]

MatrixExp[-α Table[ Sum[Signature[{i,j,k}] n[[k]],{k,1,3}],{i,1,3},{j,1,3}]]

%-%%

FullSimplify[%]

Przykład numeryczny: półkole Wignera

n=2^8

RandomReal[1,{n,n}];

m=1/2(%+Transpose[%]); (* trik symetryzujący macierz *)

ww=Eigenvalues[m]//Sort

dat=BinCounts[Most[ww]]

w1=ListPlot[dat]

Clear[a,r]

FindFit[dat,a Re[Sqrt[r^2-(x-r)^2]],{a,{r,Sqrt[n]/2}},x]

{a,r}={a,r}/.%

Sqrt[n]/2

w2=Plot[a Sqrt[r^2-(x-r)^2],{x,0,2r}]

Show[w1,w2]