% 5. 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

%
% versengés 
% dx=x*(8-x-3y) 
% dy=y*(5-y-x)
% paraméterek:
xb=8;
xd=-1;
xi=-3;
yb=5;
yd=-1;
yi=-1;
%határok:
max_pontx=8;
max_ponty=6;
min_pont=0;
lepes=0.5; %felbontás
tint=[0 10];%idõintervallum
%}

x=min_pont:lepes:max_pontx;% vizsgált pontok
y=min_pont:lepes:max_ponty;% vizsgált pontok

figure;

%
% 1.: vektormezõ
% rácspontok meghatározása x, y vektorokból:
[X,Y]=meshgrid(x,y);
% 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
%}

%2. 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');

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

%{
%az ábra 4 határoló egyeneseirõl indítok trajektóriákat:
N=100;
sxk=0.1;
syk=0.1;
sxv=0.1;
syv=y(end);

sx2k=0.1;
sy2k=0.1;
sx2v=x(end);
sy2v=0.1;

sx3k=x(end);
sy3k=0.1;
sx3v=x(end);
sy3v=y(end);

sx4k=0.1;
sy4k=y(end);
sx4v=x(end);
sy4v=y(end);

sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N), linspace(sx3k,sx3v,N), linspace(sx4k,sx4v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N), linspace(sy3k,sy3v,N), linspace(sy4k,sy4v,N)];
%}

%{
% trajektóriák indítása sûrûbben az origó környékérõl
lepes=0.01;
x=min_pont:lepes:0.5;
y=min_pont:lepes:0.5;
[X,Y]=meshgrid(x,y);
sx2=reshape(X, 1, length(x)*length(y));
sy2=reshape(Y, 1, length(x)*length(y));

%sx=sx2;
%sy=sy2;
sx=[sx, sx2];
sy=[sy, sy2];
%}

%
%trajektóriák indítása sûrûn a nyeregpont körmékérõl (2 egyenerõl)
%idõ megfordítás: f(x)=-f(x)
% a trajektória végpontja adott; a dinamika a "visszafelé zajlik le" a kezdõpont felé  
N=100;
sxk=3.43;
syk=1.52;
sxv=3.53;
syv=1.48;

sx2k=3.46;
sy2k=1.54;
sx2v=3.56;
sy2v=1.5;

sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N)];
%}

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
    [t2, xy2] = ode45(@(t, xy) -[xy(1)*(xb+xd*xy(1)+xi*xy(2)); xy(2)*(yb+yd*xy(2)+yi*xy(1))],tint,xy0); %idõmegfordítás: f(x)=-f(x)
    xy=[flipud(xy2); xy]; %  a két irányú megoldás egy folytonos trajektóriába összefûzve (az idõben visszafelé számolt megfordítás után)
    
    %h=plot(xy(:,1),xy(:,2), 'g');
    
    %
    % a szeparatrix (2 egyensúlyi pontot szeparáló görbe) kirajzolás: ha az egyik egyensúlyi ponthoz van közel a trajektória végpontja, az egyik, ha a másikhoz, akkor a másik színnel ábrázolom a pályát 
    if abs(xy(end,1)-0)<1 % a trajektória végpontjának y koordinátája megyik egyensúlyi ponthoz esik közel
        h=plot(xy(:,1),xy(:,2), 'm');
    elseif abs(xy(end,1)-8)<1
        h=plot(xy(:,1),xy(:,2), 'c');
    end
    %}  
    h.Annotation.LegendInformation.IconDisplayStyle = 'off'; %a legend-ben ne jelenjenek meg
end
%}

hold off;