Programming for Dynamical System : SIR Model by Octave

% Function file
% the rhs of the SIR system. Saved as sirrhs.m
function xp = sirrhs(t,x)
global a b;
xp= zeros(2,1);
xp(1) = -a*x(1)*x(2) ;
xp(2) = a*x(1)*x(2)-b*x(2) ;
endfunction

% Call File
% SIR model with constant population. Recovered get immunity.
% Susceptibles, S(t)=x(1)
% Infectives, I(t)=x(2)
% Recovered, R(t)
% dS/dt=-a*S*I, dI/dt=a*S*I-b*I, dR/dt=b*I
% Use that the total population N=S+I+R to solve only first two equations
global a b;
a=.01;
b=.1;
%initial conditions
x0=[99;1];
%Total population
N=100;

% time interval
t0 = 0;
tf =60;
% plot the 2-d vector field with arrows.

figure(1)
clf;
[S,I] = meshgrid(0:10:100);
%xp=sirrhs(1,[S,I]);
Sq = -a*S.*I;
Iq = a*S.*I-b*I;
quiver(S,I,Sq,Iq,1);
%print -depsc2 sir_fig1.eps;

[t x] = ode45('sirrhs',[t0 tf],x0);
%Calculate recovered
R= N-x(:,1)-x(:,2);

figure(2)
clf
plot(t,x(:,1),'b-','LineWidth',2,t,x(:,2),'r-','LineWidth',2,t,R,'g-','LineWidth',2);
ylabel('Population Size','fontsize',15);
xlabel('Time','fontsize',15);
legend('Susceptibles', 'Infectives', 'Recovered');
%print -depsc2 sir_fig2.eps;

figure(3)
clf
quiver(S,I,Sq,Iq,1);
hold on;
plot(x(:,1),x(:,2),'b-','LineWidth',2);
ylabel('Infectives','fontsize',15);
xlabel('Susceptibles','fontsize',15);
%legend('Susceptibles', 'Infectives', 'Recovered');
hold off;
%print -depsc2 sir_fig3.eps;
#######################################
#######################################
For theoretical background and output see: www.ssbtr.net --> Publication ---> Programs
----
Durjoy Majumder, Ph.D.
West Bengal State University 
Secretary, SSBTR

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