#
# Phase planes are graphed for eigenvalues of the form
#   1 pos, 1 neg
#   complex
#   imaginary
#
# Art Belmonte
# Mon, 27/May/96
# Math 308-509 [Maple V Release 4]
# 
# Section 9.1: The Phase Plane: Linear Systems
# 
# T-437/1-1: But let's try Maple first (static versus interactive phase
# portraits). Here's our constant matrix A, its eigenvalues, and the
# resulting phase portrait. Next, let's use pplane under MATLAB to
# interactive pump these out in RAPID fashion! Either way, we're
# interested in the QUALITATIVE, GRAPHICAL behavior of solution
# trajectories. Here vector x = vector 0 is an unstable saddle point
# (see Table 9.1.1, page T-436).
> with(linalg): with(DEtools): with(student):
> A:=matrix(2, 2, [3, -2, 2, -2]); evals:=eigenvals(A);

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


                            evals := -1, 2

> deq:=equate(matrix([[diff(x1(t), t)], [diff(x2(t), t)]]),
> evalm(A &* matrix([[x1(t)], [x2(t)]])));

         d                             d
 deq := {-- x1(t) = 3 x1(t) - 2 x2(t), -- x2(t) = 2 x1(t) - 2 x2(t)}
         dt                            dt

# 
> inits:={[0, 1, -2], [0, -2, 3], [0, 1.6, 2.7],
>      [0, -1.3, -2]};

   inits := {[0, -2, 3], [0, 1.6, 2.7], [0, -1.3, -2], [0, 1, -2]}

> DEplot(deq, [x1(t), x2(t)], t=-5..5, x1=-4..4, x2=-4..4, inits,
> scaling=constrained);

# T-437/1-2: ONE trajectory constructed via analytical solution and
# plot. It agrees with the forgoing.
> A:=matrix(2, 2, [3, -2, 2, -2]);

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

> x0:=matrix(2, 1, [1, -2]);

                                    [ 1]
                              x0 := [  ]
                                    [-2]

> x:=evalm(exponential(A*t) &* x0);

                      [- 5/3 exp(-t) + 8/3 exp(2 t)]
                 x := [                            ]
                      [4/3 exp(2 t) - 10/3 exp(-t) ]

> xa:=unapply(convert(convert(evalm(x), vector),
>      list), t);

xa :=

    t -> [- 5/3 exp(-t) + 8/3 exp(2 t), 4/3 exp(2 t) - 10/3 exp(-t)]

> plot([op(xa(t)), t=-5..5], 'x'=-4..4, 'y'=-4..4, scaling=constrained);

# T-437/5-1: You know the drill. Let's increase production AND explore
# "complex" territory. Here x=0 is an asymptotically stable spiral point
# (see Table 9.1.1, page T-436).
> A:=matrix(2, 2, [1, -5, 1, -3]); evals:=eigenvals(A);

                                 [1    -5]
                            A := [       ]
                                 [1    -3]


                       evals := -1 + I, -1 - I

> deq:=equate(matrix([[diff(x1(t), t)], [diff(x2(t), t)]]),
> evalm(A &* matrix([[x1(t)], [x2(t)]])));

           d                           d
   deq := {-- x2(t) = x1(t) - 3 x2(t), -- x1(t) = x1(t) - 5 x2(t)}
           dt                          dt

> inits:={[0, 1, -2], [0, -2, 3], [0, 1.6, 2.7],
>      [0, -1.3, -2]};

   inits := {[0, -2, 3], [0, 1.6, 2.7], [0, -1.3, -2], [0, 1, -2]}

> DEplot(deq, [x1(t), x2(t)], t=-5..5, x1=-4..4, x2=-4..4, inits,
> scaling=constrained);

# T-437/5-2: For MATLAB input, via Line Print mode.
> evalm(A &* matrix(2, 1, [x1, x2]));

                             [x1 - 5 x2]
                             [         ]
                             [x1 - 3 x2]

> convert(convert(", vector), list);

                        [x1 - 5 x2, x1 - 3 x2]

# T-437/6-1: Here x=0 is an stable center (see Table 9.1.1, page T-436).
# Note the elliptical trajectories.
> A:=matrix(2, 2, [2, -5, 1, -2]); evals:=eigenvals(A);

                                 [2    -5]
                            A := [       ]
                                 [1    -2]


                            evals := I, -I

> deq:=equate(matrix([[diff(x1(t), t)], [diff(x2(t), t)]]),
> evalm(A &* matrix([[x1(t)], [x2(t)]])));

          d                             d
  deq := {-- x1(t) = 2 x1(t) - 5 x2(t), -- x2(t) = x1(t) - 2 x2(t)}
          dt                            dt

> inits:={[0, 1, -2], [0, -2, 3], [0, 1.6, 2.7],
>      [0, -1.3, -2]};

   inits := {[0, -2, 3], [0, 1.6, 2.7], [0, -1.3, -2], [0, 1, -2]}

> DEplot(deq, [x1(t), x2(t)], t=-5..5, x1=-10..10,x2=-10..10,inits);

# T-437/6-2: For MATLAB input, via Line Print mode.
> evalm(A &* matrix(2, 1, [x1, x2]));

                            [2 x1 - 5 x2]
                            [           ]
                            [ x1 - 2 x2 ]

> convert(convert(", vector), list);

                       [2 x1 - 5 x2, x1 - 2 x2]

# T-437/15-1: Here vector x = [-2, 1] is an (offset) asymptotically
# stable spiral point (see Table 9.1.1, page T-436).
> A:=matrix(2, 2, [-1, -1, 2, -1]); evals:=eigenvals(A);

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


                                   1/2          1/2
                  evals := -1 + I 2   , -1 - I 2

> b:=matrix(2, 1, [-1, 5]); x0:=linsolve(A, evalm(-b));

                                   [-1]
                              b := [  ]
                                   [ 5]


                                    [-2]
                              x0 := [  ]
                                    [ 1]

> deq:=equate(matrix([[diff(x1(t), t)], [diff(x2(t), t)]]),
>      evalm(A &* matrix(2, 1, [x1(t), x2(t)]) + b));

deq :=

     d                              d
    {-- x1(t) = -x1(t) - x2(t) - 1, -- x2(t) = 2 x1(t) - x2(t) + 5}
     dt                             dt

> inits:={[0, 1, -2], [0, -2, 3], [0, 1.6, 2.7],
>      [0, -1.3, -2]};

   inits := {[0, -2, 3], [0, 1.6, 2.7], [0, -1.3, -2], [0, 1, -2]}

> DEplot(deq, [x1(t), x2(t)], t=-5..5, x1=-4..4, x2=-4..4, inits,
> scaling=constrained);

# T-437/15-2: For MATLAB input, via Line Print mode.
> [-x1-x2-1, 2*x1-x2+5];

                    [-x1 - x2 - 1, 2 x1 - x2 + 5]

>