%%  Wolf-Hase: Räuber-Beute-Modell
% System zweier gekoppelter ODES.
% Siehe auch die <https://de.wikipedia.org/wiki/Lotka-Volterra-Gleichungen Lotka-Volterra-Gleichungen>.

%%  Parameter der Evolutionsgleichungen festlegen
clear; clc; clf; close all
alpha = 8.5;    % Todesrate der Woelfe
beta  = 10;     % Vermehrungsrate der Hasen
gamma = 0.01;   % Nachwuchsrate der Woelfe (dank Beute=Hasen)
delta = 0.07;   % Todesrate der Hasen (durch die Jaeger=Woelfe)

time_a = 0;     % Startzeit
time_b = 2;     % Endzeit
tau = 1/1000;    % Schrittweite in der Zeit

t  = time_a:tau:time_b;
nt = length(t); % Anzahl der Zeitschritte

%% Anfangspopulationen
wolf(1) = 100;
hase(1) = 1000;

%% Zeititeration.
% <http://www.cae-wiki.info/wikiplus/index.php/Explizit-implizit_anschaulich Exlizites Schema> der Zeitintegration. 
% Erfordert sehr kleine Zeitschrittweite für numerische Stabilität.
wolf(nt) = 0;         % Stellt Speicher fuer alle benoetigetn Vektorelemente bereit.
hase(nt) = 0;

for m = 2:nt
    wolf(m) = wolf(m-1) -tau*alpha*wolf(m-1) + tau*gamma*wolf(m-1)*hase(m-1);
    hase(m) = hase(m-1) +tau*beta *hase(m-1) - tau*delta*wolf(m-1)*hase(m-1);
end

%% Ergebnisse anzeigen

% alt: mit plotyy
% [AX,H1,H2] = plotyy(t,wolf,t,hase);
% title('Räuber-Beute-Modell');
% xlabel('Zeit');
% set(H1,'LineStyle','--')
% set(H2,'Color','r')
% set(get(AX(1),'Ylabel'),'String','Wölfe') 
% set(get(AX(2),'Ylabel'),'String','Hasen') 
% set(get(AX(2),'Ylabel'),'Color','r') 

% neu: mit yyaxis
yyaxis left
plot(t,wolf,'--b')
ylabel('Wölfe')
hold on

yyaxis right
plot(t,hase,'-r')
ylabel('Hasen')

title('Räuber-Beute-Modell: simple explizite Zeitintegration');
xlabel('Zeit');

saveas(gcf,'wolf_hase.jpg');    % Speicher die aktuelle Grafik als jpg-File
hold off

% publish('wolf_hase.m')

%% Solve as system of ODEs

figure
tspan = [time_a time_b];
y0    = [wolf(1),hase(1)];
[t,yode] = ode45(@(t,y) odefcn(t,y,alpha,beta,gamma,delta), tspan, y0);

%% Visualisierung
% plot(t,yode)
% neu: mit yyaxis
yyaxis left
plot(t,yode(:,1),'--b')
ylabel('Wölfe')
hold on

yyaxis right
plot(t,yode(:,2),'-r')
ylabel('Hasen')

title('Räuber-Beute-Modell: Löser ode45');
xlabel('Zeit');

%% Beschreibung des ODE-Systems
% $$ \frac{dw(t)}{dt} = -\alpha w(t) + \gamma w(t)h(t) \qquad\qquad
%    \frac{dh(t)}{dt} =  \beta  w(t) - \delta w(t)h(t) $$

function dy = odefcn(t,y,alpha,beta,gamma,delta)
dy = zeros(2,1);
dy(1) = -alpha*y(1) + gamma*y(1)*y(2);
dy(2) =  beta *y(2) - delta*y(1)*y(2);
end

