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

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