#
# A predator-prey problem is studies using Phase Planes.
#
# Art Belmonte
# Mon, 27/May/96
# Math 308-509 [Maple V Release 4] (Worksheet B)
# Section 9.5: Predator-Prey Equations
# 
# T-480/3a-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/2 - y/2); G:=(x, y)->y*(x/2 - 1/4);

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


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

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

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


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

> derivs:=[dx_dt, dy_dt];

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

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

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


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

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

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

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

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

                               LF := x

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

                            LG := - 1/4 y

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

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

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

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


          [1., 1., {[1., 0]}], [-.2500000000, 1., {[0, 1.]}]

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

                           LF := 2 - x - y

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

                             LG := 3/4 y

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

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

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

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


    [-1., 1., {[1., 0]}], [.7500000000, 1., {[1., -1.750000000]}]

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

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

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

                          LG := 3/4 x - 3/8

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

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

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

                         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.
> inits:={[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]};

inits := {[0, 1/4, 1/2], [0, 1/2, 1/2], [0, 1/2, 1/4], [0, 1, 4],

    [0, 2, 4], [0, 3, 4], [0, 3, 3], [0, 4, 3], [0, 1, 3], [0, 2, 3],

    [0, 4, 1], [0, 1, 2], [0, 2, 2], [0, 3, 2], [0, 4, 2], [0, 1, 1],

    [0, 2, 1], [0, 3, 1], [0, 3/4, 3/2], [0, 1/2, 5/4], [0, 1/2, 7/4],

    [0, 1/4, 3/2], [0, 2, 1/4], [0, 9/4, 1/4], [0, 7/4, 1/4],

    [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):
Warning, computation interrupted
> 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;

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

# 
>