%% Válogatott modellek a matematikai biológiában
%% 1. gyakorlat
%% Perron-Frobenius mátrixok hatványozása
% Diffúzió a Perron-Frobenius diffúziós rendszer mátrix alakban.
% Tekintsük $u: [0,\infty]\times[0,L] \to \mathbb{R}$függvényt, mely kielégíti 
% a $u'_t(t,x) = u''_{xx}(t,x)$ diffúziós egyenletet a 
%% 
% # Carl Neumann féle ($u'_t(t,0) = u'_t(t,L) = 0$) peremfeltétellel,
% # Dirichlet féle $u(t,0) = u(t,L) = 0$ peremfeltétellel,
% # periodikus $u(t,0) = u(t,L)$ és $u'_t(t,0) = u'_t(t,L)$ peremfeltételekkel.
%% 
% Praktikus okokból csak véges időhorizonton vizsgáljuk a függvény viselkedését, 
% vagyis $u: [0,T]\times[0,L] \to \mathbb{R}$. 
% 
% Az időhorizont amelyen a dinamikát vizsgálom jelölje $N_T = T/\tau \in \mathbb 
% N$

nT = 160;
%% 
% Legyen $u_{i,j} = u(i\tau,jh)$, $i = 0,\dots,N_T, j = 0,\dots,N_L$, úgy, hogy 
% $N_T \cdot \tau = T$és $N_L \cdot h = L$.
% 
% *Figyelem!* $\mu = \frac{\tau}{h^2}$-t úgy kell megválasztanunk, hogy $0 < 
% \mu \le \frac{1}{2}$.
% 
% (Vajon mi történne, ha $\mu$-re ez nem teljesülne?)

mu = 0.6;
%% 
% A térdimenzió (x) mentén $N_L + 1$ ponttal a teret $N_L$ egyenlő részre osztottuk, 
% ezért az $A_N$ mátrix mérete is $N_L \times N_L$.
% 
% _A kódban_ $N_L$_-et jelölje:_

nL = 50;
% A diffúziós mátrix elõállítása
% A mátrixot a |toeplitz| függvénnyel könnyedén létre lehet hozni. A |toeplitz| 
% függvény diagonálisan ismétli a megadott első sort és oszlopot, lásd, |help 
% toeplitz|. _(A kódban diffúziós mátrixot jelölje |A|.)_

A_elso_sora = [1-2*mu mu zeros(1,nL-2)];
A = toeplitz(A_elso_sora,A_elso_sora);
%% 
% A peremfeltételekre külön figyelni kell a mátrix sarkaiban.
% Neumann peremfeltétel
% $A_N = \pmatrix{ 1 - \mu & \mu & 0 & \dots & 0 \cr \mu & 1 - 2\mu & \mu & 
% \dots & 0 \cr \dots \cr 0 & 0 & 0 & \dots & 1 - \mu }$. 
% 
% Az $N$alsó index az $A_N$-ben a Carl Neumann féle peremfeltételre utal. Ez 
% a peremfeltétel megörzi a kezdeti anyag/energia mennyiségét.

A_N = A;
A_N(1,1) = 1-mu;
A_N(end,end) = 1-mu;
% Periodikus peremfeltétel
% $A_p = \pmatrix{ 1 - 2\mu & \mu & 0 & \dots & \mu \cr \mu & 1 - 2\mu & \mu 
% & \dots & 0 \cr \dots \cr \mu & 0 & 0 & \dots & 1 - 2\mu }$. 

A_p = A;
A_p(1,end) = mu;
A_p(end,1) = mu;
% Dirichlet peremfeltétel
% $A_D = \pmatrix{ 1 & 0 & 0 & \dots & 0 \cr \mu & 1 - 2\mu & \mu & \dots & 
% 0 \cr \dots \cr 0 & 0 & 0 & \dots & 1 }$. 

A_D = A;
A_D(1,1) = 1;
A_D(1,2) = 0;
A_D(end,end) = 1;
A_D(end,end-1) = 0;
%% 
% Kiválasztjuk a peremfeltételt.

A = A_p;
%% 
% Számítsik ki az $A$ mátrix sajátértékeit (|D|) és sajátvektorait (|S|).

[S,D] = eig(A);
Eigenvalues = diag(D)'
%% 
% Domináns sajátérték:

[r,idx] = max(Eigenvalues)
%% 
% Domináns sajátértékhez tartozó sajátvektor (vegyük észre a sajátvektor lényegében 
% a normált "csupa 1-es"):

S(:,idx)
%% 
% Kezdeti feltétel (pl.kiindulási koncentráció eloszlás, vagy rúd hőmérsékletének 
% kezdeti eloszlása, stb.), jelölje $x_0$.

x0 = zeros(nL,1);
% x0(8:13) = 1;
x0(13:18) = 1;
%% 
% Számoljuk ki $x_T$-t anélkül, hogy az köztes időben kapott $x_1, x_2, \dots, 
% x_{T-1}$értékeket kiszámolnám, tudva, hogy $x_k = A_N^k \, x_0$.

A_ad_nT = A^nT;
% [S_T,D_T] = eig(A_ad_nT);
x_nT = A_ad_nT*x0;
%% 
% Számoljuk ki lépésenként is egyenként $x_1, x_2, \dots, x_{T-1}, x_T$-t, közben 
% azt is vizsgáljuk, hogy a mátrixhatványok sorában miképp alakulnak a sajátértékek.

x = x0; 
eigs = zeros(nL,nT);
U = zeros(nL,nT);
fig = figure(123);
fig.Visible = 'on';

Tl = tiledlayout(2,4,"TileSpacing","tight","Padding","tight","TileIndexing","rowmajor");

ax1 = nexttile(4); hold on, grid on, box on
Im0 = imagesc(x0);
clim([0,1])
axis tight
title('$x_0$','Interpreter','latex','FontSize',12)

ax2 = nexttile(5,[1,3]); hold on, grid on, box on;
ImA = imagesc(A);
axis tight
Tl = title("$L^k,~ \makebox[2cm]{$k = 1\hfill$}$",'Interpreter','latex','FontSize',12);

ax3 = nexttile(8); hold on, grid on, box on;
Imk = imagesc(x0);
clim([0,1])
axis tight
title('$x_k = L^k x_0$','Interpreter','latex','FontSize',12)

Ak = eye(nL);
xk = x0;
for k = 1:nT
    Ak = Ak*A;
    xk = A*xk;
    eigs(:,k) = eig(Ak);
    U(:,k) = xk;
    
    ImA.CData = Ak;
    Imk.CData = xk;

    Tl.String = "$L_k,~ k = " + num2str(k) + "$";
    drawnow
end
%% 
% Most ábrázolom a mátrixhatványok sajátértékeinek alakulását,

figure;
plot(1:nT, real(eigs), '.b');
xlim([1, nT]);
xlabel('hatványok');
ylabel('sajátértékek');
title('A mátrix sajátértékeinek alakulása a hatványozással.');
%%
surf(U)