%------------- forgatási transzformáció bemutatása 6. fólia ----
A = [2 6 4; ...
     1 4 7; ...
     1 1 1]
% A dot2dotc eljárás színezett kitöltéssel rajzolja meg a poligont
dot2dotc(A, 'b')
fi = 60;
B = forgat(A, fi)
dot2dotc(B, 'g') 
grid on, axis equal, hold off
%------------------------------------------------------------
% Projektor-mátrix 5. fólia (síkra vetítés: E - n'*n )
% Pithagoraszi téglák: [1 2 2], [2 3 6], [1 4 8], [4 4 7]
n = [1 2 2] 
n=n/sqrt(n*n') % egységnyi hosszúvá tesszük
E=eye(3), P=E-n'*n
P*P
det(P)  % projektormátrixok determinánsa 0
x = [2 -4 5]', x1=P*x, x2=P*x1

%-------- 5. fólia transzformációk szorzata -----------------
% inverz mátrix meghatározása lin. transzformációk szorzatával:

A = [5 2 3; 10 15 20; 40 50 60]
% 1. oszlopot egységvektorrá tevő transzformáció, T1 oszlopaiba
% a bázis-egységvektorok transzformáltjainak koordinátái kerülnek:

T1 = [1/5 0 0; -10/5 1 0; -40/5 0 1]
A1 = T1*A

% 2. és 3. sorok cseréje a pivotálás miatt:
T2 = [1 0 0; 0 0 1; 0 1 0]  % swap mátrix
A2 = T2*A1

% 2. oszlopot egységvektorrá tevő transzformáció:
T3 = [1 -0.4/34 0; 0 1/34 0; 0 -11/34 1]
A3 = T3*A2
rats(A3)  % milyen törtek a mátrix elemei?

% 3. oszlopot egységvektorrá tevő transzformáció:
T4 = [1 0 -3/40; 0 1 -18/40; 0 0 17/40]
A4 = T4*A3
rats(A4)

% Mivel T4*T3*T2*T1*A = E, ezért T4*T3*T2*T1 = inv(A)
invA = T4*T3*T2*T1 
rats(invA), invA*A
%--------------------------------------------------------
% Oldjuk meg Gauss eliminációval az A*x=b egyenletrendszert!

b = [10; 19; 4]

% Bővítsük ki az A mátrixot a b oszlopvektorral!

Ab = [A b]
% A Gauss elimináció végrehajtása:
% Az első oszlopot a főátlóbeli elem kivételével nullázuk és
% egységvektorrá alakítjuk, így a 2. és 3. egyenletből
% az x1 ismeretlent kiküszöböljük következőképp:
% a 2. sorból kivonjuk az 1. sor felét
% a 3. sorból kivonjuk az 1. sor felét
% az 1. sort elosztjuk 2-vel
% Ezeket a lépéseket végrehajtó T1 transzformáció mátrixa:
T1=[1/2 0 0; -0.5 1 0; -0.5 0 1]
Ab1=T1*Ab

% A második oszlopban a legnagyobb abszolútértékű elem a 3.
% sorban van, így felcseréljük a 2. és 3. sort a P1 mátrix-szal
% val való balról szorzás segítségével és az így kapott második
% sort rögtön elosztjuk -2-vel, ezzel a második oszlop főátlóbeli
% eleme 1 lesz:
T2 = [1 0 0; 0 0 -0.5; 0 1 0]
Ab2=T2*Ab1
% Az második oszlopot a főátlóbeli elem kivételével nullázuk és 
% egységvektorrá alakítjuk, így az 1. és 3. egyenletből
% az x2 ismeretlent kiküszöböljük következőképp:
% az 1. sorból kivonjuk a 2. sor 3-szorosát
% a 3. sorból kivonjuk a 2. sort 
% a 2. sort nem bántjuk, mert a főátlóban 1 van
% A megfelelő T2 mátrix a fenti transzformáció mátrixát az 
% e1,e2,e3 egységvektorkon végrehajtott transzformációval kapjuk:
T3=[1 -3 0; 0 1 0; 0 -1 1]
Ab3=T3*Ab2

% A harmadik oszlopot a főátlóbeli elem kivételével nullázuk és 
% egységvektorrá alakítjuk, így az 1. és 2. egyenletből
% az x3 ismeretlent kiküszöböljük következőképp:
% az 1. sorból kivonjuk a 3. sor (0.5/4.5)-szeresét
% a 2. sorból kivonjuk a 3. sor (0.5/4.5)-szeresét
% a 3. sort elosztjuk 4.5-el
% Ezzel az A együthatómátrix transzformáltja egységmátrix lett
T4= [1 0 -0.5/4.5; 0 1 -0.5/4.5; 0 0 1/4.5]
Ab4=T4*Ab3
x = Ab4(:,4)

% A négy transzformációs mátrix szorzata az A mátrix inverze lesz
invA = T4*T3*T2*T1

% ellenőrzések:
inv(A)
A*x, b

%----------------- 9-10. fólia ----------------------------------
% skaláris szorzat, vektoriális szorzat:

a = [3,4,0], b = [-4,3,0], dot(a, b)  % merőlegesek
a*b'

a = [1 0 0], b = [0 1 0], cross(a, b)

a = rand(1, 3), b = rand(1, 3), c = rand(1, 3), 
dot(cross(a, b), c), det([a; b; c])

subspace(a', b'), subspace(a, b)
acos(dot(a, b)/(norm(a)*norm(b)))

%--------------------- 12. fólia ---------------------------
A = [3 4 1; -1 2 5; 7 -2 3], 
b = [18 8 32]', rang = rank(A), det(A) 
x = A\b, A*x == b
%---------------------- 13-14. fólia ---------------------
A = xlsread('Keverés.xlsx','Mese1','b2:e5')
b = xlsread('Keverés.xlsx','Mese1','f2:f5')
det(A), x=A\b, A*x, b

%------------------ 15. fólia ------------------------

A = [10 -7 0; -3 2 6; 5 -1 5], b = [7; 4; 6]
Ab = [A b]

L1 = [1 0 0; 0.3 1 0; -0.5 0 1], Ab1 = L1*Ab
P = [1 0 0; 0 0 1; 0 1 0], P*Ab1
L2 = [1 0 0; 0 1 0; 0 0.1/2.5 1], Ab2 = L2*P*Ab1
% Mivel P inverze P és szorzat inverze az inverzek
% fordított sorrendű szorzata, és
% (P*inv(L1)*P*inv(L2))*(L2*P*L1*A*)*x = P*b  , ezért
LL = P*inv(L1)*P*inv(L2)
UU = L2*P*L1*A
% így LL*UU = P*A és LL*y = P*b majd U*x = y
LL*UU, P*A
y = LL\(P*b), x = UU\y, A*x, b

[L U PP] = lu(A)


%------------------- 18. fólia -------------------------
A = [3 4 1; -1 2 5; 2 6 6], b = [18 8 8]'
rank(A), rank([A b]), det(A)

% összefüggő eset ---- 19. fólia -------------------
b = [18 8 26]', rank(A), rank([A b])
y = pinv(A)*b, A*y, y_norma = norm(y)
% a kváziinverzzel kapott megoldás (a leg)rövidebb vektorhosszú, 
% nincs warning

% A független egyenletek számának lecsökkentése a rangig
% töröljük az utolsó egyenletet
A(3,:)=[], b(3)=[], rank(A), rank([A b])  

z = A\b, A*z, z_norma = norm(z)   % nem minimális a norma

x = pinv(A)*b, A*x, x_norma = norm(x)  % minimális normájú

v = A'/(A*A')*b  % ekkor a kváziinverzes megoldást kapjuk

% túlhatározott eset----------- 20. fólia ------------
A = [3 4 1; -1 2 5; 7 -2 3; 3 5 1], b = [18 8 32 6]' 
rank(A), rank([A b])
x = A\b, pinv(A)*b, (A'*A)\A'*b
% az utolsó a pszeudoinverz
[A*x b], tavolsag = norm(A*x-b)

%--------------- időmérés 21. fólia ---------------
n=5; 
c=[1 randi([1,2*n],[1,n-1])]; 
r=[1 randi([1,2*n],[1,n-1])];
T = toeplitz(c, r), b=[randi([1,2*n],[1,n])]'
X = [T\b inv(T)*b pinv(T)*b]
[b T*X]

m = [100 200 300 400 500];
t=[];
for n = m
 c=[1 randi([1,2*n],[1,n-1])]; 
 r=[1 randi([1,2*n],[1,n-1])];
 T = toeplitz(c, r); b=[randi([1,2*n],[1,n])]';
 tstart = tic; x = T\b; t(1,end+1) = toc(tstart);
 tstart = tic; x = inv(T)\b; t(2,end) = toc(tstart);
 tstart = tic; x = pinv(T)\b; t(3,end) = toc(tstart);
end
plot(m,t(1,:), m,t(2,:), m,t(3,:))
legend('balosztó','inverz', 'pszeudoinverz',2)
[['méret:     '; ...
  'balosztó   '; ...
  'inverz     '; ...
  'kváziinverz'] num2str([m; t],4)]

%---------------- vektornormák 23. fólia ---------------
x = (1:4)/5, 
norm1 = [norm(x,1), max(sum(abs(x)))]
norm2 = [norm(x), norm(x,2), sqrt(x*x') sqrt(dot(x,x))]
norminf = [norm(x,inf), max(abs(x))]
p=3; normp = [norm(x,p), sum(abs(x).^p)^(1/p)] % p valós

% --------------- mátrixnormák  24. fólia --------------
A=[4 5 6; 5 6 8; 1 7 7], 
norm1 = [norm(A,1) max(sum(abs(A)))] % maximális oszlopösszeg
norm2 = [norm(A,2) norm(A) max(svd(A))]  % műveletigényes!!
normInf = [norm(A,inf) max(sum(abs(A')))] % maximális sorösszeg
normFro = [norm(A,'fro') sqrt(sum(diag(A'*A)))]
% elemek négyzetösszegéből vont négyzetgyök

% ---------------- kondíciószám 24. fólia -------------
kond2 = [norm(A)*norm(inv(A)), cond(A), cond(A,2)]

% rcond(A): ha 1 közeli akkor jól kondícionált a mátrix
% ha pedig 0 közeli, akkor gyengén kondícionált
[1/(norm(A,1)*norm(inv(A),1)), rcond(A)]

% -------------- kondíciószám 26. fólia -------------
A=[4 5 6; 5 6 8; 1 7 7]
[cond(A), ...
sqrt(max(abs(eig(A'*A))))/sqrt(min(abs(eig(A'*A)))), ...
max(svd(A))/min(svd(A))]
% megjegyzés: négyzetes mátrixnál
% sort(sqrt(eig(A'*A)),'descend') és svd(A) megegyezik
% ---------- sajátérték, sajátvektor 27. fólia --------
a = poly(A), lambda = roots(a), poly(lambda)
[U Lambda] = eig(A), A*U, U*Lambda
U*Lambda*inv(U), A

% ------ sajátérték, sajátvektor 28. fólia --------
A = gallery(3)
a = poly(A), lambda = roots(a)  % karakt. polinommal
B = sym(A), syms l, det(B-l*eye(3))  

[U L] = eig(A)
U(:,1)=U(:,1)/U(1,1)
U(:,2)=U(:,2)/U(3,2)
U(:,3)=U(:,3)/U(3,3)*9

[U*L*inv(U), A]
% ------- szinguláris felbontás 29. fólia ------------
[X, S, V] = svd(A) % A(mxn) S(mxn) X(mxm) V(nxn) méretű
[A, X*S*V']  % svd felbontás
X'*X; V'*V  % X és V unitér (valós esetben ortogonális) mátrixok
% a pszeudoinverz és az svd kapcsolata:
[pinv(A) V*diag(1./diag(S))*X']


% érzékenységi elemzés 30. fólia --------------
A = gallery(5)
a = poly(A), lambda = roots(a)  % karakt. polinommal
B = sym(A), syms l, det(B-l*eye(5))

[U L] = eig(A), kondi_U = cond(U), e = diag(L)
plot(real(e), imag(e), 'r*', 0, 0, 'ko')
axis(0.1*[-1 1 -1 1]), axis square