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

[Corrigan]

Homework Statement

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.

Homework 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.

The Attempt at a Solution

My base .m-file is this:

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:

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.

 

[KataKoniK]

Hi Eric,

It seems like you have a good start on your model. To answer your question about how the function calculates the values past the first row, the function is essentially a set of instructions that tell MATLAB what to do with the input values. In this case, the input values are the time (t) and the array of populations (z). The function uses these values to calculate the derivatives of the populations (dS/dt, dI/dt, dR/dt) and returns them in the array f. The RK4 routine then uses these derivatives to update the populations in the next time step.

As for your issue with the constants becoming NaN, it seems like you have already found the mistake and corrected it. It's always a good idea to double check your equations and make sure all the variables are being used correctly. Good luck with your model! Let us know if you have any other questions.