% 6. gyakorlat
% Lotka-Volterra populáció dinamikai rendszerek
% általános alak: 
% dx=x*(xb-xd*x+xi*y)
% dy=y*(yb-yd*y+yi*x)
% x, y = fajok mennyisége
% xb, yb =pl. szaporodási/halálozási ráta ráta
% xd, yd =pl. erõforrás korlát
% xi, yi = interakció x és y közt

% bifurkáció: a rendszer minõségileg változik meg 
% pl itt: nem csak kicsit odébb megy az egyensúlyi pont, nem csak megváltoznak a paraméterei (mennyiségi, kvantitatív változások), 
% hanem minõségileg (kvalitatívan) másfajta egynsúlyi pont jön létre vagy szûnik meg
% pl: stabil csomóból nyeregpont vagy periodikus pálya lesz vagy vonzó pontból nyereg 

% Feladat: vizsgáljátok meg az alábbi rendszer viselkedését különbözõ "mu" paraméterek mellett!  
% dx=x*(1-x/2-y) 
% dy=y*(mu-y-x)
% Milyen típusú rendszerrõl van szó? (figyeld a keresztcsatolásokat)
% Legyen 0<mu<3.
% Hol (milyen mu paraméter(ek) mellett) van bifurkációs pont?
% Milyen típusú egyensúlyi helyzetek képzõdnek, tûnnek el, alakulnak át egymásba?
% Hogyan függ a rendszer viselkedése a mu paramétertõl (az y faj szaporodási rátájától)?
% Biológiailag/ökológiailag mit jelentenek az eredmények?

% paraméterek:
mu=1;
xb=1;
xd=-0.5;
xi=-1;
yb=mu;
yd=-1;
yi=-1;

%határok:
tint=[0 10];%idõintervallum
max_pontx=2.5;
max_ponty=3.5;
min_pont=0;
lepes=0.25; %felbontás
x=min_pont:lepes:max_pontx;% vizsgált pontok
y=min_pont:lepes:max_ponty;% vizsgált pontok

%az ábra alsó részén sûrûbben indítsuk a trajeltóriákat (a vizualizációhoz)
max_pontx2=1;
max_ponty2=0.5;
lepes2=0.05; %felbontás
x2=min_pont:lepes2:max_pontx2;% vizsgált pontok
y2=min_pont:lepes2:max_ponty2;% vizsgált pontok

% rácspontok meghatározása x, y vektorokból:
[X,Y]=meshgrid(x,y);
[X2,Y2]=meshgrid(x2,y2);

% több trajektória kirajzolása az X,Y pontokból indulva:
sx=[reshape(X, 1, length(x)*length(y)), reshape(X2, 1, length(x2)*length(y2))]; % a pontok x koordinátái egy sorvektorba rendezve (könnyebb végigiterálni rajta)
sy=[reshape(Y, 1, length(x)*length(y)), reshape(Y2, 1, length(x2)*length(y2))]; % a pontok y koordinátái egy sorvektorba rendezve (könnyebb végigiterálni rajta)

figure;

for mu=0:0.1:3
yb=mu;
%
% a vektormezõ gradiensének kiszámítása minden rácspontban
dx=X.*(xb+xd*X+xi*Y); 
dy=Y.*(yb+yd*Y+yi*X);
l=sqrt(dx.^2+dy.^2);
dx=dx./l; % 1-re normálás
dy=dy./l;
% a vektormezõ ábrázolása (quiver függvény):
quiver(X,Y,dx, dy, 'k');
% x és y változásának kiszámítása minden X,Y rácspontban, és a kapott vektor ábrázolása az X,Y pontból kiindulva

% dx=0 vagy dy=0 egyenesek ábrázolása
y0=(-xb-xd*x)/xi; %xb-xd*x+xi*y=0 egyenes egyenlete
x0=(-yb-yd*y)/yi;%yb-yd*y+yi*x=0 egyenes egyenlete
hold on; %trajektóriák egymásra rajzolásához
plot(x, y0,'b', zeros(1, length(y)), y, 'b', x0, y, 'r', x, zeros(1, length(x)), 'r'); %dx=0 vagy dy=0 egyenesek ábrázolása
xlim([-0.1, max_pontx]); %határok a megjelenítéshez
ylim([-0.1, max_ponty]);
title(['egyenletek: dx=x*(',num2str(xb),'+',num2str(xd),'*x+',num2str(xi),'*y); dy=y*(',num2str(yb),'+',num2str(yd),'*y+',num2str(yi),'*x)']); %ábra cím
xlabel('x mennyisége'); %tengelyfeliratok
ylabel('y mennyisége');
legend('a vektormezõ gradiense az adott pontokban', ['dx=0 <-  ',num2str(xb),'+',num2str(xd),'*x+',num2str(xi),'*y=0'],'dx=0 <-  x=0', ['dy=0 <-  ',num2str(yb),'+',num2str(yd),'*y+',num2str(yi),'*x=0'], 'dy=0 <-  y=0');

% több trajektória kirajzolása az X,Y pontokból indulva:
for i=1:length(sx) 
    xy0=[sx(i); sy(i)]; % egyenlet kiszámítása az adott kezdeti pontból (a hold on engedi a trajektóriák egymásra rajzolását)
    [t, xy] = ode45(@(t, xy) [xy(1)*(xb+xd*xy(1)+xi*xy(2)); xy(2)*(yb+yd*xy(2)+yi*xy(1))], tint, xy0); %egyenlet kiszámítása
    
    h=plot(xy(:,1),xy(:,2), 'g'); 
    h.Annotation.LegendInformation.IconDisplayStyle = 'off'; %a legend-ben ne jelenjenek meg
end

hold off;

waitforbuttonpress;
end