% 7. gyakorlat:
% általános Lotka-Volterra modell
%dx=x*(a1-b1*x+c1*y)
%dy=y*(a2-b2*y+c2*x)
% az egyensúlyban:
% a kezdeti állapottól függõen vagy az egyik vagy a másik marad meg [ 3  1 -4  1 0.5 -1] versengés

%paraméterek:
a1=3; % x szaporodási rátája %[0 4] 
b1=1; % x erõforráskorlátja (természetes halálozása) %[0.1 10]
c1=-4; % keresztcsatolás: az y hatása az x-re %[-10 10]
a2=1; % y szaporodási rátája %[0 4] 
b2=0.5; % y erõforráskorlátja (természetes halálozása) %[0.1 10]
c2=-1; % keresztcsatolás: az x hatása az y-ra %[-10 10]

%vektormezõ ábrázolása:
x=0:0.1:4;
y=0:0.1:4;
[X,Y]=meshgrid(x,y);
dx=X.*(a1-b1*X+c1*Y);
dy=Y.*(a2-b2*Y+c2*X);
%normálás:
dx=dx./ sqrt(dx.^2+dy.^2);
dy=dy./sqrt(dx.^2+dy.^2);

% 0 gradiensû egyenesek (dx=0; dy=0) ábrázolása:
dx0_1=zeros(1, length(y));
dx0_2=(-a1+b1*x)/c1;
dy0_1=zeros(1, length(x));
dy0_2=(a2+c2*x)/b2;

%ábrázolás:
figure;
%vektormezõ és a trajektória:
quiver(X,Y,dx,dy, 0.5, 'k');
hold on;
plot(dx0_1, y, 'r', x, dx0_2, 'r', x, dy0_1, 'b', x, dy0_2, 'b');

%trajektóriák kirajzolása:
%
t_null=0;
t_max=10;
N=100;
%az ábra 4 határoló egyeneseirõl indítok trajektóriákat:
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);

%az origó közelében sûrûbben veszem a pontokat
sx5k=0.01;
sy5k=0.01;
sx5v=0.01;
sy5v=0.5;

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

tolerancia=1;
for i=1:length(sx)
x_null=sx(i);
y_null=sy(i);
egyenlet=@(t,xy)[xy(1)*(a1-b1*xy(1)+c1*xy(2));xy(2)*(a2-b2*xy(2)+c2*xy(1))];
[t, xy]=ode45(egyenlet, [t_null, t_max], [x_null, y_null]);
%plot(xy(:,1), xy(:,2), 'g');

%az szerint színezem a pályákat, hogy melyik egyensúlyi ponthoz vannak közel
%
if abs(xy(end,1)-3)<tolerancia
    plot(xy(:,1), xy(:,2), 'c');
elseif abs(xy(end,2)-2)<tolerancia
    plot(xy(:,1), xy(:,2), 'm');
end
end
%}

%{
%szeparatrix ábrázolása: 2 egyenesen veszek fel pontokat közel a nyeregponthoz
t_null=0;
t_max=100;
N=100;

sx6k=0.71;
sy6k=0.58;
sx6v=0.735;
sy6v=0.565;

sx7k=0.69;
sy7k=0.575;
sx7v=0.72;
sy7v=0.56;

sx=[linspace(sx6k,sx6v,N), linspace(sx7k,sx7v,N)];
sy=[linspace(sy6k,sy6v,N), linspace(sy7k,sy7v,N)];

tolerancia=1;
for i=1:length(sx)
x_null=sx(i);
y_null=sy(i);
egyenlet=@(t,xy)[xy(1)*(a1-b1*xy(1)+c1*xy(2));xy(2)*(a2-b2*xy(2)+c2*xy(1))];
[t, xy]=ode45(egyenlet, [t_null, t_max], [x_null, y_null]);

%az szerint színezem a pályákat, hogy melyik egyensúlyi ponthoz vannak közel
% aztán visszafelé is kirajzolom a trajektóriákat: 
% a trajektória visszafelé: -f(x) 
if abs(xy(end,1)-3)<tolerancia
    egyenlet=@(t,xy)-[xy(1)*(a1-b1*xy(1)+c1*xy(2));xy(2)*(a2-b2*xy(2)+c2*xy(1))];
    [t, xyv]=ode45(egyenlet, [t_null, t_max], [x_null, y_null]);
    plot(xy(:,1), xy(:,2), 'c', xyv(:,1), xyv(:,2), 'c');
elseif abs(xy(end,2)-2)<tolerancia
    egyenlet=@(t,xy)-[xy(1)*(a1-b1*xy(1)+c1*xy(2));xy(2)*(a2-b2*xy(2)+c2*xy(1))];
    [t, xyv]=ode45(egyenlet, [t_null, t_max], [x_null, y_null]);
    plot(xy(:,1), xy(:,2), 'm', xyv(:,1), xyv(:,2), 'm');
end
end
%}

xlabel('x');
ylabel('y');
legend('vektormezõ', 'dx=0, mert x=0', 'dx=0, mert a1-b1*x+c1*y=0', 'dy=0, mert y=0', 'dy=0, mert a2-b2*x+c2*y=0', 'trajektoria');
xlim([-0.1, 4.1]);
ylim([-0.1, 4.1]);





