Homework Help: Modelling the spread of disease (in MatLab, Runge-Kutta)

1. Sep 4, 2011

Corrigan

1. The problem statement, all variables and given/known data

Hi!
I'm to build a model which solves the basic dynamic eqns in the so called SIR model in epidemology. I'm using MatLab and RK4.

There are tree populations in this model: Succeptible S, infected I and recovered R.

2. Relevant equations

The eqations are

dS/dt = -beta * S*I/N
dI/dt = beta * S*I/N - alpha*I
dR/dt = alpha*I

where alpha and beta are rates of recovery and infection (constants).

Vaccination will be introduced aswell, but I have set that to 0 for now.

3. The attempt at a solution

My base .m-file is this:

Code (Text):
clear all
close all

%Step size and number of iterations are defined
dt=1;
N=1000;         %days

P = 1000000;

%The initial conditions are entered as the first row in the array z.

z(1)=P-1;             % = S
z(2)=1;               % = I
z(3)=1;               % = R
t(1)=0;

%RK4 on function

for n=1:(N)

k1=P1f(t(n),z(n,:));
k2=P1f(t(n)+dt/2,z(n,:)+1/2*k1*dt);
k3=P1f(t(n)+dt/2,z(n,:)+1/2*k2*dt);
k4=P1f(t(n)+dt,z(n,:)+k3*dt);
z(n+1,:)=z(n,:)+1/6*(k1+2*k2+2*k3+k4)*dt;
t(n+1)=t(n)+dt;

end

%The resluts are plotted

subplot(1,2,2)
plot(t,z(:,2),'r')
axis([0 1000 0 100000])
xlabel(''),ylabel('')
subplot(1,2,1)
plot(t,z(:,3))
axis([0 1000 0 100000])
xlabel(''),ylabel('')
with the function P1f.m:

Code (Text):
function f=P1f(t,z)

%Def constants
alpha = 0.25;
sigma = 50;
T = 0.1;
beta = sigma*T;
gamma= 0;

f(1) = -beta*z(1)*z(2)/z(3) - gamma*z(1);
f(2) = beta*z(1)*z(2)/z(3) - alpha*z(2);
f(3) = gamma*z(1) + alpha*z(2);
end
The constants k in the RK4 routine all become NaN, and I don't understand why. This happenes during the first iteration.
Most of my problem is that I don't fully understand how the 'function' definition in matlab works. As you can see, I only define the RHS of the original eqations for the first 3 elements in z. How does the function calculate the values past that first row?
I know the method is sound because I've used it before, last year, and I've forgotten how it works.

Very grateful for any help,
Eric

EDIT: Major mistake was dividing by R instead of N = R + I + S in the right hand sides. Seems to be working now, as far as I can see.

Last edited: Sep 4, 2011