#
# A Predator-Prey problem is explored using Phase Plane analysis.
# 
--------------------------------------------------------------------------------
# 
# Art Belmonte
# Thu, 11/Apr/96
# Math 308-509 (Worksheet A)
# Section 9.5: Predator-Prey Equations
# 
> with(linalg): with(student): with(DEtools): with(plots):\
readlib(mtaylor):
Warning: new definition for   norm
Warning: new definition for   trace

--------------------------------------------------------------------------------
# T-480/1a-1: critical points
> F:=(x, y)->x*(3/2 - y/2); G:=(x, y)->y*(x - 1/2);\
dx_dt:=F(x, y); dy_dt:=G(x,y);\
derivs:=[dx_dt, dy_dt];\
eqs:=equate(derivs, 0); solve(eqs, {x, y});

                          F := (x,y) -> x (3/2 - 1/2 y)

                            G := (x,y) -> y (x - 1/2)

                            dx_dt := x (3/2 - 1/2 y)

                              dy_dt := y (x - 1/2)

                    derivs := [x (3/2 - 1/2 y), y (x - 1/2)]

                  eqs := {x (3/2 - 1/2 y) = 0, y (x - 1/2) = 0}

                        {x = 0, y = 0}, {y = 3, x = 1/2}
--------------------------------------------------------------------------------
# T-480/1b-1: Via Table 9.3.1, critical point (0, 0) is an unstable saddle point.
> LF:=mtaylor(F(x,y), [x=0, y=0], 2);\
LG:=mtaylor(G(x,y), [x=0, y=0], 2);\
A:=genmatrix([LF, LG], [x, y]);\
evev:=eigenvects(A, radical); evalf(");

                                   LF := 3/2 x

                                  LG := - 1/2 y

                                    [ 3/2    0  ]
                               A := [           ]
                                    [  0   -1/2 ]

               evev := [3/2, 1, {[ 1, 0 ]}], [-1/2, 1, {[ 0, 1 ]}]

         [1.500000000, 1., {[ 1., 0 ]}], [-.5000000000, 1., {[ 0, 1. ]}]
--------------------------------------------------------------------------------
# T-480/1b-2: Via Table 9.3.1, critical point (1/2, 3) is either a center or an 
# asymptotically stable spiral point. We'll examine the trajectories in the next block to 
# see which of these it is.
> LF:=mtaylor(F(x,y), [x=1/2, y=3], 2);\
LG:=mtaylor(G(x,y), [x=1/2, y=3], 2);\
A:=genmatrix([LF, LG], [x, y]);\
evev:=eigenvects(A, radical); evalf(");

                                LF := 3/4 - 1/4 y

                                 LG := 3 x - 3/2

                                     [ 0  -1/4 ]
                                A := [         ]
                                     [ 3    0  ]

evev :=

          1/2               1/2                   1/2                 1/2
  [1/2 I 3   , 1, {[ 1/6 I 3   , 1 ]}], [- 1/2 I 3   , 1, {[ - 1/6 I 3   , 1 ]}]

             [.8660254040 I, 1., {[ .2886751347 I, 1. ]}],

                 [ - .8660254040 I, 1., {[  - .2886751347 I, 1. ]}]
--------------------------------------------------------------------------------
# T-480/1cd: Trajectories; also see (e) below. It appears that (1/2, 3) is a stable center.
> p1:=DEplot2(derivs, [x, y], -5..5,\
{[0, 1/2, 1/4], [0, 1/4, 1/2], [0, 1/2, 1/2],\
[0, 3/4, 3], [0, 1/2, 13/4], [0, 1/4, 3], [0, 1/2, 11/4],\
[0, 2, 1], [0, 2, 2], [0, 2, 3]},\
x=0..4, y=0..4,\
scaling=constrained, stepsize=0.05):\
display(p1);
--------------------------------------------------------------------------------
# T-480/1e: Line print for MATLAB input. The trajectories reveal that the two 
# populations vary cyclically. Also look at the x verses t and y versus t plots to see how 
# the oscillations of the predator population (y) lags behind that of the prey (x).
> derivs;
--------------------------------------------------------------------------------
# 
>