%7. gyakorlat: 
%
%1.dinamikai rendszer, ha az egyenlet rajzol:
% 2B. egyenletek: dx=x*(1-x/k-y); dy=y*(3-2*y-x);
figure;
max_pontx=5;
max_ponty=2;
min_pont=0.;
lepes=0.1;
x=min_pont:lepes:max_pontx;
y=min_pont:lepes:max_ponty;
[X,Y]=meshgrid(x,y);

k=5;
y0=1-x/k;
x0=3-2*y;
hold on;
plot(x, y0,'b', zeros(1, length(y)), y, 'b', x0, y, 'r', x, zeros(1, length(x)), 'r');
xlim([-0.1, 6]);
ylim([-0.1, 6]);
title(['egyenletek: dx=x*(1-x/k-y); dy=y*(3-2*y-x); k=', num2str(k)]);
xlabel('x');
ylabel('y');
legend('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;

% rács pontjaiból indulok:
%sx=reshape(X, 1, length(x)*length(y)); %kezdeti pontok
%sy=reshape(Y, 1, length(x)*length(y)); %kezdeti pontok

%
N=100; %10;
%{
%egyenesen kiindulási pontok
sxk=0.8;%0;%
syk=0.3;%0.5;%
sxv=1.2;%2;%
syv=0.15;%0;%

sx2k=2;%0;%
sy2k=1.35;%2;%
sx2v=2.3;%6;%
sy2v=1.2;%0;%

%sx=linspace(sxk,sxv,N);
%sy=linspace(syk,syv,N);

sx=[linspace(sxk,sxv,N), linspace(sx2k,sx2v,N)];
sy=[linspace(syk,syv,N), linspace(sy2k,sy2v,N)];
%}
%{
%idõben vissza: dx=-dx
sx3k=4.5;%
sy3k=0.01;%0;%
sx3v=4.5;%
sy3v=0.02;%0.1;%

sx4k=1;%
sy4k=0.945;%0;%
sx4v=1;%
sy4v=0.965;%2;%

%sx=linspace(sx3k,sx3v,N);
%sy=linspace(sy3k,sy3v,N);

sx=[linspace(sx3k,sx3v,N), linspace(sx4k,sx4v,N)];
sy=[linspace(sy3k,sy3v,N), linspace(sy4k,sy4v,N)];
%}

%körben kiiindulási pontok:
%{
A=[1 2; 1/k 1];
    b=[3;1];
    mg3=A\b;
    mg1=[0;1.5];
    mg2=[5;0];
    mg0=[0;0];
   
    Nl=100; %pontok száma körvonalon
    phi=linspace(0,2*pi,Nl);
    r=0.5; %Kör sugara
    
     kor=repmat(mg3,1,length(phi))+[r*sin(phi); r*cos(phi)];
     %A populáció mérete nem lehet negatív, ezért ezeket ki kell zárnunk
     kor_jo=kor(:,kor(1,:)>=0 & kor(2,:)>=0);
     sx=kor_jo(1,:);
     sy=kor_jo(2,:);
%}
%
%figure
t0=[0 10];
for i=1:length(sx)
    x0=[sx(i); sy(i)];
    [t, y] = ode45(@(t, y) [y(1)*(1-y(1)/k-y(2)); y(2)*(3-2*y(2)-y(1))],t0,x0);
    %[t, y] = ode45(@(t, y) -[y(1)*(1-y(1)/k-y(2)); y(2)*(3-2*y(2)-y(1))],t0,x0); %az idõ visszafelé
    plot(y(:,1),y(:,2), 'g');
    hold on   
end
%}

%{
%2. biz pontokból trajekt rajzol,
%a mátrixot változtatom, hogy változik a fázisportré (sé, det, trace jelölve)
x = -2:1:2;
y = -2:1:2;
[X Y] = meshgrid(x,y);

K = 500;
fix = zeros(K,length(x)*length(y));
fiy = zeros(K,length(x)*length(y));

figure;
%trace, det görbe kirajzol
subplot(1,2,2);
tr=(-8:0.1:8);
plot(tr, tr.^2./4, 'b', tr, zeros(size(tr)), 'b', zeros(size(tr)), tr, 'b');
xlabel('trace');
ylabel('det');
hold on;
%
% paraméter változtat
hatar=3;
for a = -hatar:0.2:hatar 
    A = [0 1; -4+a^2 2*a] %egyenletrendszer mátrixos alakban
   
se=eig(A);
%adott rendszerre difegyet megoldom:
 n=1;
for i = 1:length(x)
    for j = 1:length(y)
        %
        if trace(A)>0 
           A1=-A; %idõ hátra (trajektória végpontjai a kijelölt pontok) ha elmennek a megoldások
        else
            A1=A;
        end
        %
        [t,fi] = ode45(@(t,y) A1*y, linspace(0,5,K), [x(i) y(j)]);
        fix(:,n) = fi(:,1);
        fiy(:,n) = fi(:,2);
        n = n+1;
    end
end

%plotol trajektóriákat adott pontokból:
subplot(1,2,1);
plot(fix,fiy);
hold on;
plot(X,Y,'x');
xlim([2*min(x) 2*max(x)]);
ylim([2*min(y) 2*max(y)]);
title(['sajátértékek: ', num2str(se(1)), ', ', num2str(se(2))]);
xlabel('x');
ylabel('y');
hold off

%trace det kirajzol:
subplot(1,2,2);
plot(trace(A), det(A), 'r*');
title(['trace(A): ', num2str(trace(A)), ', det(A):', num2str(det(A))]);

waitforbuttonpress;
%pause(0.1)
end
%}

%3.
% Hopf-bifurkáció:
%{
% Idõ változtatása adott mu esetén, folyamatosan rajzol
mu =0.6;%1.8
x0= [0 1]; 
N=40;
    t0=[0 N];
    [t, y] = ode45(@(t,y) [1/8 - y(1)+y(1)^2*y(2); mu - y(1)^2*y(2)],t0,x0); %az egyenletben a t és az y a változók a mu fix paraméter
 for i=1:length(t) %eddig rajzolom ki a folyamatot
    subplot(2, 2, 1);%x és y akakulása egymás függvényében
    plot(y(1:i,1),y(1:i,2))
    title('dx/dt=1/8-x+x^2*y dy/dt=mu - x^2*y');
    xlim([0 7]);
    xlabel('x');
    ylim([0 7]);
    ylabel('y');
    
    subplot(2, 2, 2);%y irányú elmozdulás az idõben
    plot(t(1:i),y(1:i,2))
    xlim([0 N]);
    xlabel('t');
    ylim([0 7]);
    ylabel('y');
    
    subplot(2, 2, 3); %x irányú elmozdulás az idõben
    plot(y(1:i,1),t(1:i))
    ylim([0 N]);
    ylabel('t');
    xlim([0 7]);
    xlabel('x');
    
    pause(0.1);
end
%}

%{
% a mû változtatásával kialakuló bifurkáció (adott idõ)
t0 = [0 400];
x0= [0 1];
figure
for mu = 2:-0.05:0 
    mu;
[t, y] = ode45(@(t,y) [1/8 - y(1)+y(1)^2*y(2); mu - y(1)^2*y(2)],t0,x0);

subplot(2, 2, 1);
plot(y(:,1),y(:,2))
 title(['dx/dt=1/8-x+x^2*y dy/dt=mu - x^2*y, mu: ', num2str(mu)]);
    xlim([0 7]);
   xlabel('x');
    ylim([0 7]);
    ylabel('y');

subplot(2, 2, 2);
plot(t, y(:,2));  
    xlim([0 120]);
    xlabel('t');
    ylim([0 7]);
    ylabel('y');

subplot(2, 2, 3);
plot(y(:,1), t)   
    ylim([0 120]);
    ylabel('t');
    xlim([0 7]);
    xlabel('x');

pause(0.1);
%waitforbuttonpress
end
%}

%{
%% Ha kör mentén indítjuk a rendszert: az egyensúlyi pont körüli viselkedés
figure;
t0 = [0 400]; 
mu=0.65;
x00=[1.6, 0.5];
%x00=[1/8+mu mu/(1/8+mu)^2]; %egyensúlyi pont: mindét egyenletet kielégíti,
r0=10;%0.4; %kör sugara, ahonnan pontokat veszünk (ezen pontokból indítjuk majd a rendszert)

%{
%kör kirajzolása pirossal:
x=0:0.01:2*pi;
x0=repmat(x00',1,length(x))+r0*[cos(x);sin(x)];
plot(x0(1,:),x0(2,:),'r');
%}

 imax=70; %ennyi pont
for i=0:imax %sok pont egy körön
    x0=x00+r0*[cos(i*2*pi/imax),sin(i*2*pi/imax)];   %adott körön lévõ pont koordinátái
    [t, y] = ode45(@(t,y) [1/8 - y(1)+y(1)^2*y(2); mu - y(1)^2*y(2)],t0,x0);
    hold on;
    plot(y(:,1),y(:,2), x0(1), x0(2), 'bx');%  trajeltóriák és pontok ábrázolása
    xlim([-20 20]);
    ylim([-20 20]);
    pause(0.1);
end
%}

