Ode 45
Ode 45
Ode 45
9.2
We consider the predator to be a fox and a prey to be either rabbits or turkeys. One could also consider different species of fish such as sharks and bass.
One Predator, and One Prey with Constant Birth and Death Rate. Suppose a fox population only eats rabbits and has population of x(t) at time t. With no rabbits the fox population will die out according to the differential equation x' = -dx where d is the birth rate minus the death rate. The rabbit population y(t) will, without any foxes, continue to grow for some time according to the differential equation y' = by where b is the birth rate minus the death rate.
We want to model the case where both populations are positive. In the above differential equation for the foxes we need an additional term on the right side that will be positive when there are rabbits. More precisely, this additional term should be large if either x(t) or y(t) are large. There are many ways to do this, but a good option is modify the birth rate minus the death rate for the fox populations to be -d + ey. The e is a positive constant because the rabbits are a source of food for the foxes. Thus, the differential equation for the foxes is x' = (-d + ey)x and x(0) = x0. Fox Equation
A similar term must be inserted into the birth rate minus the death rate for the rabbits. But, in this case the rabbit population should decrease, and so, we will insert -cx. y' = (b - cx)y and y(0) = y0. Rabbit Equation
These two differential equations are coupled, and they are known as the Lotka-Volterra system of differential equations. Solutions to such systems are usually difficult to find, and so, numerical methods are used. We will discuss this below. Fortunately, one can gain some qualitative knowledge about the solutions. There should be oscillations in the population. If the rabbit population is very large, then there will be plenty of food for the fox population and it should increase. As the fox population increases, more of the rabbits are consumed. So, the food supply for the fox population decreases, and the fox population will have to decrease. One expects there is a perfect balance of the populations where neither population changes. This means both x' and y' must be zero. This is called the steady state or equilibrium solution. There is one positive steady state solution x = b/c from the rabbit equation, and y = d/e from the fox equation.
In fact, the graph of the points (x(t), y(t)) must be a closed curve. In order to see this, consider the composed function H(x(t), y(t)), where H(x, y) = b ln(x) - c x + d ln(y) ey and x(t) and y(t) satisfy the Lotka-Volterra equation. The time derivative of this composed function is computed by the chain rule as follows H' = Hx x' + Hy y' = (b/x -c) (-d + e y) x + (d/y - e) (b - c x) y = 0. Thus, H(x(t), y(t)) must be a constant. With some additional effort one can show that H(x, y) = constant defines a simple closed curve. Any numerical method should generate such a simple closed curve.
9.3
Method of Solution.
In this lesson we will use the Matlab command ode45 to solve our systems of
differential equation. This command is a robust implementation for systems of differential equations, which uses a variable step size method and the fourth and fifth order Runge-Kutta method.
9.4
Matlab Implementation.
In the following calculations we have attempted to solve the second predator-prey
model
The steady state solutions are given by x' = -.5x + .01xy = 0 and y' = .5y - .01xy = 0.
The nonzero steady state solutions are x = 50 and y = 50. Our numerical solution should orbit about the steady state solution. The first graph below does not quite do this. This is because some of the error options for ode45 need to be made smaller. The default value for relative error is .001, and this is what was used in the first graph below. The Matlab command odeset('RelTol', 1e-5) was used to "tune" the ode45 method. Adjust the rf.m file by replacing line 6, the line with ode45, by the two following lines: options = odeset('RelTol', 1e-5); [t y] = ode45('yprf',[to tf],yo,options);
Here more time steps will be required; use the whos command to find the number of time steps n = 133 with RelTol = 1e-3, and n = 289 with RelTol = 1e-5.
m-file yprf.m: function yprf =yprf(t,y) yprf(1) =-.5*y(1) + .01*y(1)*y(2); yprf(2) = .5*y(2) -.01*y(1)*y(2); yprf = [yprf(1) yprf(2)]';
m-file rf.m: %your name, your student number, lesson number clear; to = 0; tf =50; yo = [80 100]; [t y] = ode45('yprf',[to tf],yo); plot(y(:,1),y(:,2)) title('your name, your student number, lesson number') ylabel('rabbits') xlabel('fox') %plot(t,y(:,1),t,y(:,2)) %xlabel('time') %ylabel('rabbits and fox')
100
80
rabbits
60
40
20
100
80
rabbits
60
40
20
In the yprf.m file the x(t), fox, is associated with the symbol y(1), and y(t), rabbits, is associated with the symbol y(2). The right sides of the differential equations are given by yprf(1) and yprf(2). In rf.m the yo = [80 100] is a 1x2 array of initial populations where x(0) = 80 and y(0) = 100. The output [t y] from the rf.m file has three column vectors where the first column is for all the time values and the second and third columns are for the corresponding values of x(t) and y(t). These values can be viewed by typing [t y] at the Matlab prompt. Also the following graph of x(t) versus t and y(t) versus t was generated by typing plot(t, y(:,1),t,y(:,2)) at the Matlab prompt.
120
100
80
60
40
20
0 0 5 10 15 20 25 time 30 35 40 45 50
9.5
Numerical Experiments.
In the following calculations we have attempted to solve the fox and rabbit
problem with harvesting of the fox population. The fox equation must be modified as follows to include a harvesting rate h(t) fox per unit time. x' = -.5x + .01xy - h(t) y' = .5y - .01xy . In the calculation indicated in the following graph h(t) = 1, and the graph spirals outward.
your name, your student number, lesson number 250
200
150
rabbits
100 50 0 0
50
100 fox
150
200
250
9.6
Additional Calculations.
Consider the following fox and rabbit model with no harvesting x' = -.5x + .01xy, y' = .4y - .01xy and
with variable initial conditions [x(0) y(0) ] = [(100 +10*(1+.S)) 90 ], [(100 +10*(1+.S)) 1 ] and [(100 +10*(1+.S)) 40]. (a). (b). (c). Find a steady state solution where both populations are positive. Modify the yprf.m file. Insert you name and student number into the rf.m file. Execute rf.m for each set of initial conditions. (d). Print the yprf.m and rf.m files and the graphs.