#           COMPETING SPECIES
#     A NONLINEAR AUTONOMOUS SYSTEM
# 
# W. E. Boyce
# Rensselaer Polytechnic Institute
# This file was acquired from the RPI ftp site: ftp.rpi.edu
# 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.
# 
# ---------------
# 
# This file has been slightly modified to work in Maple V R4
# 
# 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 hand sides of the two differential equations as f and
# g, respectively and the entire DEs as F and G.
# 
> f := x(t)*(3-x(t)-y(t));
> g := y(t)*(2-y(t)-x(t)/2);
> F := diff(x(t),t) = f;
> G := diff(y(t),t) = g;

                     f := x(t) (3 - x(t) - y(t))


                   g := y(t) (2 - y(t) - 1/2 x(t))


                     d
                F := -- x(t) = x(t) (3 - x(t) - y(t))
                     dt


                   d
              G := -- y(t) = y(t) (2 - y(t) - 1/2 x(t))
                   dt

# 
# (a)  To find the critical points, we must solve the equations  f=0 
# and  g=0 simultaneously.
# 
> crit := solve({f=0,g=0},{x(t),y(t)});

crit := {x(t) = 0, y(t) = 0}, {x(t) = 0, y(t) = 2},

    {x(t) = 3, y(t) = 0}, {x(t) = 2, y(t) = 1}

# 
# 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(t),y(t)},
> t=0..1,x=0..4,y=0..4,arrows=THIN,title=`Figure 1`);

# 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( [ x(0)=0.5, y(0)=0.4*i], i=0..6);
> init2 := seq( [ x(0)=3, y(0)=0.4*i], i=0..6);
> init3 := seq( [ x(0)=0.4*i, y(0)=0.5], i=0..6);
# 
# The following command superimposes the trajectories passing through
# the given sets of points on the direction field shown in Figure 1.

init1 := [x(0) = .5, y(0) = 0], [x(0) = .5, y(0) = .4],

    [x(0) = .5, y(0) = .8], [x(0) = .5, y(0) = 1.2],

    [x(0) = .5, y(0) = 1.6], [x(0) = .5, y(0) = 2.0],

    [x(0) = .5, y(0) = 2.4]


init2 := [x(0) = 3, y(0) = 0], [x(0) = 3, y(0) = .4],

    [x(0) = 3, y(0) = .8], [x(0) = 3, y(0) = 1.2],

    [x(0) = 3, y(0) = 1.6], [x(0) = 3, y(0) = 2.0],

    [x(0) = 3, y(0) = 2.4]


init3 := [x(0) = 0, y(0) = .5], [x(0) = .4, y(0) = .5],

    [x(0) = .8, y(0) = .5], [x(0) = 1.2, y(0) = .5],

    [x(0) = 1.6, y(0) = .5], [x(0) = 2.0, y(0) = .5],

    [x(0) = 2.4, y(0) = .5]

> DEplot( {F,G},{x(t),y(t)}, t=0..5, {init1,init2,init3} );

# (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(t),y(t)},t=0..5,{[x(0)=1/2,y(0)=1/2],[x(0)=3,y(0)=2]},
> scene=[t,x(t)]);
Error, (in DEtools/DEplot/CheckDE) invalid subscript selector
# 
# 
# 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):
Warning, new definition for norm
Warning, new definition for trace
# 
# 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(subs(x(t)=x,y(t)=y,[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 5    - 1/4]}]

# 
# 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.
# 
> C2:=multiply(C1,[x(t),y(t)] );

             C2 := [-2 x(t) - 2 y(t), - 1/2 x(t) - y(t)]

> DEplot({diff(x(t),t)=C2[1],diff(y(t),t)=C2[2]},[x(t),y(t)],
> t=0..1,{[x(0)=0,y(0)=0]},x=-1..1,y=-1..1,arrows=THIN,title=`Figure
> 4`);

> 
# 
# 
# 
# 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.
# 
#