MatLab Assignment 6 -- Turn in on April 24

-- MatLab Functions: ode23, ode45 to Solve Systems of Differential Equations
If you have any question or problem with MatLab, please send me an email at: chenm or call me at: 953-7896.


Consider the problem of modeling populations of two species, one of which is predator (say, wolf), whose population at time t is y1(t), feeding on the other, which is called a prey (say, rabbit), whose population is y2(t). We may assume that the prey always has an adequate food supply and its birth rate at any time is proportional to the number of prey alive at that time; that is, birth rate(prey)=k4y2(t). The death rate of the prey depends on both the number of prey and predators alive at that time. For simplicity, we may assume the death rate(prey)=k3y1(t)y2(t). The birth rate of the predator, on the other hand, depends on its food supply, y2(t), as well as on the number of predators available for reproduction purposes. For this reason, we may assume that the birth rate(predator)=k2y1(t)y2(t). The death rate of the predator will be taken as simply proportional to the number of predators alive at the time; that is, death rate(predator)=k1y1(t).
The rates of change in the predator and prey populations can be expressed as:
     y'1=-k1 y1+k2y1y2,
     y'2=k4y2-k3y1y2,
     y1(0)=y01, y2(0)=y02.
Recall that algebraically we can solve this system of two nonlinear differential equations in terms of y1 and y2, but not y1(t) and y2(t) in t explicitly. MatLab has a function called ode45.m (same as ode23) which solves this system numerically. Type help ode45 to learn more about this function. Here is an example which uses ode45.m to solve the following system
     y'1=-0.5y1+0.0006y1y2,
     y'2=3y2-0.002y1y2, 
     y1(0)=100,y2(0)=1000, for t in [0,20].
  1. First we create a MatLab function called funrw.m which evaluates the values of y1(t) and y2(t) for a given t. File funrw.m looks like the following:
  2.      function uv=funrw(t,y)
         k1=0.5; k2=.0006; k3=.002; k4=3;
         uv(1)=-k1*y(1)+k2*y(1)*y(2);
         uv(2)=k4*y(2)-k3*y(1)*y(2);
  3. Now we use MatLab function ode45.m to solve the system numerically with the initial condition: y1=100 and y21000, and then plot the populations of predator y1(t), and prey y2(t) with respect to time t.
  4.    clf
       clear
       [tv yv1]=ode45('funwr',[0 20],[100;1000]);
       plot(tv,yv1(:,1),':',tv,yv1(:,2),'--')
       axis([0 20 0 10000])
       text(4,7500,'y_1(0)=100')
       text(4,4000,'y_2(0)=1000')
       title('.. Wolves Population, -- Rabbits Population')
    
    Function ode45.m generates a time vector tv and a two-column vector yv whose two columns contain values of y1(t) and y2(t), respectively. Graphs of y1 and y2 are given at the end of the assignment

    Problem: Consider the system of differential equations with the following initial value conditions:

        (a) y1(0)=100, y2(0)=1000;
        (b) y1(0)=200, y2(0)=1000;
        (c) y1(0)=400, y2(0)=1000.
        (d) y1(0)=800, y2(0)=1000.
    
     
  5. Now let us consider the problem of modeling the populations of two species (say, fox and wolf) which compete for same food supply. Let the populations of the species alive at time t be y1(t) and y2(t), respectively. We may assume that the birth rate of each of the species is simply proportional to the number of the species alive at that time and the death rate of each species depends on the population of both species. Now let us assume that the population of a particular pair of species is described by the equations:
  6.  
        y'1=4y1-0.0003y12-0.0004y1y2
        y'2=2y2-0.0002y1y2-0.0001y2y22
        y1(0)=y01,y2(0)=y02 for t in [0,12].
    Problem: Use ode45.m to solve the system with the following initial conditions:
        (a) y1(0)=1000, y2(0)=1000;
        (b) y1(0)=2000, y2(0)=2000;
        (c) y1(0)=3000, y2(0)=3000;
        (d) y1(0)=4000, y2(0)=4000.
    
    For each case, plot the solutions y1(t) and y2(t) for t in [0,12], and answer the following questions from the graphs or from the equations.
Here is the graph: