Assignment given 28 April 2009

This assignment was given on 28 April 2009. The work is due 05 May 2009.

So far, our SIR model of disease has been discrete in time and space. Simple discrete models are useful for trying to understand a problem; if a discrete model turns out to be useful, it is often important to see if it can be turned into a continuous model; in particular, we are tempted to replace finite differences by derivatives.

We are going to use MATLAB's ODE solver ode23 for this assignment. Let's go over an example of using this solver.

Our example ODE is to be solved from an initial time T0 = 10 to a final time TMAX = 100. There are N = 2 equations, for variables U(T) and V(T). At our initial time T0, we have U(T0) = 17 and V(T0) = 33. The ODE is

        dU/dt = - V
        dV/dt =   U + U^2
      

To use ode23, we must write a MATLAB function that evaluates the differential equation. This function will receive as input the current time T, and the values of the variables U and V. The variables will be stored in a vector. The function must compute the derivatives and return them as a vector. MATLAB wants this output vector to be a column vector, so we have to write it with two indices.

      function dydt = example ( t, y )

        dydt(1,1) = - y(2)
        dydt(2,1) =   y(1) + y(1) * y(1)

        return
      end
      

Now we need to call ode23, giving it the name of our function, the time range, and the initial values. The solution is returned as a list of times TOUT(:), and a table of variable values YOUT(:,1:2). The number of values MATLAB returns depends on how hard it has to work to compute the answer and follow the solution curve.

        t0 = 10;
        tmax = 100;
        u0 = 17;
        v0 = 33;

        [ tout, yout ] = ode23 ( @example, [ t0, tmax ], [ u0, v0 ] );

        plot ( tout, yout(:,1), 'g', tout, yout(:,2), 'r' );
      

Question 1: our discrete SIR model produced, at each time step, quantities S(T), I(T), and R(T), which were either counts of the number of patients in each class, or normalized values that added to 1. If we want to study a large population, we don't care so much about the individual patients, but just these overall numbers.

we suppose that S, I and R can be thought of as continuous functions of time. If we know the initial values of these quantities, and if we can provide rules for how they change, we have an initial value problem for three ordinary differential equations.

        d S(t) /dt = - tau I(t) S(t)
        d I(t) /dt =   tau I(t) S(t) - I(t) / k
        d R(t) /dt =                   I(t) / k
      

Solve the system of differential equations using ode23. We start at time t = 0, and we want to run til tmax = 50. Our initial values are S(0) = 0.995, I(0) = 0.005, R(0) = 0. Use tau = 0.8 and k = 4.

Turn in: Plot the three functions S(t), I(t) and R(t), for t from 0 to 50.

Question 2: Let us recall the predator prey equations which we worked on earlier, using finite differences:

        B(I+1) = B(I) + DT * (    2 * B(I) - 0.001 * B(I) * F(I) )
        F(I+1) = F(I) + DT * ( - 10 * F(I) + 0.002 * B(I) * F(I) )
      

If we subtract B(I) from both sides of equation 1, and divide by DT and take the limit as DT goes to zero, we get a differential equation. Doing the same thing for the second equation, we have:

        dB/dt =    2 * B(I) - 0.001 * B(I) * F(I)
        dF/dt = - 10 * F(I) + 0.002 * B(I) * F(I)
      

Write a MATLAB program that calls ode23 to solve this system, from t = 0 to t = 5, using initial values b(0) = 5000, f(0) = 100.

Turn in: Plot the functions B(t) and F(t) for t from 0 to 5.