%% Interaktyvus pratimas nr. 6. Modelis "plėšrūnas-auka"
% Kitaip  žinomas kaip Lotkos-Volteros modelis šiš modelis aprašo dviejų gyvūnų 
% populiacijų (plešrūnų ir žoliaėdžių) tarpusavio sąveiką natūralioje aplinkoje. 
% Dažniausiai modelio aprašyme kaip plėsrūnai figūruoja lapės, o kaip žoliaėdžiai 
% - kiškiai (arba triūšiai, zuikiai). Daromos prielaidos, kad kiškių maisto resursai 
% (žolė, medžio žievė) yra neriboti, o lapių maistas yra kiškiai. Tokią ekosistemą 
% galima aprašyti diferencialinių lygčių sistema:
% 
% $$\left\lbrace \begin{array}{ll}u\prime \;\left(t\right)= & \left(a-b\cdot 
% v\left(t\right)\right)u\left(t\right)\\v\prime \left(t\right)= & \left(-c+d\cdot 
% u\left(t\right)\right)v\left(t\right)\end{array}\right.,$$ 
% 
% čia $u\left(t\right)$ yra kiškių populiacijos funkcija nuo laiko $t$, o $v\left(t\right)$ 
% yra lapių populiacijos funkcija. Koeficientai šioje sistemoje yra suprantamos 
% kaip tikimybės:
%% 
% * $a$ - tikimybė kad kiškių populiacija didės (kiškiai dauginasi);
% * $b$ - tikimybė kad kiškių populiacija mažės (lapės valgo kiškius);
% * $c$ - tikimybė kad lapių populiacija mažės (lapės miršta nuo bado);
% * $d$ - tikimybė kad lapių populiacija didės (lapės turi užtektinai maisto, 
% lapės dauginasi).
%% 
% Akivaizdu, kad šie koeficientai priklauso intervalui $\left\lbrack 0;1\right\rbrack$. 
% Būtent šiuos koeficientus gali keisti šio kurso klausytojai interaktyviame pratime. 
% Taip pat galite keisti populiacijų dydžius pradiniu momentu ir laiko intervalo 
% galą $T$, t.y. stebimas laiko intervalas yra $\left\lbrack 0;T\right\rbrack$. 
% Pirmajame grafike galite stebėti kaip keičiasi abiejų populiacijų dydžiai pasirinktame 
% laiko intervale.
% 
% Modelyje taip pat nėra nagriinėjama gyvūnų migracija, t.y. nei kiškiai, nei 
% lapės nemigruoja iš nagrinėjamos ekosistemos. Taip pat į ekosistemą nemigruoja 
% jokie kiti givūnai. Susidomėję šią sistemą skaitytojai gali patys surasti kitų 
% modelių, kuriuose šie veiksniai gali būti įtraukiami į nagrinėjamą ekosistemą.

a = 0.56;
b = 0.28;
c = 0.75;
d = 0.25;

t_end = 100;

t_int = [0 t_end];

rabbits = 2;
foxes = 3;

IC = [rabbits foxes];
%%
[t, X] = ode45(@(t, X) preditor_prey(t, X, a, b, c, d), t_int, IC);

Rabbits = X(:, 1);
Foxes = X(:, 2);
%%
plot(t, Rabbits);
hold on;
plot(t, Foxes);
legend('Kiškiai', 'Lapės')
xlabel('Laikas t')
ylabel('Populiacijos dydis')
hold off;
%% 
% Reikia paminėti, kad ši sistema turi ypatinguosius taškus:koordinačių pradžios 
% taškas ir taškas $u\left(0\right)=\frac{a}{b}$ , $v\left(0\right)=\frac{c}{d}$. 
% Aišku, kad pirmasis taškas yra neįdomus, nes aprašo ekosistemą, kurioje nėra 
% gyvūnų. Jeigu nagrinėjamas kitas ypatingasis taškas, tai ekosistema yra balansuota: 
% kiškių ir lapių populiacijos yra pastovios arba labai mažai kinta. Tą galite 
% matyti ir patys tinkamai parinkę parametrus. Jeigu pradžios taškas nėra ypatingasis, 
% tai galime matyti kaip keičiasi populiacijų fazinės kreivių pavidalas. Šias 
% kreives pavaizduosime antrame grafike, kuriame kiškių ir lapių populiacijų dydžiai 
% yra atitinkamai ašyse Ox ir Oy. Taip galime stebėti apiejų populiacijų dydžių 
% pokyčius.

plot(Rabbits, Foxes)
xlabel('Kiškių populiacijos dydis')
ylabel('Lapių populiacijos dydis')
%%
function dxdt = preditor_prey(t, X, a, b, c, d)
    dx1 = (a - b*X(2))*X(1);
    dx2 = (-c + d*X(1))*X(2);
    dxdt = [dx1; dx2];
end
%% Pabandykite patys
%% 
% # Remiantis pateiktu kodu pabandykite realizuoti Gudvino modelį, kurį taip 
% pat nagrinėjome paskaitoje.