Maple Tutorial: Discrete Dynamical Systems -- (Adapted from VLA Module)

> with(LinearAlgebra):

Predator-Prey Model : Foxes and Rabbits

In this section we consider a dynamical systems model that simulates the interactions between two populations, a population of foxes (the predators) and a population of rabbits (the prey).  We denote the fox and rabbit populations at time k by x[k] = MATRIX([[F[k]], [R[k]]]) ,  where k is in years and where F[k] is the number of foxes after k years and R[k] is the number of rabbits (in hundreds) after k years. Also we suppose that:  

F[k+1] = Float(6, -1)*F[k]+Float(4, -1)*R[k]

R[k+1] = -Float(125, -3)*F[k]+Float(12, -1)*R[k]

The term Float(6, -1)*F[k] in the first equation says that in the absence of rabbits for food only 60% of the foxes would survive each year. On the other hand, the presence of rabbits increases the number of foxes, which is represented by the term Float(4, -1)*R[k] . Looking at the second equation, we see the term Float(12, -1)*R[k] , which represents a 20% growth of rabbits per year if there were no foxes. Finally, the effect of the foxes on the rabbits is represented by the term -Float(125, -3)*F[k] ; thus, the greater the number of foxes, the more the rabbit population will decrease in a given year.

Here is the transition matrix M for the fox and rabbit populations:

> M := Matrix([[6/10,4/10],[-125/1000,12/10]]);

M := Matrix([[3/5, 2/5], [(-1)/8, 6/5]])

Let's find the  eigenvalues and eigenvectors of M and use them to find an eigenvector decomposition of the points x[k] :

> Eigenvectors(M);

Vector[column]([[11/10], [7/10]]), Matrix([[4/5, 4], [1, 1]])


x[k] = a[1]*(11/10)^k*MATRIX([[Float(8, -1)], [1]])+a[2]*(7/10)^k*MATRIX([[4], [1]])

Since one eigenvalue is larger than 1 and the other is smaller than 1, the origin is a saddle point for this dynamical system. More specifically, since x[k] is approximately equal to a[1]*(11/10)^k*MATRIX([[Float(8, -1)], [1]]) for large values of k,  both the fox and rabbit populations grow at a yearly rate of 10% (except in the special case when a[1] = 0 ), and the populations are in the ratio of 8 foxes to every 1000 rabbits. Note, moreover, that these conclusions are independent of the initial population x[0] . Let's look at the trajectory and the component plots in which we start with 50 foxes and 1500 rabbits:

> x0 := Vector([50,15]);

x0 := Vector[column]([[50], [15]])


In the above plot, observe that the fox population decreases dramatically at first, as does the rabbit population, and the rabbits recover first, before both populations eventually get larger and larger. Note the contrasting situation below where we start with 5 foxes and 1000 rabbits. Give a plausible explanation for why these different initial conditions lead to such different trajectories.

> x0 := Vector([5,10]);

x0 := Vector[column]([[5], [10]])