% 1. gyakorlat: mátrix hatványozás
%%
%Perron-Frobenius mátrixok hatványozása -> diffúzió
%a Perron-Frobenius diffúziós rendszer mátrix alakban:
a=; %0.5>=mu>=0 % PF kritériuma mu>0.5-nél a rendszer elszáll
b=; %1-2mu, PF kritériuma
%mx=[b a 0 0; a b a 0; 0 a b a; 0 0 b a]; %kis mx-ot kézzel is be lehet írni, de a nagyot már nem
meret=20;
hatvanyok=100;%800;
mx=; %a diffúziós  mátrix elõállítás
%... %Neumann (fix, zero-flux) határ, ~ anyagmegmaradás
%... %Neumann (fix, zero-flux) határ, ~ anyagmegmaradás
[sajatvectorok, sajatertekek]=eig(mx); %sajátérték számítás

cc_vect_init=zeros(meret,1); %kezdeti vektor, amire a mátrix hat
cc_vect_init(8:13)=1; % kiindulási koncentráció eloszlás
mxhatv=mx^hatvanyok; %hatványozom a mátrixot
[hatvsajatvectorok, hatvsajatertekek]=eig(mxhatv);
cc_vect_final=mxhatv*cc_vect_init; %végsõ cc

%%
%lépésenként számolva a folyamatot:
cc_vect=cc_vect_init; 
se=zeros(meret,hatvanyok); %ide tárolom minden iterációban a sajátértékeket
figure;

for sz=1:hatvanyok 
mxhatv=; %hatványozom a Perron-Frobenius mátrixot
sajatertekek=; %a hatványmátrix sajátértékei
se(:,sz)=sajatertekek; %sajátértékeket mindig eltárolom
cc_vect=mx*cc_vect; %aktuális koncentrácó értékek

%a diffúziós folyamat ábrázolása 
subplot(6,4,[13:15,17:19,21:23]);
imagesc(mxhatv) %ábrázolom a mátrixot
caxis([0,1]); %a színskála határainak fixálása (egyébként automatikusan módosul a skála a legnagyobb és legkisebb érték közé)
title(['hatvány: ', num2str(sz)])
subplot(6,4,[4,8,12]);
imagesc(cc_vect_init) %kezdeti cc értékek
caxis([0,1]);
title('kezdeti koncentráció értékek')
subplot(6,4,[16,20,24]);
imagesc(cc_vect) %ábrázolom a cc változásokat
caxis([0,1]);
title(['koncentráció értékek ', num2str(sz), ' lépés diffúzió után'])

waitforbuttonpress; %kattintásra lép egyet a ciklus
pause(0.1) %késleltetés az iterációk között
end

figure;
plot(1:hatvanyok, se, '.b'); %ábrázolom a sajátértékeket, ezek jellemzik a mátrixot (az abszolútértékben maximális leginkább) 
xlim([1, hatvanyok]);
xlabel('hatványok');
ylabel('sajátértékek');
title('A mátrix sajátértékeinek alakulása a hatványozással.');

% legnagyobb sé-hez csupa 1-es sv tartozik, ez lesz domináns
% hatványozással~ idõ múlással, ezért van kiegyenlítõdés a folyamat elõre haladtával
% a kis sé-hez tartozó sv gyakran vált elõjelet, oszcillál, de ez hamar
% eltûnik, lecseng a hatványozásnál
%a<0.5, akkor max sé 1 (anyagmegmaradásos esetben), minimális -1-nél nagyobb, a max dominál hatványozásnál a többi eltûnik hatványozással
%(ha a=0.5, akkor +- -ban 1 erõs hullám szõrösödik ez 2 tag marad meg legtovább hatványozással, szimmetria 
%a>0.5, akkor max sé 1 (anyagmegmaradásos esetben), minimális -1-nél kisebb, ez száll el hatványozással, oszcillálva, ez a domináns
%

%%
% stabilizálódó korfa: a mátrix határozza meg
L= %Leslie mátrix
% más mátrixok: stabil populáció, népesség robbanás, elöregedõ/fogyó közösség 
pop_kezd=[1 0 0 0]' %mi történik más kezdeti értéknél?
se=
hatvanyok=40;

figure;
for i=1:hatvanyok
    i
   L_hatv= %aktuális generáció
   pop_akt= %populáció aktuális értéke
   se_hatv= %sajátértékek
   
    subplot(6,4,[13:15,17:19,21:23]);
    imagesc(L_hatv) %ábrázolom a mátrixot 
    title(['generáció: ', num2str(i)])
    subplot(6,4,[4,8,12]);
    imagesc(pop_kezd) %kezdeti értékek
    title('kezdeti populáció (korfa)')
    subplot(6,4,[16,20,24]);
    imagesc(pop_akt) %ábrázolom a vektor értékek változását
    title(['korfa ', num2str(i), ' generáció után'])
   
   waitforbuttonpress;
   pause(0.1)
end

%%
% nem minden mátrix egyenlítõdik ki!
mx=; %hatványozáskor ciklikusan magába visszatérõ mx
vect_init=; %kezdeti érték
se=; %sajátértékek

figure;

for i=1:20
    i
   mxhatv= %aktuális hatvány
   vect_akt= %a vektor aktuális értéke
   se_hatv= %sajátértékek
   
    subplot(6,4,[13:15,17:19,21:23]);
    imagesc(mxhatv) %ábrázolom a mátrixot
    caxis([0,1]); 
    title(['hatvány: ', num2str(i)])
    subplot(6,4,[4,8,12]);
    imagesc(vect_init) %kezdeti értékek
    title('kezdeti értékek')
    subplot(6,4,[16,20,24]);
    imagesc(vect_akt) %ábrázolom a vektor értékek változását
    title(['új értékek ', num2str(i), ' lépés után'])
    
   waitforbuttonpress;
   pause(0.1)
end

%%
%az zaj (epsilon) pozíciója fontos
%nagyobb rendszerekben egyre jobban: "curse of dimensions"
n = 10;
epsilon=10^(-3); %valamilyen kis érték (zaj a rendszerben)
M=; %mátrix
%különbözõ pozíciókban kis epsilon zaj a rendszerben:
%... %adj epsilon zajt a rendszerhez az különbözõ pozíciókban
[v, s] =; %a rendszer sajátértékei
s=diag(s)
figure;
plot(real(s), imag(s), 'o'); %sajátértékek ábrázolása: képzetes rész a valós rész függvényében 
xlabel('Valós rész');
ylabel('Képzetes rész');

%{
%a rendszer egyenletének megoldása:
syms x_sym %szimbolikus változó "analitikus" megoldáshoz
M_sym=sym(M);
char_polinom=charpoly(M_sym, x_sym) %a rendszer karakterisztikus polinomja 
x_megoldas=solve(char_polinom); %a rendszer megoldása
double(x_megoldas) %megoldások számként kiíratva
%}
