%3. gyakorlat: Lotka-Volterra rendszer
%dinamikai rendszer, az egyenlet rajzol:
%egyenletek: dx=x*(8-x-2y); dy=y*(5-x-y);
figure;
%ábrázolási határok:
max_pontx=9;%20;%
max_ponty=6;%20;%
min_pont=0.;

%
lepes=0.1;% lépésköz
x=min_pont:lepes:max_pontx;%
y=min_pont:lepes:max_ponty;%
%}

%
% rácsot képez x,y alapján:
[X,Y]=meshgrid(x,y);
%a rácspontokban pontokban kiszámolom a vektormezõ gradiensét
dx=X.*(8-X-2*Y); 
dy=Y.*(5-X-Y);
dx=dx./ sqrt(dx.^2+dy.^2); %vektormezõ lenormálva
dy=dy./sqrt(dx.^2+dy.^2);
%1.vektormezõ kirajzolás:
quiver(X,Y,dx, dy, 'k');
% az x y rácspontokon az x és az y irányba változást számolja ki és az
% eredményül kapott vektort ábrázolja az adott pontból kiindulva,
%}

%2. egyenesek kirajzolása, ahol dx=0 vagy dy=0
y0=4-x/2; % 8-x-2*y=0 egyenes egyenlete
y1=5-x; % 5-x-y=0 egyenes egyenlete 
hold on; %egymásra akarom rajzolni a trajektóriákat
plot(x, y0,'b', zeros(1, length(y)), y, 'b', x, y1, 'r', x, zeros(1, length(x)), 'r'); %ábrázolom azokat az egyeneseket ahol dx=0 vagy dy=0
xlim([-0.1, max_pontx]); %ábrázolás határa
ylim([-0.1, max_ponty]);
title('egyenletek: dx=x*(8-x-2*y); dy=y*(5-x-y)'); %ábra cím
xlabel('x faj mennyisége'); %tengelyfeliratok
ylabel('y faj mennyisége');
%legend('dx=0, mert x=0', 'dx=0, mert 8-x-2*y=0', 'dy=0, mert y=0', 'dy=0, mert 5-x-y=0'); %jelmagyarázat
legend('vektormezõ gradiense adott pontokban', 'dx=0, mert x=0', 'dx=0, mert 8-x-2*y=0', 'dy=0, mert y=0', 'dy=0, mert 5-x-y=0');

%
%3. az X, Y rács pontjaiból indulok, egyszerû, de lassú:
sx=reshape(X, 1, length(x)*length(y)); %kezdeti x pontok sorba rendezve
sy=reshape(Y, 1, length(x)*length(y)); %kezdeti y pontok sorba rendezve


%4.a egy egyenesrõl indítok trajektóriákat:
N=100; %hány pontot számolok ki egy egyenes mentén
sxk=0;
syk=1;
sxv=1;
syv=0;
sx=linspace(sxk,sxv,N);
sy=linspace(syk,syv,N);

%
%4.b 2 egyenesrõl indítok trajektóriákat:
N=100;
sxk=0;
syk=1;
sxv=1;
syv=0;

sx2k=0;
sy2k=6;
sx2v=9;
sy2v=0;

%2 pont készletet fûzök egymáshoz
sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N)];

%
%4.c az ábra 4 határoló egyeneseirõl indítok trajektóriákat:
sxk=0.1;
syk=0.1;
sxv=0.1;
syv=max_ponty;

sx2k=0.1;
sy2k=0.1;
sx2v=max_pontx;
sy2v=0.1;

sx3k=max_pontx;
sy3k=0.1;
sx3v=max_pontx;
sy3v=max_ponty;

sx4k=0.1;
sy4k=max_ponty;
sx4v=max_pontx;
sy4v=max_ponty;

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)];

%
%5. [0,0] környékérõl indítok trajektóriákat:
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));

%
%6. teljes szeparatrix kirajzolása: nagyobb max_pontok is lehetnek
sx=[sx, sx2];
sy=[sy, sy2];
%}

%7. idõben vissza: f(x)=-f(x)
%a trajektória végpontjait adom meg, a futás során visszafele halad a folyamat a kiindulási hely felé
%{
sxk=1.95;
syk=3;
sxv=2;
syv=2.95;

sx2k=1.97;
sy2k=3.05;
sx2v=2.05;
sy2v=3;

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

%
tint=[0 10]; %futási idõsáv
%minden kiindulási pontra kiszámolom a diffegyenletet
for i=1:length(sx) 
    xy0=[sx(i); sy(i)]; %mindig az adott x,y pontból számolom ki a diffegyenletet (a hold on be van kapcsolva így mindent egyszerre rajzolok ki)
    [t, xy] = ode45(@(t, xy) [xy(1)*(8-xy(1)-2*xy(2)); xy(2)*(5-xy(1)-xy(2))],...
        tint, xy0); %az egyenletet kiszámolom
    %[t2, xy2] = ode45(@(t, xy) -[xy(1)*(8-xy(1)-2*xy(2)); xy(2)*(5-xy(1)-xy(2))],tint,xy0); %az idõ visszafelé: f(x)=-f(x)
    %xy=[flipud(xy2); xy]; %2 kimenetet (idõ oda, vissza) 1 folytonos tajeltóriává fûzöm össze (az idõ visszafelét elõször megfordítom)
    
    %plot(xy(:,1),xy(:,2), 'g');
    
    %
    %a szeparatrix ábrázolása: ha az egyik stabil pontba fut be a mg akkor az egyik, ha a másikba, akkor a másik színnel ábrázolom a trajektóriát
    if abs(xy(end,1)-0)<1 %melyik egyensúlyi helyzethez van adott közel (az x tengelyen) a trajektória végpontja
        plot(xy(:,1),xy(:,2), 'g');
    elseif abs(xy(end,1)-8)<1
        plot(xy(:,1),xy(:,2), 'c');
    end
    %}  
end
%}

hold off;
