#
# Phase planes are graphed for eigenvalues which are
#   1 pos, 1 neg
#   complex
#   imaginary
#
--------------------------------------------------------------------------------
# 
# Art Belmonte
# Thu, 28/Mar/96 (REVISED)
# Math 308-509
# Section 9.1: The Phase Plane: Linear Systems
# 
--------------------------------------------------------------------------------
# 
> with(vec_calc): with(DEtools):
Warning: new definition for   D
Error, (in a) wrong number (or type) of parameters in function readlib

--------------------------------------------------------------------------------
# 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).
> A:=matrix(2, 2, [3, -2, 2, -2]); evals:=eigenvals(A);\
DEplot2(A, [x1, x2], -5..5,\
{[0, 1, -2], [0, -2, 3], [0, 1.6, 2.7], [0, -1.3, -2]},\
x1=-4..4, x2=-4..4, scaling=constrained);

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

                                 evals := 2, -1
--------------------------------------------------------------------------------
# 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]);\
x0:=matrix(2, 1, [1, -2]);\
x:=evalm(exponential(A*t) &* x0);\
xa:=unapply(convert(convert(evalm(x), vector),\
     list), t);\
plot([op(xa(t)), t=-5..5], 'x'=-4..4, 'y'=-4..4, scaling=constrained);

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

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

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

    xa := t -> [- 5/3 exp(- t) + 8/3 exp(2 t), 4/3 exp(2 t) - 10/3 exp(- t)]
--------------------------------------------------------------------------------
# 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);\
DEplot2(A, [x1, x2], -5..5,\
{[0, 1, -2], [0, -2, 3], [0, 1.6, 2.7], [0, -1.3, -2]},\
x1=-4..4, x2=-4..4, scaling=constrained);

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

                            evals := - 1 + I, - 1 - I
--------------------------------------------------------------------------------
# T-437/5-2: For MATLAB input, via Line Print mode.
> evalm(A &* matrix(2, 1, [x1, x2]));\
convert(convert(", vector), list);

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

                             [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);\
DEplot2(A, [x1, x2], -5..5,\
{[0, 1, -2], [0, -2, 3], [0, 1.6, 2.7], [0, -1.3, -2]},\
x1=-10..10, x2=-10..10, scaling=constrained);

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

                                 evals := I, - I
--------------------------------------------------------------------------------
# T-437/6-2: For MATLAB input, via Line Print mode.
> evalm(A &* matrix(2, 1, [x1, x2]));\
convert(convert(", vector), list);

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

                            [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);\
b:=matrix(2, 1, [-1, 5]); x0:=linsolve(A, -b);\
evalm(A &* matrix(2, 1, [x1, x2]) + b);\
derivs:=convert(convert(", vector), list);\
DEplot2(derivs, [x1, x2], -5..5,\
{[0, 1, -2], [0, -2, 3], [0, 1.6, 2.7], [0, -1.3, -2]},\
x1=-4..4, x2=-4..4, scaling=constrained);

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

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

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

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

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

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

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