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

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


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

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

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


                         dy_dt := y (x - 1/2)

> derivs:=[dx_dt, dy_dt];

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

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

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


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

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

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

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

                             LF := 3/2 x

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

                            LG := - 1/2 y

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

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

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

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


     [-.5000000000, 1., {[0, 1.]}], [1.500000000, 1., {[1., 0]}]

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

                          LF := 3/4 - 1/4 y

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

                           LG := 3 x - 3/2

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

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

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

                1/2      [       1/2   ]
evev := [1/2 I 3   , 1, {[1/6 I 3   , 1]}],

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

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

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

    [0, 1/2, 1/4], [0, 1/4, 1/2]}

> 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/1e: MatLab input
> derivs;
>