#
# A competing species problem is studied using Phase Planes.
#
# Art Belmonte
# Mon, 27/May/96
# Math 308-509 [Maple V Release 4]
# 
# Section 9.4: Competing Species
# 
# T-471/6a-1: critical points
> restart;
> with(linalg): with(student): with(DEtools): with(plots):
Warning, new definition for norm
Warning, new definition for trace
> readlib(mtaylor):
> F:=(x, y)->x*(1 - x + y/2); G:=(x, y)->y*(5/2 - 3/2*y + 1/4*x);

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


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

> dx_dt:=F(x, y); dy_dt:=G(x,y);

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


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

> derivs:=[dx_dt, dy_dt];

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

> eqs:=equate(derivs, 0); solve(eqs, {x, y});

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


   {y = 0, x = 0}, {x = 0, y = 5/3}, {y = 0, x = 1}, {y = 2, x = 2}

> deq:=equate([diff(x(t), t), diff(y(t), t)],
>      derivs(t));

        d
deq := {-- x(t) = x(t) (1 - x(t) + 1/2 y(t)),
        dt

    d
    -- y(t) = y(t) (5/2 - 3/2 y(t) + 1/4 x(t))}
    dt

# T-471/6a-2: slanted nullclines
> F(x,y)/x=0; solve(", {y});

                          1 - x + 1/2 y = 0


                            {y = -2 + 2 x}

> n1:=rhs(op("));

                            n1 := -2 + 2 x

> G(x,y)/y=0; solve(", {y});

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


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

> n2:=rhs(op("));

                          n2 := 5/3 + 1/6 x

> p0:=plot({n1, n2}, x=0..4, y=0..4, thickness=3):
> display(p0);

# 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);

                               LF := x

> LG:=mtaylor(G(x,y), [x=0, y=0], 2);

                             LG := 5/2 y

> A:=genmatrix([LF, LG], [x, y]);

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

> evev:=eigenvects(A, radical); evalf(");

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


          [1., 1., {[1., 0]}], [2.500000000, 1., {[0, 1.]}]

# 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);

                         LF := 1 - x + 1/2 y

> LG:=mtaylor(G(x,y), [x=1, y=0], 2);

                             LG := 11/4 y

> A:=genmatrix([LF, LG], [x, y]);

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

> evev:=eigenvects(A, radical); evalf(");

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


     [-1., 1., {[1., 0]}], [2.750000000, 1., {[1., 7.500000000]}]

# 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);

                             LF := 11/6 x

> LG:=mtaylor(G(x,y), [x=0, y=5/3], 2);

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

> A:=genmatrix([LF, LG], [x, y]);

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

> evev:=eigenvects(A, radical); evalf(");

         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);

                          LF := -2 x + 2 + y

> LG:=mtaylor(G(x,y), [x=2, y=2], 2);

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

> A:=genmatrix([LF, LG], [x, y]);

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

> evev:=eigenvects(A, radical); evalf(");

                      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. (This takes LONG to
# execute!)
> inits:={[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]};

inits := {[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, 15/8, 15/8], [0, 15/8, 17/8], [0, 17/8, 15/8],

    [0, 17/8, 17/8], [0, 3, 3], [0, 4, 4]}

> p1:=DEplot(deq, [x(t), y(t)], t=-5..5, x=0..4, y=0..4,
>      inits, 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;
# 
>