% =========================================================================
% Biological_Models_Prey_Predator_Lotka_Volterra.m
%  This is a script-function that simulates the Lotka-Volterra pray-predator model.
% 
%  Pray - Predator Biological Models  
%
% This m-file describes a prey-predator systems.
% The general form of the system's equations is:
%
%  dN/dt = F(N,P) 
%  dP/dt = G(N,P)
%
% =========================================================================
% Different models can be considered with different functional forms for F
% and G. 
% Here we consider the following models:
% (a) Lotka-Volterra model
%       dN/dt = N(a-bP) 
%       dP/dt = P(cN-d)
%
% Model parameters are a, b, c, d
%
% -- First we draw a phase portrait for the pray-predator model. 
%    Results are shown at each time step. 
% -- Second, we plot the populations of pray and predator for the time
% period
%
% Equations are solved using a numerical ODE solver. 
%
% 04/15
% =========================================================================
 
function Lotka_Volterra()
 
% Clears history and command window 
close all; clear all;
 
% -- run parameters
pausetime = 0.1;  % Shows solutions at each time step. 
iterations = 1;  % Sets initial interation count to 1;
runtime = 20;    % Duration time of simulation.
 
% -- Equation parameter values
a = 1.0;
b = 2.0;
c = 2.0;
d = 1.0;
 
% -- Initial conditions for x and y
initialx = 1.0;
initialy = 1.0;
 
fprintf('----------------------------------\nLotkaVolterra Predetor Prey model \n\nMatlab code written by James Adams\n----------------------------------')
%fprintf('\n\nParameter values set,')
%fprintf('\n\nalpha = %2.6f \nbeta = %2.6f \ngamma = %2.6f \ndelta = %2.6f ',alpha,beta,gamma,delta)
%fprintf('\n\nTo plot a time series graph for x, set the choice parameter in the code to 1.') 
%fprintf('\nSet choice to 2 to plot a time series graph for y.\n\nSimulation will run in ')
 
% Solves equations using numerical ODE solver 45 (nonstiff runge kutta)
deq1=@(t,x) [x(1)*(a - b*x(2)); x(2)*(c*x(1)-d)];
[t,sol] = ode45(deq1,[0 runtime],[initialx initialy]);
 
arraysize = size(t);  % Sets time array size for the for loop.
 
% -- Solutions are plotted at each time step
for i = 1 : max(arraysize) 
    plot(sol(iterations,1),sol(iterations,2),'.','color',[rand; rand; rand],'markersize',14,'MarkerFaceColor','b');                               
    hold on
    title(['Lotka-Volterra Equations       t = ' num2str(t(iterations))],'fontsize',12)
    xlabel('N','fontsize',12)
    ylabel('P','fontsize',12)
    axis([min(sol(:,1)) max(sol(:,1)) min(sol(:,2)) max(sol(:,2))])
   
    %text(0.1,0.5,'Time Series graph will be shown at the end of the simulation')
    
    iterations = iterations + 1;   % Adds 1 to the iteration count. 
    pause(pausetime)
end
 
% -- Graph of x and y time series
figure (2);
hold on;
%title(['Lotka-Volterra time series for x and y,  run time = ' num2str(max(t)) ' seconds '],'fontsize',12)
%plot(t(:,1),sol(:,1),'b.','markersize',10)
plot(t(:,1),sol(:,1),'b')
plot(t(:,1),sol(:,2),'r','markersize',10)
%axis([0 max(t(:,1)) 0 max(sol(:,2),sol(:,1))])
%axis([min(t(:,1)) max(t(:,1)) min(sol(:,2)) max(sol(:,2))])
legend('Pray','Predator');
xlabel('Time Steps');
ylabel('Population');
hold off;
