Programming for Dynamical System : Lotka-Volterra Prey Predator by Octave

Euler Method

%forward Euler

K(1)=5.0; %initial krill population
a=1.0; %(a-bW) is the intrinsic growth rate of the krill
b=0.5;

W(1)=1.0;%initial whale population
c=0.75; %(-c+dK) is the intrinsic growth rate of whales
d=0.25;

tinit=0.0;
tfinal=50;
n=5000;%no. of time steps
dt=(tfinal-tinit)/n;%time step size
%T=[tinit:dt:tfinal]; %create vector of discrete solution times

%Execute forward Euler to solve at each time step
for i=1:n
    K(i+1)=K(i)+dt*K(i)*(a-b*W(i));
    W(i+1)=W(i)+dt*W(i)*(-c+d*K(i));
end;

%Plot Results...
%S1=sprintf('IC: W0=%g, K0=%g',W(1),K(1));

figure(1);
plot(K,W);
title('Phase Plane Plot');
xlabel('Krill');
ylabel('Whales');
%legend(S1,0)
%grid;
%init=1;

init=1;
figure(2);
%clg;
plot((init:tfinal),K(init:tfinal),'r-',(init:tfinal),W(init:tfinal),'g-.');
legend('Krill','Whales');
xlabel('time');
ylabel('whales and krill');

----- ##### ------ ##### ----- ##### ----- ##### ----- ##### ------ #####
----- ##### ------ ##### ----- ##### ----- ##### ----- ##### ------ #####

Runge-Kutta Method

% function file VolterraLotka
function res = VolterraLotka(x,t)
c1 = 1;  c2 = 2;  c3 = 1; c4 = 1;
res = [(c1-c2*x(2))*x(1); (c3*x(1)-c4)*x(2)];
endfunction




% VolterraLotkaCall File
% creating a vector field
x = 0:0.2:2.6; %define the x values to be examined
y=0:0.2:2.0; %define the y values to be examined
n=length(x); m=length(y);
Vx=zeros(n,m); Vy=Vx; %create zero vectors for vector field

for i=1:n
    for j=1:m
        v=VolterraLotka([x(i),y(j)],0); %compute the vector
        Vx(i,j)=v(1); Vy(i,j) =v(2);
        endfor
        endfor
       
        figure(1);
        quiver(x,y,Vx',Vy',3);
        axis([min(x),max(x),min(y),max(y)]);
        grid on; xlabel('prey');ylabel('predator');
       
% with the function VolterraLotka we can solve the differential eqn.
% for times 0<= t <=15 with the
% initial conditions x(0)=2 and y(0)=1.
%In the code below we choose to use 100 different values for
% the time t to display the solution. Observe the algorithm lsode() internally uses % considerably more points.

       
t= linspace(0,15,100);
XY=lsode('VolterraLotka',[2,1],t); %COMMAND
       
figure(2);
plot(XY(:,1),XY(:,2));
axis([min(x),max(x),min(y),max(y)]); hold on
quiver(x,y,Vx',2);
grid on; xlabel('prey');ylabel('predator'); hold off
               
figure(3);
plot(t,XY)
xlabel('time'); legend('prey','predator'); axis([0,15,0,3]); grid on;

For Background : See www.ssbtr.net ---> Publication ---> Program

-----
Durjoy Majumder, Ph.D.
West Bengal State University
Secretary, SSBTR