tmcv12

echo on
format compact
format short
% Pr 35: Ho/35
% Operator dany matici A v bazi U:
load Diagmat1
A
A =
     3    -2    -1
     4    -3     2
     8    -7     6
% 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]
P =
     2     3     1
     3     4     2
     1     1     2
% Matice reprezentujici operator A v bazi V:
X = inv(P)*A*P
X =
    9.0000    6.0000   16.0000
   -6.0000   -4.0000  -12.0000
   -1.0000    0.0000    1.0000
% Neboli reseni maticove rovnice P*X=A*P:
X = P\A*P
X =
    9.0000    6.0000   16.0000
   -6.0000   -4.0000  -12.0000
   -1.0000         0    1.0000
% ***************************************
pause;
% 2) DIAGONALIZOVANI MATIC
% % Pr. 13.4/KaSk s.123 (matice nediagonalizovatelna v komplexnim oboru)
load Diagmat2
A
A =
     2     1     1
     2     1    -2
    -1     0    -2
n = length(A)
n =
     3
% Charakteristicky polynom:
p = poly(A)
p =
    1.0000   -1.0000   -5.0000   -3.0000
pause;
% Vlastni cisla jako jeho koreny:
lambda = roots(p)
lambda =
    3.0000
   -1.0000
   -1.0000
% Zkontrolujeme presnost korenu
pk = polyval(p,lambda)
pk =
  1.0e-013 *
   -0.2931
   -0.0089
   -0.0089
if max(abs(pk)) > 100*eps
   lambda = round(lambda)  % Zkusime zaokrouhlit
lambda =
     3
    -1
    -1
   pk = polyval(p,lambda)
pk =
  1.0e-013 *
   -0.2931
    0.0355
    0.0355
   if max(abs(pk)) > 100*eps
      warning('Problem s presnosti korenu: osetrete rucne')
Warning: Problem s presnosti korenu: osetrete rucne
> In D:\Vesely\TEXTY\VYUKA\ESF\2004-2005\Pmzmii-2005\cviceni\tmcv12.m at line 34
% 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
lambda = eval(solve('x^3-x^2-5*x-3'))
lambda =
     3
    -1
    -1
return
   end
end
pause;
% Algebraicke nasobnosti:
an1 = 1
an1 =
     1
an2 = 2
an2 =
     2
% Matice soustav pro hledani vlastnich vektoru prislusnych
% k jednotlivym vlastnim cislum:
A1 = lambda(1)*eye(n)-A
A1 =
     1    -1    -1
    -2     2     2
     1     0     5
A2 = lambda(2)*eye(n)-A
A2 =
    -3    -1    -1
    -2    -2     2
     1     0     1
pause
% Geometricke nasobnosti:
gn1 = n-rank(A1)
gn1 =
     1
gn2 = n-rank(A2)
gn2 =
     1
an1==gn1 & an2==gn2
ans =
     0
% Neni diagonalizovatelna, nebot an2~=gn2
% ***************************************
pause;
% Pr. 13.5/KaSk s.124 (matice diagonalizovatelna dokonce v realnem oboru)
load Diagmat3
A
A =
     3     0     6
     0    -3     0
     5     0     2
n = length(A)
n =
     3
% Charakteristicky polynom:
p = poly(A)
p =
     1    -2   -39   -72
pause;
% Vlastni cisla jako jeho koreny:
lambda = roots(p)
lambda =
   8.0000          
  -3.0000 + 0.0000i
  -3.0000 - 0.0000i
% Zkontrolujeme presnost korenu
pk = polyval(p,lambda)
pk =
  1.0e-012 *
  -0.9805          
   0.0426 + 0.0000i
   0.0426 - 0.0000i
if max(abs(pk)) > 100*eps
   lambda = round(lambda)  % Zkusime zaokrouhlit
lambda =
     8
    -3
    -3
   pk = polyval(p,lambda)
pk =
     0
     0
     0
   if max(abs(pk)) > 100*eps
   end
end
pause;
% Algebraicke nasobnosti:
an1 = 1
an1 =
     1
an2 = 2
an2 =
     2
% Matice soustav pro hledani vlastnich vektoru prislusnych
% k jednotlivym vlastnim cislum:
A1 = lambda(1)*eye(n)-A
A1 =
     5     0    -6
     0    11     0
    -5     0     6
A2 = lambda(2)*eye(n)-A
A2 =
    -6     0    -6
     0     0     0
    -5     0    -5
pause
% Geometricke nasobnosti:
gn1 = n-rank(A1)
gn1 =
     1
gn2 = n-rank(A2)
gn2 =
     2
an1==gn1 & an2==gn2
ans =
     1
% 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 =
    1.2000
         0
    1.0000
V1'*V1  % Kontrola (ne)ortonormality
ans =
    2.4400
V2 = null(A2,'r')
V2 =
     0    -1
     1     0
     0     1
V2'*V2  % Kontrola (ne)ortonormality
ans =
     1     0
     0     2
V = [V1,V2] % Sloupce jsou vlastni vektory
V =
    1.2000         0   -1.0000
         0    1.0000         0
    1.0000         0    1.0000
V'*V        % Neni ortonormalni
ans =
    2.4400         0   -0.2000
         0    1.0000         0
   -0.2000         0    2.0000
D = diag(lambda)
D =
     8     0     0
     0    -3     0
     0     0    -3
% Kontrola A*V = V*D
A*V - V*D
ans =
     0     0     0
     0     0     0
     0     0     0
% Diagonalizace
inv(V)*A*V
ans =
    8.0000         0    0.0000
         0   -3.0000         0
         0         0   -3.0000
pause
% Rozklad se stejnymi vlasnimi cisly, ale pripadne jinymi vlastnimi
% vektory primo pomoci prikazu 'eig':
help eig

 EIG    Eigenvalues and eigenvectors.
    E = EIG(X) is a vector containing the eigenvalues of a square 
    matrix X.
 
    [V,D] = EIG(X) produces a diagonal matrix D of eigenvalues and a
    full matrix V whose columns are the corresponding eigenvectors so
    that X*V = V*D.
 
    [V,D] = EIG(X,'nobalance') performs the computation with balancing
    disabled, which sometimes gives more accurate results for certain
    problems with unusual scaling. If X is symmetric, EIG(X,'nobalance')
    is ignored since X is already balanced.
 
    E = EIG(A,B) is a vector containing the generalized eigenvalues
    of square matrices A and B.
 
    [V,D] = EIG(A,B) produces a diagonal matrix D of generalized
    eigenvalues and a full matrix V whose columns are the
    corresponding eigenvectors so that A*V = B*V*D.
 
    EIG(A,B,'chol') is the same as EIG(A,B) for symmetric A and symmetric
    positive definite B.  It computes the generalized eigenvalues of A and B
    using the Cholesky factorization of B.
    EIG(A,B,'qz') ignores the symmetry of A and B and uses the QZ algorithm.
    In general, the two algorithms return the same result, however using the
    QZ algorithm may be more stable for certain problems.
    The flag is ignored when A and B are not symmetric.
 
    See also CONDEIG, EIGS.

 Overloaded methods
    help lti/eig.m
    help sym/eig.m

pause
[V1,D1]=eig(A)
V1 =
    0.7682   -0.7071         0
         0         0    1.0000
    0.6402    0.7071         0
D1 =
     8     0     0
     0    -3     0
     0     0    -3
V1'*V1     % Neni ortonormalni
ans =
    1.0000   -0.0905         0
   -0.0905    1.0000         0
         0         0    1.0000
% Diagonalizace
inv(V1)*A*V1
ans =
    8.0000         0         0
    0.0000   -3.0000         0
         0         0   -3.0000
pause
% Jen vlastni cisla spocteme pomoci
eig(A)
ans =
     8
    -3
    -3
% ******************************************************
pause;
% Pr. 13.6/KaSk s.128
% (symetricka matice ortogonalne diagonalizovatelna dokonce v realnem oboru)
load Diagmat4
A
A =
     8    -2     2
    -2     5     4
     2     4     5
n = length(A)
n =
     3
% Charakteristicky polynom:
p = poly(A)
p =
     1   -18    81     0
pause;
% Vlastni cisla jako jeho koreny:
lambda = roots(p)
lambda =
         0
    9.0000
    9.0000
% Zkontrolujeme presnost korenu
pk = polyval(p,lambda)
pk =
  1.0e-012 *
         0
    0.1279
    0.1279
if max(abs(pk)) > 100*eps
   lambda = round(lambda)  % Zkusime zaokrouhlit
lambda =
     0
     9
     9
   pk = polyval(p,lambda)
pk =
     0
     0
     0
   if max(abs(pk)) > 100*eps
   end
end
pause;
% Algebraicke nasobnosti:
an1 = 1
an1 =
     1
an2 = 2
an2 =
     2
% Matice soustav pro hledani vlastnich vektoru prislusnych
% k jednotlivym vlastnim cislum:
A1 = lambda(1)*eye(n)-A
A1 =
    -8     2    -2
     2    -5    -4
    -2    -4    -5
A2 = lambda(2)*eye(n)-A
A2 =
     1     2    -2
     2     4    -4
    -2    -4     4
pause
% Geometricke nasobnosti:
gn1 = n-rank(A1)
gn1 =
     1
gn2 = n-rank(A2)
gn2 =
     2
an1==gn1 & an2==gn2
ans =
     1
% 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 =
    0.3333
    0.6667
   -0.6667
V1'*V1  % Kontrola ortonormality
ans =
    1.0000
V2 = null(A2)
V2 =
    0.9428         0
   -0.2357    0.7071
    0.2357    0.7071
V2'*V2  % Kontrola ortonormality
ans =
    1.0000    0.0000
    0.0000    1.0000
V = [V1,V2] % Sloupce jsou vlastni vektory
V =
    0.3333    0.9428         0
    0.6667   -0.2357    0.7071
   -0.6667    0.2357    0.7071
V'*V        % Je ortonormalni
ans =
    1.0000    0.0000   -0.0000
    0.0000    1.0000    0.0000
   -0.0000    0.0000    1.0000
D = diag(lambda)
D =
     0     0     0
     0     9     0
     0     0     9
% Kontrola A*V = V*D
A*V - V*D
ans =
  1.0e-015 *
         0         0    0.2220
         0   -0.4441         0
         0    0.4441   -0.8882
% Ortogonalni (unitarni) diagonalizace
V'*A*V
ans =
         0         0         0
         0    9.0000    0.0000
         0    0.0000    9.0000
pause;
% Rozklad se stejnymi vlasnimi cisly, ale pripadne jinymi vlastnimi
% vektory primo pomoci prikazu 'eig':
[V1,D1]=eig(A)
V1 =
    0.3333    0.2981    0.8944
    0.6667    0.5963   -0.4472
   -0.6667    0.7454         0
D1 =
     0     0     0
     0     9     0
     0     0     9
V1'*V1
ans =
    1.0000   -0.0000    0.0000
   -0.0000    1.0000         0
    0.0000         0    1.0000
% Ortogonalni (unitarni) diagonalizace
V1'*A*V1
ans =
    0.0000   -0.0000    0.0000
   -0.0000    9.0000         0
    0.0000         0    9.0000
% ******************************************************
pause;
% Ukazka nalezeni singularniho rozkladu pomoci
% prikazu 'svd' (Singular Value Decomposition)
help svd

 SVD    Singular value decomposition.
    [U,S,V] = SVD(X) produces a diagonal matrix S, of the same 
    dimension as X and with nonnegative diagonal elements in
    decreasing order, and unitary matrices U and V so that
    X = U*S*V'.
 
    S = SVD(X) returns a vector containing the singular values.
 
    [U,S,V] = SVD(X,0) produces the "economy size"
    decomposition. If X is m-by-n with m > n, then only the
    first n columns of U are computed and S is n-by-n.
 
    See also SVDS, GSVD.

 Overloaded methods
    help sym/svd.m

pause;
m=3, n=5
m =
     3
n =
     5
% Nahodne vygenerujeme celociselnou matici:
A = -2+round(10*rand(m,n))
A =
     2     2     2    -1     4
     7     7     6     0     1
     7    -1    -2     0     0
[U,D,V] =svd(A)
U =
   -0.3211    0.1698    0.9317
   -0.9046    0.2363   -0.3549
   -0.2804   -0.9567    0.0777
D =
   12.6460         0         0         0         0
         0    6.7208         0         0         0
         0         0    3.5929         0         0
V =
   -0.7067   -0.6998   -0.0214    0.0198   -0.1000
   -0.5293    0.4390   -0.1944   -0.1985    0.6707
   -0.4356    0.5462   -0.1172    0.1687   -0.6853
    0.0254   -0.0253   -0.2593    0.9353    0.2382
   -0.1731    0.1362    0.9385    0.2388    0.1169
pause
% Kontrola ortonormality matic U,V:
U'*U
ans =
    1.0000   -0.0000    0.0000
   -0.0000    1.0000    0.0000
    0.0000    0.0000    1.0000
V'*V
ans =
    1.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    1.0000         0   -0.0000    0.0000
    0.0000         0    1.0000    0.0000   -0.0000
    0.0000   -0.0000    0.0000    1.0000   -0.0000
    0.0000    0.0000   -0.0000   -0.0000    1.0000
% Kontrola rozkladu A = U*D*V'
A - U*D*V'
ans =
  1.0e-013 *
    0.0133         0         0   -0.0089    0.0355
    0.0089   -0.0089   -0.0178   -0.0117    0.0555
    0.0266    0.0377    0.0089    0.0737   -0.2742
pause
% Neboli zapis diagonalizace:
U'*A*V
ans =
   12.6460    0.0000    0.0000    0.0000   -0.0000
   -0.0000    6.7208    0.0000    0.0000    0.0000
    0.0000   -0.0000    3.5929         0   -0.0000
% je prave matice D, na jejiz diagonale jsou tzv. singularni cisla.
D
D =
   12.6460         0         0         0         0
         0    6.7208         0         0         0
         0         0    3.5929         0         0
pause;
% Jen singularni cisla pocita jednoduchy prikaz:
svd(A)
ans =
   12.6460
    6.7208
    3.5929
% Jsou to odmocniny z kladnych vlastnich cisel matice A'*A:
lambda = eig(A'*A)
lambda =
   -0.0000
    0.0000
   12.9086
   45.1692
  159.9222
sqrt(lambda(lambda>100*eps))
ans =
    3.5929
    6.7208
   12.6460
% ********************************
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
A =
     2     1     1
     2     1    -2
    -1     0    -2
n = length(A)
n =
     3
[U,D,V] =svd(A)
U =
   -0.5937    0.4851    0.6421
   -0.7936   -0.4851   -0.3673
    0.1333   -0.7276    0.6729
D =
    3.3028         0         0
         0    3.0000         0
         0         0    0.3028
V =
   -0.8804    0.2425   -0.4075
   -0.4200   -0.0000    0.9075
    0.2201    0.9701    0.1019
% Matice U a V jsou ruzne
pause
% Kontrola ortonormality matic U,V:
U'*U
ans =
    1.0000   -0.0000    0.0000
   -0.0000    1.0000    0.0000
    0.0000    0.0000    1.0000
V'*V
ans =
    1.0000   -0.0000   -0.0000
   -0.0000    1.0000   -0.0000
   -0.0000   -0.0000    1.0000
% Kontrola rozkladu A = U*D*V'
A - U*D*V'
ans =
  1.0e-014 *
         0         0    0.0444
    0.0222    0.0333   -0.1776
         0    0.0111   -0.0888
pause
% Neboli zapis diagonalizace:
U'*A*V
ans =
    3.3028   -0.0000   -0.0000
   -0.0000    3.0000    0.0000
    0.0000    0.0000    0.3028
% je prave matice D, na jejiz diagonale jsou tzv. singularni cisla.
D
D =
    3.3028         0         0
         0    3.0000         0
         0         0    0.3028
diary off