%% 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.