(adsbygoogle = window.adsbygoogle || []).push({}); 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:

with the function P1f.m: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('')

The constants k in the RK4 routine all become NaN, and I don't understand why. This happenes during the first iteration.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

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.

**Physics Forums | Science Articles, Homework Help, Discussion**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

Can you offer guidance or do you also need help?

Draft saved
Draft deleted

**Physics Forums | Science Articles, Homework Help, Discussion**