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