1. Not finding help here? Sign up for a free 30min 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!

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
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?



Similar Discussions: Modelling the spread of disease (in MatLab, Runge-Kutta)
  1. Runge kutta (Replies: 0)

  2. Help with Matlab (Replies: 0)

  3. Matlab Problem (Replies: 0)

  4. MATLAB: Optimization (Replies: 0)

Loading...