% 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.1; %0.5>=mu>=0 % PF kritériuma mu>0.5-nél a rendszer elszáll
b=1-2*a; %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=diag(b*ones(meret,1),0)+diag(a*ones(meret-1,1),1)+diag(a*ones(meret-1,1),-1); %a diffúziós  mátrix elõállítás
mx(1,1)=1-a; %Neumann (fix, zero-flux) határ, ~ anyagmegmaradás
mx(end,end)=1-a; %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);
figure;

for sz=1:hatvanyok 
mxhatv=mx^sz; %hatványozom a Perron-Frobenius mátrixot
sajatertekek=eig(mxhatv);
se(:,sz)=sajatertekek; %sajátértékeket mindig eltárolom
cc_vect=mx*cc_vect;

%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=[0 1 1 1; 1 0 0 0; 0 1 0 0; 0 0 1 0] %általános Leslie mátrix
L=[0 1/2 0 1/2; 1 0 0 0; 0 1 0 0; 0 0 1 0] % L mx: irreducibilis, Pl: Csertölgy csak minden 2. évben terem
%L=[0 0 2 6; 1/2 0 0 0; 0 1/2 0 0; 0 0 1/3 0] % L mx: primitív, Pl: idõsebb egyedek szaporodnak csak 
%L=[0 1 0.75 0.5; 0.5 0 0 0; 0 1 0 0; 0 0 0.5 0] %sok születés, fiatal egyedek szaporábbak, gyengébb életkilátások: hagyományos társadalom 
%L=[0 1 0.75 0.5; 1 0 0 0; 0 1 0 0; 0 0 0.5 0] %sok születés, fiatal egyedek szaporábbak, jobb életkörülmények: népességrobbanás 
%L=[0 0.25 0.25 0.25; 1 0 0 0; 0 1 0 0; 0 0 1 0] %kevés születés, bár sokáig aktív egyedek: elöregedõ társadalom 
R=L(2,1)*L(1,2)+L(3,2)*L(2,1)*L(1,3)+L(4,3)*L(3,2)*L(2,1)*L(1,4)%stabil populáció, ha R=s1*b1+s2*s1*b2+s3*s2*s1*b3=1 <-> domináns se=1
pop_kezd=[1 0 0 0]'
%pop_kezd=rand(4,1) %tetszõleges kezdeti pouláció -> ez után is stabilizálódás 
se=eig(L)
hatvanyok=40;

figure;
for i=1:hatvanyok
    i
   L_hatv=L^i
   pop_akt=L_hatv*pop_kezd
   se_hatv=eig(L_hatv)
   
    subplot(6,4,[13:15,17:19,21:23]);
    imagesc(L_hatv) %ábrázolom a mátrixot 
    title(['generáció: ', num2str(i)])
    colorbar;
    subplot(6,4,[4,8,12]);
    imagesc(pop_kezd) %kezdeti értékek
    title('kezdeti populáció (korfa)')
    colorbar;
    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'])
    colorbar;
   
   waitforbuttonpress;
   pause(0.1)
end

%%
% nem minden mátrix egyenlítõdik ki!
mx=[0 0 0 1; 1 0 0 0; 0 1 0 0; 0 0 1 0] %hatványozáskor ciklikusan magába visszatérõ mx %Pl:cserebogarak 4 évesen szaporodnak 1szer, a populáció szinkronban 
vect_init=[1; 2; 3; 4];
se=eig(mx)
hatvanyok=20;

figure;

for i=1:hatvanyok
    i
   mxhatv=mx^i
   vect_akt=mxhatv*vect_init
   se_hatv=eig(mxhatv)
   
    subplot(6,4,[13:15,17:19,21:23]);
    imagesc(mxhatv) %ábrázolom a mátrixot
    caxis([0,1]); 
    title(['hatvány: ', num2str(i)])
    colorbar;
    subplot(6,4,[4,8,12]);
    imagesc(vect_init) %kezdeti értékek
    title('kezdeti értékek')
    colorbar;
    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'])
    colorbar;

    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=diag(ones(n,1)) + diag(ones(n-1,1),1);
%különbözõ pozíciókban kis epsilon zaj a rendszerben:
%M(n,1) = epsilon;
M(n,n-1) = epsilon;
[v, s] = eig(M);
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
%}
