% 1) REPREZENTACE MATICOVEHO OPERATORU V JINE BAZI
echo on
format compact
format short
% Pr 35: Ho/35
% Operator dany matici A v bazi U:
load Diagmat1
A
% Nova baze V: v1=2*u1+3*u2+u3; v2=3*u1+4*u2+u3; v3=u1+2*u2+2*u3
% Matice P prechodu U<-V:
P=[2 3 1;3 4 2;1 1 2]
% Matice reprezentujici operator A v bazi V:
X = inv(P)*A*P
% Neboli reseni maticove rovnice P*X=A*P:
X = P\A*P
% ***************************************
pause;
% 2) DIAGONALIZOVANI MATIC
% % Pr. 13.4/KaSk s.123 (matice nediagonalizovatelna v komplexnim oboru)
load Diagmat2
A
n = length(A)
% Charakteristicky polynom:
p = poly(A)
pause;
% Vlastni cisla jako jeho koreny:
lambda = roots(p)
% Zkontrolujeme presnost korenu
pk = polyval(p,lambda)
if max(abs(pk)) > 100*eps
   lambda = round(lambda)  % Zkusime zaokrouhlit
   pk = polyval(p,lambda)
   if max(abs(pk)) > 100*eps
      warning('Problem s presnosti korenu: osetrete rucne')
% Misto numerickeho vypoctu lze spocist symbolicky bez rizika
% nepresnosti pomoci solve a pak symbolicky vyraz vyhodnotit
% numericky pomoci eval:
% lambda = eval(solve('x^3-x^2-5*x-3'))
      keyboard
   end
end
pause;
% Algebraicke nasobnosti:
an1 = 1
an2 = 2
% Matice soustav pro hledani vlastnich vektoru prislusnych
% k jednotlivym vlastnim cislum:
A1 = lambda(1)*eye(n)-A
A2 = lambda(2)*eye(n)-A
pause
% Geometricke nasobnosti:
gn1 = n-rank(A1)
gn2 = n-rank(A2)
an1==gn1 & an2==gn2
% Neni diagonalizovatelna, nebot an2~=gn2
% ***************************************
pause;
% Pr. 13.5/KaSk s.124 (matice diagonalizovatelna dokonce v realnem oboru)
load Diagmat3
A
n = length(A)
% Charakteristicky polynom:
p = poly(A)
pause;
% Vlastni cisla jako jeho koreny:
lambda = roots(p)
% Zkontrolujeme presnost korenu
pk = polyval(p,lambda)
if max(abs(pk)) > 100*eps
   lambda = round(lambda)  % Zkusime zaokrouhlit
   pk = polyval(p,lambda)
   if max(abs(pk)) > 100*eps
      warning('Problem s presnosti korenu: osetrete rucne')
      keyboard
   end
end
pause;
% Algebraicke nasobnosti:
an1 = 1
an2 = 2
% Matice soustav pro hledani vlastnich vektoru prislusnych
% k jednotlivym vlastnim cislum:
A1 = lambda(1)*eye(n)-A
A2 = lambda(2)*eye(n)-A
pause
% Geometricke nasobnosti:
gn1 = n-rank(A1)
gn2 = n-rank(A2)
an1==gn1 & an2==gn2
% Matice je diagonalizovatelna, nebot an1==gn1 a an2=gn2
pause;
% Nalezneme baze vektorovych prostoru prislusnych vlastnim cislum:
% Baze algoritmem prevodu na schodovity tvar (nemusi byt ortonormalni)
V1 = null(A1,'r')
V1'*V1  % Kontrola (ne)ortonormality
V2 = null(A2,'r')
V2'*V2  % Kontrola (ne)ortonormality
V = [V1,V2] % Sloupce jsou vlastni vektory
V'*V        % Neni ortonormalni
D = diag(lambda)
% Kontrola A*V = V*D
A*V - V*D
% Diagonalizace
inv(V)*A*V
pause
% Rozklad se stejnymi vlasnimi cisly, ale pripadne jinymi vlastnimi
% vektory primo pomoci prikazu 'eig':
help eig
pause
[V1,D1]=eig(A)
V1'*V1     % Neni ortonormalni
% Diagonalizace
inv(V1)*A*V1
pause
% Jen vlastni cisla spocteme pomoci
eig(A)
% ******************************************************
pause;
% Pr. 13.6/KaSk s.128
% (symetricka matice ortogonalne diagonalizovatelna dokonce v realnem oboru)
load Diagmat4
A
n = length(A)
% Charakteristicky polynom:
p = poly(A)
pause;
% Vlastni cisla jako jeho koreny:
lambda = roots(p)
% Zkontrolujeme presnost korenu
pk = polyval(p,lambda)
if max(abs(pk)) > 100*eps
   lambda = round(lambda)  % Zkusime zaokrouhlit
   pk = polyval(p,lambda)
   if max(abs(pk)) > 100*eps
      warning('Problem s presnosti korenu: osetrete rucne')
      keyboard
   end
end
pause;
% Algebraicke nasobnosti:
an1 = 1
an2 = 2
% Matice soustav pro hledani vlastnich vektoru prislusnych
% k jednotlivym vlastnim cislum:
A1 = lambda(1)*eye(n)-A
A2 = lambda(2)*eye(n)-A
pause
% Geometricke nasobnosti:
gn1 = n-rank(A1)
gn2 = n-rank(A2)
an1==gn1 & an2==gn2
% Matice je diagonalizovatelna, nebot an1==gn1 a an2=gn2
pause;
% Nalezneme baze vektorovych prostoru prislusnych vlastnim cislum:
% Ortonormalni baze algoritmem singularniho rozkladu
V1 = null(A1)
V1'*V1  % Kontrola ortonormality
V2 = null(A2)
V2'*V2  % Kontrola ortonormality
V = [V1,V2] % Sloupce jsou vlastni vektory
V'*V        % Je ortonormalni
D = diag(lambda)
% Kontrola A*V = V*D
A*V - V*D
% Ortogonalni (unitarni) diagonalizace
V'*A*V
pause;
% Rozklad se stejnymi vlasnimi cisly, ale pripadne jinymi vlastnimi
% vektory primo pomoci prikazu 'eig':
[V1,D1]=eig(A)
V1'*V1
% Ortogonalni (unitarni) diagonalizace
V1'*A*V1
% ******************************************************
pause;
% Ukazka nalezeni singularniho rozkladu pomoci
% prikazu 'svd' (Singular Value Decomposition)
help svd
pause;
m=3, n=5
% Nahodne vygenerujeme celociselnou matici:
A = -2+round(10*rand(m,n))
[U,D,V] =svd(A)
pause
% Kontrola ortonormality matic U,V:
U'*U
V'*V
% Kontrola rozkladu A = U*D*V'
A - U*D*V'
pause
% Neboli zapis diagonalizace:
U'*A*V
% je prave matice D, na jejiz diagonale jsou tzv. singularni cisla.
D
pause;
% Jen singularni cisla pocita jednoduchy prikaz:
svd(A)
% Jsou to odmocniny z kladnych vlastnich cisel matice A'*A:
lambda = eig(A'*A)
sqrt(lambda(lambda>100*eps))
% ********************************
pause;
% Pomoci 'svd' lze unitarne diagonalizovat i nediagonalizovatelne
% ctvercove matice. Musime vsak pouzit ruzne matice U,V.
% Uzijeme matici z prikladu 13.4:
load Diagmat2
A
n = length(A)
[U,D,V] =svd(A)
% Matice U a V jsou ruzne
pause
% Kontrola ortonormality matic U,V:
U'*U
V'*V
% Kontrola rozkladu A = U*D*V'
A - U*D*V'
pause
% Neboli zapis diagonalizace:
U'*A*V
% je prave matice D, na jejiz diagonale jsou tzv. singularni cisla.
D