#
# A Competing species system is analyzed using Phase Planes
#
# 
--------------------------------------------------------------------------------
# 
# Art Belmonte
# Tue, 09/Apr/96
# Math 308-509
# Section 9.4: Competing Species
# 
> with(linalg): with(student): with(DEtools): with(plots):\
readlib(mtaylor):
Warning: new definition for   D
Warning: new definition for   Int
Warning: new definition for   Sum

--------------------------------------------------------------------------------
# T-471/6a-1: critical points
> F:=(x, y)->x*(1 - x + y/2); G:=(x, y)->y*(5/2 - 3/2*y + 1/4*x);\
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 - x + 1/2 y)

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

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

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

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

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

        {x = 0, y = 0}, {x = 0, y = 5/3}, {y = 0, x = 1}, {x = 2, y = 2}
--------------------------------------------------------------------------------
# T-471/6a-2: slanted nullclines
> F(x,y)/x=0; solve(", {y});\
n1:=rhs(op("));\
G(x,y)/y=0; solve(", {y});\
n2:=rhs(op("));\
p0:=plot({n1, n2}, x=0..4, y=0..4, thickness=3):\
display(p0);

                                1 - x + 1/2 y = 0

                                 {y = - 2 + 2 x}

                                 n1 := - 2 + 2 x

                             5/2 - 3/2 y + 1/4 x = 0

                                {y = 5/3 + 1/6 x}

                                n2 := 5/3 + 1/6 x
--------------------------------------------------------------------------------
# T-471/6b-1: Via Table 9.3.1, critical point (0, 0) is an unstable improper node.
> 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 := 5/2 y

                                      [ 1   0  ]
                                 A := [        ]
                                      [ 0  5/2 ]

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

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

                               LF := 1 - x + 1/2 y

                                  LG := 11/4 y

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

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

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

                                  LF := 11/6 x

                          LG := - 5/2 y + 25/6 + 5/12 x

                                    [ 11/6    0  ]
                               A := [            ]
                                    [ 5/12  -5/2 ]

             evev := [-5/2, 1, {[ 0, 1 ]}], [11/6, 1, {[ 52/5, 1 ]}]

    [-2.500000000, 1., {[ 0, 1. ]}], [1.833333333, 1., {[ 10.40000000, 1. ]}]
--------------------------------------------------------------------------------
# T-471/6b-4: Via Table 9.3.1, critical point (2, 2) is an asymptotically stable 
# improper node.
> LF:=mtaylor(F(x,y), [x=2, y=2], 2);\
LG:=mtaylor(G(x,y), [x=2, y=2], 2);\
A:=genmatrix([LF, LG], [x, y]);\
evev:=eigenvects(A, radical); evalf(");

                               LF := - 2 x + 2 + y

                             LG := - 3 y + 5 + 1/2 x

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

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

                                 1/2             1/2
                   [- 5/2 - 1/2 3   , 1, {[ 1 - 3   , 1 ]}]

                  [-1.633974596, 1., {[ 2.732050808, 1. ]}],

                      [-3.366025404, 1., {[ -.732050808, 1. ]}]
--------------------------------------------------------------------------------
# T-471/6cd-1: 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, 3/4, 1/4], [0, 1, 1/4], [0, 5/4, 1/4],\
[0, 1/3, 4/3], [0, 1/3, 5/3], [0, 1/3, 2],\
[0, 2-1/8, 2-1/8], [0, 2-1/8, 2+1/8],\
     [0, 2+1/8, 2-1/8], [0, 2+1/8, 2+1/8],\
[0, 3, 3], [0, 4, 4]},\
x=0..4, y=0..4,\
scaling=constrained, stepsize=0.05):\
display(p1);
--------------------------------------------------------------------------------
# T-471/6cd-2: Nullclines superimposed upon trajectories.
> display([p0, p1]);
--------------------------------------------------------------------------------
# T-471/6e: 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;
--------------------------------------------------------------------------------
# 
>