#
# A Predator-prey problem is analyzed using Phase Planes.
# 
--------------------------------------------------------------------------------
# 
# Art Belmonte
# Thu, 11/Apr/96
# Math 308-509 (Worksheet B)
# 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/3a-1: critical points
> F:=(x, y)->x*(1 - x/2 - y/2); G:=(x, y)->y*(x/2 - 1/4);\
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 (1 - 1/2 x - 1/2 y)

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

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

                            dy_dt := y (1/2 x - 1/4)

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

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

               {y = 0, x = 0}, {y = 0, x = 2}, {x = 1/2, y = 3/2}
--------------------------------------------------------------------------------
# T-480/3b-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 := x

                                  LG := - 1/4 y

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

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

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

                                LF := - x + 2 - y

                                   LG := 3/4 y

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

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

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

                           LF := - 1/4 x + 1/2 - 1/4 y

                                LG := 3/4 x - 3/8

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

                                 1/2                        1/2
        evev := [- 1/8 + 1/8 I 11   , 1, {[ - 1/6 + 1/6 I 11   , 1 ]}],

                             1/2                        1/2
            [- 1/8 - 1/8 I 11   , 1, {[ - 1/6 - 1/6 I 11   , 1 ]}]

[ - .1250000000 + .4145780988 I, 1., {[  - .1666666667 + .5527707984 I, 1. ]}],

  [ - .1250000000 - .4145780988 I, 1., {[  - .1666666667 - .5527707984 I, 1. ]}]
--------------------------------------------------------------------------------
# T-480/3cd: Trajectories; also see (e) below.
> p1:=DEplot2(derivs, [x, y], -5..5,\
{[0, 1/2, 1/4], [0, 1/4, 1/2], [0, 1/2, 1/2],\
[0, 7/4, 1/4], [0, 2, 1/4], [0, 9/4, 1/4],\
[0, 3/4, 3/2], [0, 1/2, 7/4], [0, 1/4, 3/2], [0, 1/2, 5/4],\
[0, 1, 1], [0, 2, 1], [0, 3, 1], [0, 4, 1],\
[0, 1, 2], [0, 2, 2], [0, 3, 2], [0, 4, 2],\
[0, 1, 3], [0, 2, 3], [0, 3, 3], [0, 4, 3],\
[0, 1, 4], [0, 2, 4], [0, 3, 4], [0, 4, 4]},\
x=0..4, y=0..4,\
scaling=constrained, stepsize=0.05):\
display(p1);
--------------------------------------------------------------------------------
# T-480/3e: Line print for MATLAB input. It appears that any trajectory which starts 
# in the first quadrant goes to (2, 2) as t approaches infinity.
> derivs;\

--------------------------------------------------------------------------------
# 
>