%2. gyakorlat: diffúzió

dt = 0.1; %idõbeli lépésköz
dx = 1; %térbeli lépésköz
D=1; %diffúzió sebessége %ezt változtassuk most
nt = 1500; %idõtartomány vége, eddig futtatok

%
x = 0:dx:500; %viszgált pontok egy egyenesen
nx = length(x); %vizsgált pontok száma
u = zeros(nx,1); % kezdeti értékek 0-ák

u(250:255)=1; %csúcs indul, sûrûség függvény, eloszlásfüggvény deriváltja
%u(260:265)=1; %két csúcs induláskor
%u(250:end)=1;% eloszlásfüggvény,

figure;
%grid on;

for t = 1:nt %minden idõpillanatban
    u2 = u; %u2=értékek az elõzõ idõplillanatban; u=új értékek
    for i = 2:nx-1 %minden pontban kiszámolom a hullám aktuális értéké   
       %minden pontra kiszámoljuk az új értékét:               
       u(i) = u2(i) + D*dt/dx^2*(u2(i+1)-2*u2(i)+u2(i-1)); %dU/dt=D*dU^2/dx^2 %diffúzió, szétfolyik
    end
    
    plot(u); %aktuális értékeket kirajzol
    grid(gca, 'minor'); % sûrû rácsot tesz a képre
    ylim([-0.5 1.5]); %ábrázolás határai
    xlim([100, 400])
    title(['lépés: ', num2str(t)]); %kiírítja az ábrára, hogy hányadik lépésnél vagyok
    %waitforbuttonpress; %gombnyomásra lép egyet a ciklus
     pause(0.01); %ilyen idõközönként rajzolja ki az új állapotot
end
%}

%
% a rendszer elszáll, ha:
%D*2*dt/dx^2>1
%Perron-Frobenius-ként nézve:
a=D*dt/dx^2;
b=1-D*2*dt/dx^2;
meret=4;
hatvanyok=300;
se=zeros(1,hatvanyok);
PFmx=diag(b*ones(meret,1),0)+diag(a*ones(meret-1,1),1)+diag(a*ones(meret-1,1),-1) %Perron-Flobenius mátrix elõállítás
%PFmx(1,1)=0.9;
%PFmx(end,end)=0.9;
%PFmx
sajatertekek=eig(PFmx) %sajátérték számítás
se(1)=max(abs(sajatertekek)); % (abszolútértékben) maximális sé jellemzi a mátrixot, eltárolom egy vektorban
inv(PFmx)
%
for sz=2:hatvanyok 
PFmxhatv=PFmx^sz; %hatvanyozom a mátrixot
sajatertekek=eig(PFmxhatv);
se(sz)=max(abs(sajatertekek)); %max sajátértéket mindig eltárolom
end
PFmxhatv
%eig(PFmxhatv)
%inv(PFmxhatv)
figure;
plot(se); %ábrázolom a maximális sajátértékeket, ezek jellemzik a mátrixot
ylim([0 10]);
xlim([1, hatvanyok]);

%}