# COMPETING SPECIES # A NONLINEAR AUTONOMOUS SYSTEM # # W. E. Boyce # Rensselaer Polytechnic Institute # 04/15/94 # # PURPOSE. To show how a two-dimensional nonlinear autonomous system can be # investigated systematically. # # EXAMPLE. Consider the system of differential equations # # dx/dt = x(3 - x - y), dy/dt = y(2 - y - x/2). # # We can interpret x(t) and y(t) in this system as the populations of two # species that compete with each other in a closed ecological system. We are # interested primarily in the long-term behavior of x(t) and y(t). In # particular, at what levels, if any, can they coexist in a stable # relationship? # # (a) Find all critical points (equilibrium solutions). # # (b) Plot a direction field on a rectangle in the first quadrant that # includes all of the critical points. # # (c) Plot several trajectories of the system. # # (d) Plot the time response of x versus t. # # (e) Find the linear system that approximates the given system near the # critical point in the positive first quadrant. # # (f) Find the eigenvalues and eigenvectors of the coefficient matrix of this # linear system, and determine its stability characteristics. # # # SOLUTION. First we bring up the package DEtools; the colon after the # command suppresses the screen output, i.e., the list of commands in the # package. # > with(DEtools): -------------------------------------------------------------------------------- # # There are several ways to enter the system of equations. We will enter the # right side of each equation, calling these expressions f and g, # respectively. # > f := x*(3-x-y); g := y*(2-y-x/2); f := x (3 - x - y) g := y (2 - y - 1/2 x) -------------------------------------------------------------------------------- # # (a) To find the critical points, we must solve the equations f=0 and g=0 # simultaneously. # > crit := solve({f,g},{x,y}); crit := {x = 0, y = 0}, {x = 0, y = 2}, {y = 0, x = 3}, {y = 1, x = 2} -------------------------------------------------------------------------------- # # We observe that there is one critical point for which both x and y are # positive. This critical point corresponds to coexistence. For the other # critical points one or both species become extinct. # # (b) Next we plot a direction field on a rectangle that includes all of the # critical points. We choose the square x=0..4, y=0..4. The following # command produces the direction field. Note that the third argument in the # command gives a t interval, and this argument is required, even though the # direction field is independent of t. # > DEplot([f,g],[x,y],t=0..1,x=0..4,y=0..4,arrows=THIN,title=`Figure 1`); -------------------------------------------------------------------------------- # # ** Maple V Graphics ** # # From the direction field it appears that the critical point (2,1) is an # asymptotically stable node, and that the other critical points are unstable. # # (c) In order to plot trajectories we must assign some initial points. The # following command is one way to do this. # > init1:= seq([0,0.5,.4*i],i=0..6); init2:=seq([0,3,0.4*i],i=0..6); init3:=seq([0,0.4*i,0.5],i=0..6); init1 := [0, .5, 0], [0, .5, .4], [0, .5, .8], [0, .5, 1.2], [0, .5, 1.6], [0, .5, 2.0], [0, .5, 2.4] init2 := [0, 3, 0], [0, 3, .4], [0, 3, .8], [0, 3, 1.2], [0, 3, 1.6], [0, 3, 2.0], [0, 3, 2.4] init3 := [0, 0, .5], [0, .4, .5], [0, .8, .5], [0, 1.2, .5], [0, 1.6, .5], [0, 2.0, .5], [0, 2.4, .5] -------------------------------------------------------------------------------- # # The following command superimposes the trajectories passing through the # given sets of points on the direction field shown in Figure 1. # > DEplot([f,g],[x,y],t=0..5,{init1,init2,init3},arrows=THIN,title=`Figure 2`); -------------------------------------------------------------------------------- # # ** Maple V Graphics ** # # (d) To plot x versus t we can use a command such as the following. # Note the presence of the scene option. The initial conditions were chosen # arbitrarily. # > DEplot([f,g],[x,y],t=0..5,{[0,1/2,1/2],[0,3,2]},scene=[t,x],title=`Figure 3`); # # ** Maple V Graphics ** # Both solutions appear to be approaching the value two, consistent with the # phase portrait. # # In a similar way one can also plot y versus t. # # (e) To find the linear system valid near the critical point (2,1) it is # helpful to bring up the linalg package. Since this package contains a great # many commands, we save space by ending the next command with a colon. If # you wish to see a list of the commands in the linalg package, replace the # colon by a semicolon. # > with(linalg): -------------------------------------------------------------------------------- # # To find the required linear system we need to find the matrix of partial # derivatives of f and g with respect to x and y, and then evaluate the # elements of this matrix at the critical point (2,1). The following three # commands perform these steps. # > A := vector([f,g]); A := [ x (3 - x - y), y (2 - y - 1/2 x) ] -------------------------------------------------------------------------------- # # Note that there is only a space separating the components of the vector A. # Thus you need to be careful in reading output such as this. # > B := jacobian(A,[x,y]); [ 3 - 2 x - y - x ] B := [ ] [ - 1/2 y 2 - 2 y - 1/2 x ] -------------------------------------------------------------------------------- > C1 := subs(x=2,y=1,op(B)); [ -2 -2 ] C1 := [ ] [ -1/2 -1 ] -------------------------------------------------------------------------------- # # Note the presence of the `op' in the last command. # # (f) The analysis of this linear system involves finding the eigenvalues and # eigenvectors of the matrix C1. Both can be found at the same time by use of # the command `eigenvects'. The option `radical' directs Maple to express (if # possible) the eigenvalues and eigenvectors in terms of radicals. # > q1 := eigenvects(C1,'radical'); 1/2 1/2 q1 := [- 3/2 + 1/2 5 , 1, {[ 1, - 1/4 - 1/4 5 ]}], 1/2 1/2 [- 3/2 - 1/2 5 , 1, {[ 1, - 1/4 + 1/4 5 ]}] -------------------------------------------------------------------------------- # # To interpret Maple's last response observe that it consists of two lists. # In each list the first item is an eigenvalue, the second item is the # multiplicity of the eigenvalue, and the third item is the corresponding # eigenvector. In reading the eigenvectors note that the first component is # one, with the following space separating the components. # # To obtain a numerical answer we can use `evalf', as shown in the following # command. # > evalf(q1[1]); evalf(q1[2]); [-.381966011, 1., {[ 1., -.8090169945 ]}] [-2.618033989, 1., {[ 1., .3090169945 ]}] -------------------------------------------------------------------------------- # # Since the eigenvalues are both negative, the critical point is an # asymptotically stable node of the linear system, and the same is true for # the nonlinear system. Further, trajectories approach the critical point # tangent to the eigenvector associated with the larger eigenvalue (-0.38... # in this case). # # To see tis graphically we can plot a direction field for the linear system. # > DEplot(C1,[x,y],t=0..1,{[0,0,0]},x=-1..1,y=-1..1,arrows=THIN,title=`Figure 4`); -------------------------------------------------------------------------------- # # ** Maple V Graphics ** # # Compare this direction field with the portion of Figure 1 near the point # (2,1). # # Of course, one can proceed in the same way to obtain a linear approximation # near each of the critical points of the original system. # #