1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
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)

  1. Sep 4, 2011 #1
    1. The problem statement, all variables and given/known data

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

    %RK4 on function

    for n=1:(N)

    %The resluts are plotted

    axis([0 1000 0 100000])
    axis([0 1000 0 100000])
    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);
    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,

    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
  2. jcsd
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted