%gyak 5 Lotka-Volterra
%{
% 1. egyenletek: dx=x*(1-2*x-y); dy=y*(3-2*y-x);
max_pont=4.; %a vizsgált pontok megadás
min_pont=0;
lepes=0.1;
x=min_pont:lepes:max_pont;
y=min_pont:lepes:max_pont;
figure;
[X,Y]=meshgrid(x,y); %elkészíti a hálót a rácspontok kiszámításához
%
dx=X.*(1-2*X-Y);
dy=Y.*(3-2*Y-X);
dx=dx./ sqrt(dx.^2+dy.^2); %vektormezõ lenormálva
dy=dy./sqrt(dx.^2+dy.^2);
quiver(X,Y,dx, dy, 'k');
%
%quiver(X,Y,X.*(1-2*X-Y),Y.*(3-2*Y-X), 20, 'k'); 
% (a 2:=vektorok hosszára vonatkozik, az alapértelmezett hossz 2-szerese)
% az x y rácspontokon a 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,
hold on; %az eddigi ábrára rajzolok a továbbiakban is
y0=1-2*x; %dx egyenese 1-2*x-y=0
x0=3-2*y; %dy egyenese 3-2*y-x=0
plot(x, y0,'b', zeros(1, length(x)), y, 'b', x0, y, 'r', x, zeros(1, length(y)), 'r'); 
% kirajzol: dx=0 ha 1-2*x-y=0; x=0 kékkel; dy=0 ha 3-2*y-x=0; y=0 pirossal
%ábrazolás:
xlim([-0.1, max_pont]);
ylim([-0.1, max_pont]);
title('egyenletek: dx=x*(1-2*x-y); dy=y*(3-2*y-x)');
xlabel('x');
ylabel('y');
legend(' vektormezõ gradiense adott pontokban', 'dx=0, mert x=0', 'dx=0, mert 1-2*x-y=0', 'dy=0, mert y=0', 'dy=0, mert 3-2*y-x=0');
hold off;
%}

% A különbözõ paraméterek változtatása mellett, vizsgálható, hol alakulnak
% ki egyensúlyi pontok, és, hogy a diffrenciál egyenlet által leírt
% vektormezõ, milyen irányú és méretû lesz.
%{
% 2.A egyenletek: dx=b*x*(1-2*x)-x*y; dy=y*(3-2*y-x);
figure;
max_pont=4;
min_pont=0;
lepes=0.1;
for b=0.1:lepes:max_pont %ciklus, amiben a paramétert változtatom
x=min_pont:lepes:max_pont;
y=min_pont:lepes:max_pont;
[X,Y]=meshgrid(x,y);
%
dx=b*X.*(1-2*X)-X.*Y; %dx=X.*(b-b*2*X-Y);
dy=Y.*(3-2*Y-X);
dx=dx./ sqrt(dx.^2+dy.^2);
dy=dy./sqrt(dx.^2+dy.^2);
quiver(X,Y,dx, dy, 0.5, 'k');
%
%quiver(X,Y,X.*(b-b*2*X-Y),Y.*(3-2*Y-X), 2, 'k'); %irányokba vektorok (irányt 2 koord-dal ad meg)
y0=b*(1-2*x);
x0=3-2*y;
hold on;
plot(x, y0,'b', zeros(1, length(x)), y, 'b', x0, y, 'r', x, zeros(1, length(y)), 'r');
xlim([-0.1, max_pont]);
ylim([-0.1, max_pont]);
title(['egyenletek: dx=b*x*(1-2*x)-x*y; dy=y*(3-2*y-x); b=', num2str(b)]);
xlabel('x');
ylabel('y');
legend(' vektormezõ gradiense adott pontokban', 'dx=0, mert x=0', 'dx=0, mert b-b*2*x-y=0', 'dy=0, mert y=0', 'dy=0, mert 3-2*y-x=0');
hold off;
%pause(0.1);
waitforbuttonpress; %gombnyomásra lép egyet a szimuláció
end
%}

%{
% 2.A.2 egyenletek: dx=x*(b-2*x)-x*y; dy=y*(3-2*y-x);
figure;
max_pont=8;
min_pont=0;
lepes=0.1;
for b=0.1:lepes:max_pont %ciklus, amiben a paramétert változtatom
x=min_pont:lepes:max_pont;
y=min_pont:lepes:max_pont;
[X,Y]=meshgrid(x,y);
%
dx=X.*(b-2*X)-X.*Y; %dx=X.*(b-2*X-Y);
dy=Y.*(3-2*Y-X);
dx=dx./ sqrt(dx.^2+dy.^2);
dy=dy./sqrt(dx.^2+dy.^2);
quiver(X,Y,dx, dy, 0.5, 'k');
%
%quiver(X,Y,X.*(b-2*X-Y),Y.*(3-2*Y-X), 2, 'k'); %irányokba vektorok (irányt 2 koord-dal ad meg)
y0=(b-2*x);
x0=3-2*y;
hold on;
plot(x, y0,'b', zeros(1, length(x)), y, 'b', x0, y, 'r', x, zeros(1, length(y)), 'r');
xlim([-0.1, 4]);
ylim([-0.1, 4]);
title(['egyenletek: dx=x*(b-2*x)-x*y; dy=y*(3-2*y-x); b=', num2str(b)]);
xlabel('x');
ylabel('y');
legend(' vektormezõ gradiense adott pontokban', 'dx=0, mert x=0', 'dx=0, mert b-2*x-y=0', 'dy=0, mert y=0', 'dy=0, mert 3-2*y-x=0');
hold off;
%pause(0.1);
waitforbuttonpress; %gombnyomásra lép egyet a szimuláció
end
%}

%
% 2B. egyenletek: dx=x*(1-x/k-y); dy=y*(3-2*y-x);
figure;
max_pont=5;
min_pont=0.;
lepes=0.1;
for k=0.1:0.1:max_pont
x=min_pont:lepes:max_pont;
y=min_pont:lepes:max_pont;
[X,Y]=meshgrid(x,y);
%
dx=X.*(1-X/k-Y);
dy=Y.*(3-2*Y-X);
dx=dx./ sqrt(dx.^2+dy.^2);
dy=dy./sqrt(dx.^2+dy.^2);
quiver(X,Y, dx, dy, 'k');
%
%quiver(X,Y,X.*(1-X/k-Y),Y.*(3-2*Y-X), 2, 'k'); %irányokba vektorok (irányt 2 koord-dal ad meg)
y0=1-x/k;
x0=3-2*y;
hold on;
plot(x, y0,'b', zeros(1, length(x)), y, 'b', x0, y, 'r', x, zeros(1, length(y)), 'r');
xlim([-0.1, max_pont]);
ylim([-0.1, max_pont]);
title(['egyenletek: dx=x*(1-x/k-y); dy=y*(3-2*y-x); k=', num2str(k)]);
xlabel('x');
ylabel('y');
legend(' vektormezõ gradiense adott pontokban', 'dx=0, mert x=0', 'dx=0, mert 1-x/k-y=0', 'dy=0, mert y=0', 'dy=0, mert 3-2*y-x=0');
hold off;
waitforbuttonpress;
end
%}

%{
% 2C. egyenletek: dx=x*(1-2*x-y*g); dy=y*(3-2*y-x);
figure;
max_pont=4;
min_pont=0;
lepes=0.1;
for grec=0.1 :lepes:max_pont
    g=1/grec;
x=min_pont:lepes:max_pont;
y=min_pont:lepes:max_pont;
[X,Y]=meshgrid(x,y);
%
dx=X.*(1-2*X-Y*g);
dy=Y.*(3-2*Y-X);
dx=dx./ sqrt(dx.^2+dy.^2);
dy=dy./sqrt(dx.^2+dy.^2);
quiver(X,Y, dx, dy, 'k');
%
%quiver(X,Y,X.*(1-2*X-Y*g),Y.*(3-2*Y-X), 2, 'k'); %irányokba vektorok (irányt 2 koord-dal ad meg)
y0=(1-2*x)/g;
x0=3-2*y;
hold on;
plot(x, y0,'b', zeros(1, length(x)), y, 'b', x0, y, 'r', x, zeros(1, length(y)), 'r');
xlim([-0.1, max_pont]);
ylim([-0.1, max_pont]);
title(['egyenletek: dx=x*(1-2*x-y*g); dy=y*(3-2*y-x); 1/g=', num2str(grec)]);
xlabel('x');
ylabel('y');
legend(' vektormezõ gradiense adott pontokban', 'dx=0, mert x=0', 'dx=0, mert 1-2*x-y*g=0', 'dy=0, mert y=0', 'dy=0, mert 3-2*y-x=0');
hold off;
waitforbuttonpress;
end
%}

%{
% 3. egyenletek: dx=x*(1-x/k-y*g); dy=y*(3-2*y-x);
figure;
max_pont=5;
min_pont=0.;
lepes=0.1;
%paraméterek: k és g, esetek: k<3 g<2/3; k<3 g>2/3; k>3 g<2/3; k>3 g>2/3
k=[2 4];
g=[1 0.5];
for i=1:length(k)
    for j=1:length(g)
x=min_pont:lepes:max_pont;
y=min_pont:lepes:max_pont;
[X,Y]=meshgrid(x,y);
%
dx=X.*(1-X/k(i)-Y*g(j));
dy=Y.*(3-2*Y-X);
dx=dx./ sqrt(dx.^2+dy.^2);
dy=dy./sqrt(dx.^2+dy.^2);
quiver(X,Y, dx, dy, 'k');
%
%quiver(X,Y,X.*(1-X/k-Y*g),Y.*(3-2*Y-X), 2, 'k'); %irányokba vektorok (irányt 2 koord-dal ad meg)
y0=(1-x/k(i))/g(j);
x0=3-2*y;
hold on;
plot(x, y0,'b', zeros(1, length(x)), y, 'b', x0, y, 'r', x, zeros(1, length(y)), 'r');
xlim([-0.1, max_pont]);
ylim([-0.1, max_pont]);
title(['egyenletek: dx=x*(1-x/k-y*g); dy=y*(3-2*y-x); k=', num2str(k(i)), ', g=', num2str(g(j))]);
xlabel('x');
ylabel('y');
legend(' vektormezõ gradiense adott pontokban', 'dx=0, mert x=0', 'dx=0, mert 1-x/k-y*g=0', 'dy=0, mert y=0', 'dy=0, mert 3-2*y-x=0');
hold off;
waitforbuttonpress;
    end
end
%}

%paraméterek és sajátértékek kapcsolata:
%{
%A.)
% S: 1, 2 %melyik ponthoz melyik egyenletek kellenek (egyenket, jobboldal)
% Q: 2 3
% P: 1 4
B=(0.5:0.01:3); %végignézett b értékek
SE=zeros(2, size(B, 2)); %sajátértékeket itt tárolom
e1=1; %melyik egyenletek nézem?
e2=2;
for d=1:size(B, 2)
b=B(1,d); %aktuális b érték
egyenlet=[-2*b -1; -1 -2; 1 0; 0 1]; %egyenletek
jobboldal=[-b; -3; 0; 0]; %jobb oldalaik
A=[egyenlet(e1,:); egyenlet(e2,:)]; %egyemlet
v=[jobboldal(e1); jobboldal(e2)]; %egyenlet jobb oldala
mg=A\v; %egyenlet megoldása
x=mg(1); %adott b mellett aktuális x 
y=mg(2); %és y
J=[(b-4*x*b-y) (-x); (-y) (3-4*y-x)]; %Jacobi mátrix
sajatertek=eig(J); %J sajátértékei
SE(:,d)=sajatertek;
end
figure; %ábrázolás
plot(B, SE(1,:), 'b.', B, SE(2,:), 'r.', B, 0, 'g');
legend('sajátérték 1', 'sajátérték 2');
xlabel('b paraméter értékei');
ylabel('sajátértékek');
%}

%{
%B.)
% S: 1, 2
% Q: 2 3
% P: 1 4
K=(3:0.01:6); 
SE=zeros(2, size(K, 2));
e1=1;
e2=2;
for d=1:size(K, 2)
k=K(1,d);
egyenlet=[-1/k -1; -1 -2; 1 0; 0 1];
jobboldal=[-1; -3; 0; 0];
A=[egyenlet(e1,:); egyenlet(e2,:)]; %egyemlet
v=[jobboldal(e1); jobboldal(e2)]; %egyenlet jobb oldala
mg=A\v; %egyenlet megoldása
x=mg(1); %adott b mellett aktuális x 
y=mg(2); %és y
J=[(1-2*x/k-y) (-x); (-y) (3-4*y-x)]; %Jacobi mátrix
sajatertek=eig(J); %J sajátértékei
SE(:,d)=sajatertek;
end
figure;
plot(K, SE(1,:), 'b.', K, SE(2,:), 'r.', K, 0, 'g');
legend('sajátérték 1', 'sajátérték 2');
xlabel('K paraméter értékei');
ylabel('sajátértékek');
%}

%{
%C.)
% S: 1, 2
% Q: 2 3
% P: 1 4
G=(0:0.01:(2/3)); 
SE=zeros(2, size(G, 2));
e1=2;
e2=1;
for d=1:size(G, 2)
g=G(1,d);
egyenlet=[-2 -g; -1 -2; 1 0; 0 1];
jobboldal=[-1; -3; 0; 0];
A=[egyenlet(e1,:); egyenlet(e2,:)]; %egyemlet
v=[jobboldal(e1); jobboldal(e2)]; %egyenlet jobb oldala
mg=A\v; %egyenlet megoldása
x=mg(1); %adott b mellett aktuális x 
y=mg(2); %és y
J=[(1-4*x-y*g) (-x*g); (-y) (3-4*y-x)]; %Jacobi mátrix
sajatertek=eig(J); %J sajátértékei
SE(:,d)=sajatertek;
end
imag(SE)
figure;
plot(G, SE(1,:), 'b.', G, SE(2,:), 'r.', G, 0, 'g');
legend('sajátérték 1', 'sajátérték 2');
xlabel('gamma paraméter értékei');
ylabel('sajátértékek');
%}