MGM.09 Eigensense: Diagonalizable Matrices,
Matrix Exponential, Matrix Powers and Dynamical Systems
TUTORIALS

Needs["LinearAlgebra`MatrixManipulation`"];
Needs["Graphics`Colors`"];
Off[General::"spell"];
Off[General::"spell1"];
Off[ParametricPlot3D::"ppcom"];
Off[ParametricPlot::"ppcom"];
Off[Plot::"plnr"];
SetOptions[Limit, Analytic -> True];
$Pre = Chop;

CMView = {2.7, 1.6, 1.2};

Vector Drawer

\(\(\( (*\ \(\(:\)\(Date : \ \ Copyright\ 2001\)\), \ 
      Math\ Everywhere, \ \(\(Inc\)\(.\)\)\ 
*) \)\(
\)\( (*
Original\ Conception\ to\ create\ a\ vector\ drawer\ is\ due\ to\ Bill\ Davis . 
Current\ code\ is\ by\ Bruce\ \(\(Carpenter\)\(.\)\)
*) \)\(
\)\( (*\ \ \(\(\(:\)\(Drawbacks : \ \ Since\ the\ vectors\ are\ drawn\ out\ of\)\)\ 

      \ \ \ \ \ context\), \ 
      the\ shapes\ of\ the\ arrowheads\ are\ only\ 
\ \ \ \ \ correct\ with\ AspectRatio\ set\ to\ \(\(Automatic\)\(.\)\)
*) \)\)\)

BeginPackage["Vector3D`","Graphics`Colors`"];

Axes3D::usage = "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, y, and z running from -a/3 to a, and with axes labels b units beyond the tips of the axes.  Axes3D[a] is Axes3D[a, a/8].";

Perpend::usage = "Perpend[a] returns a unit vector which is perpendicular to the vector a.";

Arrow::usage = "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. Arrow[a, Tail -> tail] gives a vector from tail to a.";
	
Vector::usage = "Vector[a, b] gives a 2D or 3D vector from a to b.";

ArrowHead::usage = "ArrowHead[a, vec] produces an arrowhead with its tip placed at the point a, pointing in the direction of the vector vec.";

Tail::usage = "Tail -> point puts the tail of the vector at point.";

Norm::usage = "Norm[vec] returns the Euclidean length of vec.";

Normalize::usage = "Normalize[vec] returns a unit vector in the direction of vec, provided vec is nonzero.";

Aperture::usage = "Aperture is the ratio of the radius of the base to the length of the head of a vector.";

TipSize::usage = "TipSize specifies an absolute size for the inner tip length of the head of a vector.";

TipRatio::usage = "TipRatio is the ratio of the inner length of the arrowhead to the length of the shaft of a vector.";

HeadRatio::usage = "HeadRatio is the ratio of the outer length of the arrowhead to the length of the shaft of a vector.";

HeadSize::usage = "HeadSize specifies an absolute size for the head of a vector.";

ScaleFactor::usage = "ScaleFactor specifies the amount to scale a vector in length and may be either a number or a pure function.";

ZeroVectorPointSize::usage = "ZeroVectorPointSize is the size of the point used to represent a zero vector.";

Options[Vector] = Options[Arrow] = Options[ArrowHead] =
{HeadRatio -> 0.18, TipRatio -> 0.14, Aperture -> 0.3, ScaleFactor -> 1, HeadType -> Polygon, ShaftQ -> True, EdgesQ -> True, VectorColor -> RGBColor[0,0,1], ShaftWidth -> 0.005, ZeroVectorPointSize -> 0.01};
		
Begin["`Private`"];

SetOptions[ParametricPlot,AspectRatio -> Automatic];
SetOptions[Plot,AspectRatio -> Automatic];
SetOptions[Graphics,AspectRatio -> Automatic];

Axes3D[u_,v_] := 
      Graphics3D[{
            {Blue,Line[{{-u/3,0,0},{u,0,0}}]},
                    Text["x",{u + v,0,0}],
            {Blue,Line[{{0,-u/3,0},{0,u,0}}]},
                    Text["y",{0, u + v,0}],
            {Blue,Line[{{0,0,-u/3},{0,0,u}}]},
                    Text["z",{0,0,u + v}]}]; 
Axes3D[u_] := Axes3D[u,u/8];
				
Norm[a_?VectorQ] := N[Sqrt[a.a]];

Normalize[a_?VectorQ] := If[a.a == 0, a, a/Norm[a]];

Perpend[a_?VectorQ] := Normalize[{a[[2]],-a[[1]]}]/;Length[a]==2;

Perpend[a_?VectorQ]:=
    Normalize[If[a[[1]]==0,{1,0,0},{a[[2]],-a[[1]],0}]]/;
      Length[a]==3;
      
Base3D = Table[{Cos[2 k Pi/8.], Sin[2 k Pi/8.]},{k,0,8}]; 

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, tip},
   {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, TipRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]; tr = 0.75 hr];
   If[NumberQ[ts], tr = ts/Norm[diff]];
   Graphics[{If[sq,{vc,Thickness[sw],Line[{from,tip - tr diff}]},{}],
      With[{trans = ap hr Norm[diff] Perpend[diff]},
      {vc, Thickness[sw], ht[{tip,
                   (tip - hr diff + trans),
                   (tip - tr diff),
                   (tip - hr diff - trans),tip}]}]}]
] /; Length[from]==Length[to]==2

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, eq, sw, hs, hr, tr,ap, sf, vc, zvps, diff = to - from, tip},
   {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics3D[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]];
   tr = hr/2;
   Graphics3D[{If[sq,{vc,Thickness[sw],Line[{from,tip-tr diff}]},{}], If[eq, {}, EdgeForm[]], {SurfaceColor[vc],
      With[{temp = Perpend[diff]},
      Polygon[Append[#,tip]]& /@ Partition[(tip - hr diff + #)&/@((ap hr Norm[diff] Base3D) . {temp,Normalize[Cross[temp,diff]]}),2,1]]}}]
] /; Length[from]==Length[to]==3


Arrow[a_, Tail->b_, opts___] := Vector[b,a+b,opts]

Arrow[a_, opts___] := Vector[Table[0,{Length[a]}],a,opts]

ArrowHead[a_, b_, opts___] := Vector[a-b,a,ShaftQ->False,opts];

End[];

EndPackage[];
Null

Graphics`Colors`GosiaGreen = 
        RGBColor[0,  0.392187,  0];

ThreeAxes[u,v]

ThreeAxes[u_, v_] := 
        Graphics3D[{{Blue, Line[{{-u, 0, 0}, {u, 0, 0}}]}, 
            Text["x", {u + v, 0, 0}], {Blue, 
              Line[{{0, -u, 0}, {0, u, 0}}]}, 
            Text["y", {0, u + v, 0}], {Blue, 
              Line[{{0, 0, 0}, {0, 0, u}}]}, 
            Text["z", {0, 0, u + v}]}];
ThreeAxes[u_] := ThreeAxes[u, (u/8)];
ThreeAxes::"usage" = "ThreeAxes[a,b] makes a standard cartesian \
axis graphics object with x, y, and z running from -a to a, and with axis \
labels b units beyond the tips of the axes.  ThreeAxes[a] is \
ThreeAxes[a,a/8].";

Needs["Statistics`ContinuousDistributions`"];

Clear[matrixplucker];
matrixplucker[expr_, 
        vars_] := Simplify[
          Outer[D, {expr}, 
            vars]][[1]]

T.1) Continuous dynamical systems 
Systems of linear homogeneous differential equations {x'[t],y'[t]} = A.{x[t],y[t]} solved with the matrix exponential {x[t],y[t]} = (e^(A t)).starter

T.1.a.i) Coefficient matrix of a system of linear homogeneous differential equations

Take this system of linear homogeneous differential equations
		x'[t]  =  0.3 x[t]  -  0.6 y[t],
		y'[t]  =  0.1 x[t]  +  0.8 y[t] 
Come up with a matrix A (called the coefficient matrix) to write this system in matrix form
		{x'[t], y'[t]} = A . {x[t], y[t]}..

Look at the system of linear homogeneous differential equations
		x'[t]  =  0.3 x[t]  -  0.6 y[t],
		y'[t]  =  0.1 x[t]  +  0.8 y[t].
Now make this matrix:

\(A\  = \ {{0.3, \(-0.6\)}, {0.1, 0.8}}; 
MatrixForm[A]\)

See what you get when you hit this matrix on {x[t], y[t]}:

\(Clear[x, y, t]; 
	A . {x[t], y[t]}\)

Because
         x'[t]  =  0.3 x[t]  -  0.6 y[t],
         y'[t]  =  0.1 x[t]  +  0.8 y[t] ,
This tells you that:
         {x'[t], y'[t]} = A . {x[t], y[t]}.
Check:

ColumnForm[Thread[{x'[t], y'[t]}  ==  A . {x[t], y[t]}]]

This checks and you're out of here.

T.1.a.ii) Slamming out a solution of a system of linear homogeneous differential equations  
 {x'[t],y'[t]} = A.{x[t],y[t]} with the matrix exponential

Stay with the same setup as in part i)
Use the matrix exponential to slam out formulas {x[t], y[t]} for solutions of
                     x'[t]  =  0.3 x[t]  -  0.6 y[t],
	          y'[t]  =  0.1 x[t]  +  0.8 y[t]
	          x[0] =  xstarter 
                     y[0] =   ystarter
in terms of xstarter and ystarter.

The system is
          x'[t]  =  0.3 x[t]  -  0.6 y[t],
          y'[t]  =  0.1 x[t]  +  0.8 y[t]
The matrix form of this system is:
        {x'[t], y'[t]} = A . {x[t], y[t]} 
where the coefficient matrix A is:

A = (({{0.3, -0.6}, {0.1, 0.8}}));
MatrixForm[A]

Getting formulas for solutions {x[t],y[t]} of
                     {x'[t], y'[t]} = A . {x[t], y[t]} 
in terms of xstarter and ystarter , is really easy.
You just hit the matrix exponential (E^(A t)) on {xstarter, ystarter}
and put
          {x[t], y[t]} = (E^(A t)) . {xstarter, ystarter}:

Clear[x, y, xstarter, ystarter];
{x[t_], y[t_]} = MatrixExp[A t].{xstarter, ystarter}

Check whether
         {x'[t],y'[t]}=A . {x[t], y[t]}:

Expand[{x'[t], y'[t]}]

Expand[A . {x[t], y[t]}]

Good.
Check whether
          {x[0], y[0]} = {xstarter, ystarter}

{x[0], y[0]}

Everything checks.

T.1.a.iii) Why that worked

Explain why that worked. 
And while you're at it explain this:
When you go with a 2D diagonalizable coefficient matrix A and put
          {x'[t],y'[t]}=A . {x[t], y[t]}, 
then you can slam out formulas for the solution {x[t], y[t]} that starts at a given point
            {x[0], y[0]}  =  {a, b} 
by putting
            {x[t], y[t]} = (E^(A t)) . {a, b}.

Just check it:

Step 1: Check whether {x'[t],y'[t]}=A . {x[t], y[t]}:

Take the formula
          {x[t], y[t]} = (E^(A t)) . {a, b}
and differentiate with respect to t and see what you get:
          {x'[t],y'[t]}=\(D[{x[t], y[t]}, 
            t]\ 
			 = \ \(D[E\^\(A\ t\) . {a, b}, t]\ 
				 = \ 
            A . E\^\(A\ t\) . {a, b}\)\)

So {x[t], y[t]} = (E^(A t)) . {a, b} gives
           {x'[t],y'[t]} = A . (E^(A t)) . {a, b}
and this is the same as
           {x'[t],y'[t]}=A . {x[t], y[t]}
because
           {x[t], y[t]} = (E^(A t)) . {a, b}

Step 2: Check whether {x[0], y[0]} = {a, b}:

Take the formula
           {x[t], y[t]} = (E^(A t)) . {a, b}
and plug in t = 0 and see what you get:
{x[0], y[0]} =  
          (E^(A 0)) . {a, b}  =  
            IdentityMatrix . {a, b}  =  [a, b} 
 
The result: Everything checks.

T.1.b.i) Using eigenvalues of the coefficient matrix to predict how solutions plot out

Here's a 2D system of homogeneous linear differential equations:

Clear[x, y, t];
linearsystem={x'[t],y'[t]}=={-0.8 x[t] + 0.1 y[t], 0.9 x[t] - 1.75 y[t]};

ColumnForm[Thread[linearsystem]]

The coefficient matrix of this system is:

A=(({{-0.8, 0.1}, {0.9, -1.75}}));
MatrixForm[A]

The eigenvalues of A are:

Eigenvalues[A]

Use this information to predict how solutions {x[t],y[t]} of this linear system plot out as t gets large.
Illustrate your answer with a good plot.

Given a starting point {a, b}, the solution {x[t], y[t]} with {x[0], y[0]} = {a, b} is
                {x[t], y[t]} = (E^(A t)) . {a, b}.
Take another look at the eigenvalues of A::

Eigenvalues[A]

Both negative. 
According to the Basics, this tells you that no matter what {a, b} is, you are guaranteed that   
         {x[t], y[t]} = (E^(A t)) . {a, b} -> {0, 0} 
as t gets large..
The upshot: As t gets large, all solutions of this linear system are attracted to {0,0}.
See it happen for four solutions starting at random points:

Clear[trajectoryplot, eigenplot, scaler, eigenvector, pointer, 
          starter, thigh, k, j];

trajectoryplot[thigh_, starter_]  :=  
          ParametricPlot[MatrixExp[A t] . starter, {t, 0, thigh}, 
            PlotStyle -> {{Thickness[0.012], CadmiumOrange}}, 
            DisplayFunction -> Identity];


\(pointer[thigh_, starter_]\  := \ 	ArrowHead[
            MatrixExp[A\ \ thigh]\  . 
              starter, \((D[MatrixExp[A\ t], t] /. t -> thigh)\)\  . 
              starter, \ 
			\ \ \ \ \ \ \ \ HeadSize -> 0.6, 
            VectorColor -> Black, Aperture -> 0.4];\)

starterplots  =  
          Table[Graphics[{Navy, PointSize[0.03], Point[starter[j]]}], {j, 1, 
              4}]; 

starter[1] =  {Random[Real, {4, 10}], 
          Random[Real, {4, 10}]};
starter[
          2] =  {Random[Real, {-10, -4}], 
          Random[Real, {-10, -4}]};
starter[
          3] =  {Random[Real, {-10, -4}], 
          Random[Real, {4, 10}]};
starter[4] =  {Random[Real, {4, 10}], 
            Random[Real, {-10, -4}]};

Clear[solutionplots, thigh];
solutionplots[thigh_] := Show[Table[trajectoryplot[thigh, starter[k]], {k, 1, 4}] ,
	Table[pointer[thigh, starter[k]], {k, 1, 4}],starterplots,PlotRange -> All,PlotLabel->"Solutions for 0< t"<thigh,DisplayFunction -> Identity];
Show[solutionplots[0.2], 
        DisplayFunction -> $DisplayFunction];
Show[
        solutionplots[0.5], 
        DisplayFunction -> $DisplayFunction];
Show[
        solutionplots[2.0], 
        DisplayFunction -> $DisplayFunction];
Show[
        solutionplots[4.0], 
        DisplayFunction -> $DisplayFunction];
Show[
        solutionplots[10.0], 
        DisplayFunction -> $DisplayFunction];
Show[
        solutionplots[50.0], 
        DisplayFunction -> $DisplayFunction];

Sucked down the drain at {0,0}.
 And you knew this in advance because both eigenvalues of the coefficient matrix A are negative.

T.1.b.ii) Another linear system with complex eigenvalues

Here's a new 2D system of homogeneous linear differential equations:

Clear[x, y, t];
linearsystem={x'[t],y'[t]}=={0.8 x[t] + 1.7 y[t], -2.9 x[t] - 
                1.75 y[t]};

ColumnForm[
      Thread[linearsystem]]

The coefficient matrix of this system is:

A=(({{0.8, 1.7}, {-2.9, -1.75}}));
MatrixForm[A]

The eigenvalues of A are:

Eigenvalues[A]

Use this information to predict how solutions {x[t],y[t]} of this linear system plot out as t gets large.
Illustrate your answer with a good plot.

Given a starting point {a, b}, the solution {x[t], y[t]} with {x[0], y[0]} = {a, b} is
                {x[t], y[t]} = (E^(A t)) . {a, b}.
Take another look at the eigenvalues of A::

Eigenvalues[A]

This are of the form p + I q and p - I q with p < 0 and q ≠ 0. 
Because p < 0 this tells you that no matter what {a, b} is, you are guaranteed that   
         {x[t], y[t]} = (E^(A t)) . {a, b} -> {0, 0} 
as t gets large..
Because q ≠ 0 this tells you that no matter what {a, b} is, you are guaranteed that   
         {x[t], y[t]} = (E^(A t)) . {a, b} swirls its way into {0,0}
as t advances..
The upshot: As t gets large, all solutions of this linear system swirl their way to {0,0}.
See it happen for four solutions starting at random points:

Clear[trajectoryplot, eigenplot, scaler, pointer, eigenvector, 
          starter, thigh, k, j];

trajectoryplot[thigh_, starter_]  :=  
          ParametricPlot[MatrixExp[A t] . starter, {t, 0, thigh}, 
            PlotStyle -> {{Thickness[0.012], CadmiumOrange}}, 
            DisplayFunction -> Identity];


\(pointer[thigh_, starter_]\  := \ 	ArrowHead[
            MatrixExp[A\ \ thigh]\  . 
              starter, \((D[MatrixExp[A\ t], t] /. t -> thigh)\)\  . 
              starter, \ 
			\ \ \ \ \ \ \ \ HeadSize -> 1, 
            VectorColor -> Black, Aperture -> 0.4];\)

starterplots  =  
          Table[Graphics[{Navy, PointSize[0.03], Point[starter[j]]}], {j, 1, 
              4}]; 

starter[1] =  {Random[Real, {4, 10}], 
          Random[Real, {4, 10}]};
starter[
          2] =  {Random[Real, {-10, -4}], 
          Random[Real, {-10, -4}]};
starter[
          3] =  {Random[Real, {-10, -4}], 
          Random[Real, {4, 10}]};
starter[4] =  {Random[Real, {4, 10}], 
            Random[Real, {-10, -4}]};

Clear[solutionplots, thigh];
solutionplots[thigh_] := Show[Table[trajectoryplot[thigh, starter[k]], {k, 1, 4}] ,
	Table[pointer[thigh, starter[k]], {k, 1, 4}],starterplots,PlotRange -> All,PlotLabel->"Solutions for 0< t"<thigh,DisplayFunction -> Identity];
ranger  =  13;
Show[solutionplots[1.0], 
          PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
          DisplayFunction -> $DisplayFunction];


Show[solutionplots[2.0], 
          PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
          DisplayFunction -> $DisplayFunction];

Show[solutionplots[7.0], 
          PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
          DisplayFunction -> $DisplayFunction];

Swirling down the drain at {0,0}.
And you knew this in advance because the eigenvalues of the coefficient matrix A are of the form p + I q and p - I q with p < 0 and q ≠ 0.

T.2) Discrete dynamical systems: Customer dynamics solved with matrix powers

T.2.a.i) Setting up a the discrete difference equations

Southern Ag Equipment, importers of specialized carrot harvesters, has been through a lengthy process of computerizing their sales records through the Macintosh software QuickBook Pro from Intuit.

Their customer base consists of the 1500 active carrot farmers in the US.
Through computer analysis, they have noticed: 
▪ 71% of their customers who order in year k also order in year k + 1
▪ 45% of their customers who did not order in year k do order in year k + 1.
Put
          x[k] = fraction of their customers who do order in year k 
          y[k] =  fraction of their customers who do not order in year k. 
Explain why this gives
         x[k + 1] = 0.71 x[k] + 0.45 y[k] 
 and
         y[k + 1] = 0.29 x[k] + 0.55 y[k].

The relationship
            x[k + 1] = 0.71 x[k] + 0.45 y[k]  
makes sense because 71% of  their customers who ordered in year k also order in year k+1 and  
45% of their customers who did not order in year k do order in year k + 1.
The relationship
             y[k + 1] = 0.29 x[k] + 0.55 y[k].
makes sense because 29% of  their customers who ordered in year k also do not order in year k+1 and  
55% of their customers who did not order in year k do not order in year k + 1.

T.2.a.ii) Solving the discrete difference equations with matrix powers

Stay with everything just as it was in part i) above.
Use the relationships
         x[k + 1] = 0.71 x[k] + 0.45 y[k] and
         y[k + 1] = 0.29 x[k] + 0.55 y[k]

to explain why

{x[k], y[k]}=  .(((({{0.71`, 0.45`}, {0.29`, 0.55`}})))^k).{x[0],y[0]}

The relationships
         x[k + 1]=0.71 x[k]+0.45 y[k] and
         y[k + 1] = 0.29 x[k] + 0.55 y[k]
are the same as
       {x[k + 1], y[k + 1]}=(({{0.71, 0.45`}, {0.29, 0.55`}})).{x[k], y[k]}.
To see why just multiply it out:

Clear[x, y, k];
(({{0.71, 0.45}, {0.29, 0.55}})).{x[k], y[k]}

This gives:

{x[1], y[1]}=(({{0.71, 0.45`}, {0.29, 0.55`}})).{x[0], y[0]}.


         {x[2], y[2]}=(({{0.71, 0.45`}, {0.29, 0.55`}})).{x[1], y[1]} 

                            = (({{0.71, 0.45`}, {0.29, 0.55`}})).(({{0.71, 0.45`}, {0.29, 0.55`}})).{x[0], y[0]}
                
                            = (((({{0.71, 0.45`}, {0.29, 0.55`}})))^2).{x[0], y[0]}.
                               
                 
         {x[3], y[3]}=(({{0.71, 0.45`}, {0.29, 0.55`}})).{x[2], y[2]} 

                          = (({{0.71, 0.45`}, {0.29, 0.55`}})).(((({{0.71, 0.45`}, {0.29, 0.55`}})))^(2  )).{x[0], y[0]}
                 
                          = (((({{0.71, 0.45`}, {0.29, 0.55`}})))^3).{x[0], y[0]}
                
This leads to the general pattern:
            {x[k], y[k]} =  (((({{0.71, 0.45`}, {0.29, 0.55`}})))^k).{x[0], y[0]}.

T.2.a.iii) The ultimate ratio

Stay with the same set up as in part i) and remember that
      x[k] = fraction of their customers who order in year k 
      y[k] =  fraction of their customers who do not order in year k.  
Predict the ultimate ratio x[k]/y[k] as k gets large and interpret the result.

You know
      {x[k], y[k]} = (A^k) . {x[0], y[0]}.
with 
      A= (({{0.71, 0.45}, {0.29, 0.55}})).
Enter A and check its eigenvalues and eigenvectors:

A =   (({{0.71, 0.45}, {0.29, 0.55}}));
Clear[eigenvalue];
{eigenvalue[1], eigenvalue[2]}  =  
          Eigenvalues[A]

Note that | eigenvalue[1] |  > | eigenvalue[2] | .
Check the eigenvectors:

Clear[eigenvector];
{eigenvector[1], eigenvector[2]}  =  Eigenvectors[A]

Remember from the Basics:
 If | eigenvalue[1] |  > | eigenvalue[2] | , then the ultimate direction of (A^k) . X as k gets large is in the direction of eigenvector[1] (unless X is a multiple of eigenvector[2]).
Here 
         X = {x[0], y[0]}
cannot be a multiple of eigenvector[2] because {x[0], y[0]} cannot have negative entries.
So  the ultimate ratio x[k]/y[k] as k gets large is approximately:

UltimateRatio  =  0.840571/0.541701

In the long run,Southern Ag will have about 1.55 customers who order for each customer who doesn't.

T.2.a.iv) Limiting values

Keep everything the same as in the parts above. Give the limiting values of x[k] and y[k] and interpret the result.

The point to remember here is that
      x[k] = fraction of their customers who order in year k 
      y[k] =  fraction of their customers who do not order in year k. 
So 
          x[k] + y[k] = 1.
 You also know that as k gets large  x[k]/y[k] closes in on:

UltimateRatio

So in the long run, 
          x[k] ≈  UltimateRatio  y[k].
Plug this into
           x[k] + y[k] = 1
 to get
          (UltimateRatio  y[k]) + y[k] ≈ 1.
 This is the same as
          (UltimateRatio   +  1) y[k] ≈ 1.
 This tells you that in the long run (for k large),  y[k] is approximately:

ylongrun  =  1/(1  +  UltimateRatio)

And this tells you that in the long run (for k large),  x[k] is approximately:

xlongrun  =  1  -  ylongrun

The interpretation: 
In the long run, about 61% of Southern Ag's customers will order in a given year and about 39% will not.

See it happen by starting with x[0] = 0 and y[0] = 1, the situation when Southern Ag opened for business:

A  =  {{0.71, 0.45}, {0.29, 0.55}};
starter  =  {0, 1};
Clear[x, y, k];
ColumnForm[
      Table[{{x[k], y[k]}, "=", MatrixPower[A, k] . {0, 1}}, {k, 0, 
          10}]]

Compare:

{xlongrun, ylongrun}

The model predicts that after about 4 years, Southern Ag's  long term market penetration will be approximately achieved.

T.3) Discrete dynamical systems: Linear homogeneous difference equations solved with matrix powers

T.3.a.i) Linear difference equation in matrix form

Take the linear homogeneous difference equation
		y[k + 2] - 0.5 y[k + 1] + 0.27 y[k] = 0, 
Come up with a matrix A (called the coefficient matrix) to write this difference equation in matrix form
		{y[k], y[k + 1]} = (A^k) . {y[0], y[1]}
for all k = 0, 1, 2, 3, 4, 5,  .   .   . .

Look at the difference equation:
        y[k + 2] = 0.5 y[k + 1] - 0.27 y[k];
Now make this matrix:

A = (({{0, 1}, {0.5, -0.27}}));
MatrixForm[
      A]

See what you get when you hit this matrix on {y[k], y[k + 1]}:

Clear[y, k];
A.{y[k], y[k + 1]}

Because
          y[k + 2] = 0.5 y[k + 1] - 0.27 y[k],
This tells you that:
         {y[k + 1], y[k + 2]} = A . {y[k], y[k + 1]}
for all k's.

Plug in k = 0 to get
             {y[1], y[2]} = A . {y[0], y[1]}.
             
Plug in k = 1 into
            {y[k + 1], y[k + 2]} = A . {y[k], y[k + 1]}
to get
              {y[2], y[3]} = A . {y[1], y[2]} = A . A . {y[0], y[1]}

Plug in k = 2 into
            {y[k + 1], y[k + 2]} = A . {y[k], y[k + 1]}
to get
              {y[3], y[4]} = A . {y[2], y[3]} = A . A . A . {y[0], y[1]}

Continuing on in this way, you are led to the matrix form
             {y[k], y[k + 1]} = (A^k) . {y[0], y[1]}
of the given difference equation.

T.3.a.ii) Using the matrix to slamming out y[2], y[3], y[4], y[5], y[6], y[7], y[8], y[9] and y[10] in terms of y[0] and y[1]

Stay with the same setup as in part i)
Use the matrix form  of the difference equation to slam out  
           y[2], y[3], y[4], y[5], y[6], y[7], y[8], y[9] and y[10]
in terms of y[0] and y[1].

The difference equation is
        y[k + 2] = 0.5 y[k + 1] - 0.27 y[k];
The matrix form is:
        {y[k], y[k + 1]} = (A^k) . {y[0], y[1]}where A is:

MatrixForm[A]

Use the matrix form
         {y[k], y[k + 1]} = (A^k) . {y[0], y[1]}
with k = 1 to see that
          {y[1], y[2]} 
is given by:

\(\(\ k\  = \ 1; 
Clear[y]; 
MatrixPower[A, k] . {y[0], y[1]}\)\)

Use the matrix form
         {y[k], y[k + 1]} = (A^k) . {y[0], y[1]}
with k = 2 to see that
          {y[2], y[3]} 
is given by:

\(\(\ k\  = \ 2; 
MatrixPower[A, k] . {y[0], y[1]}\)\)

Here comes {y[3], y[4]} in terms of y[0] and y[1]:

\(\(\ k\  = \ 3; 
MatrixPower[A, k] . {y[0], y[1]}\)\)

Here comes {y[4], y[5]} in terms of y[0] and y[1]:

\(\(\ k\  = \ 4; 
MatrixPower[A, k] . {y[0], y[1]}\)\)

Here comes {y[5], y[6]} in terms of y[0] and y[1]:

k  =  5;
MatrixPower[A, k] . {y[0], y[1]}

Here comes {y[6], y[7]} in terms of y[0] and y[1]:

k  =  6;
MatrixPower[A, k] . {y[0], y[1]}

Here comes {y[7], y[8]} in terms of y[0] and y[1]:

k  =  7;
MatrixPower[A, k] . {y[0], y[1]}

Here comes {y[8], y[9]} in terms of y[0] and y[1]:

k  =  8;
MatrixPower[A, k] . {y[0], y[1]}

Here comes {y[9], y[10]} in terms of y[0] and y[1]:

k  =  9;
MatrixPower[A, k] . {y[0], y[1]}

If you want a list, you can get it

Clear[k];
\(	ColumnForm[
      Table[{y[k]\ , \((MatrixPower[A, k] . {y[0], y[1]})\)[\([1]\)]}, {k, 2, 
          20}]]\)

T.3.a.iii) Using eigenvalues to size up the situation

Stay with the same setup as in part i) and look at the coefficient matrix A for the
difference equation y[k + 2] - 0.5 y[k + 1] + 0.27 y[k] = 0:

MatrixForm[A]

Check the eigenvalues of A:

Eigenvalues[A]

Use this eigenvalue information to determine what happens to y[k] as k gets large.

Look at the eigenvalues again:

Eigenvalues[A]

\(The\ absolute\ values\ of\ both\ eigenvalues\ \ are\ less\ than\ 1. \)
\(One\ of\ the\ Basics\ \(says : 
	
              If\ A\ is\ any\ diagonalizable\ matrix\ and\ the\ absolute\ values\ of\ the\ of\ all\ the\ \ \ \ \ \ \ \ \ \ \ eigenvalues\ of\ A\ are\ less\ than\ 1\), \ then\ no\ matter\ what\ X\ \ is, \ 
          you\ are\ guaranteed\ that\)

(A^k) X -> {0, 0}.

So no matter what {y[0], y[1]} are,  
      you are guaranteed that
              {y[k], y[k + 1]}=(A^k) . {y[0], y[1]} ->{0, 0} as k gets large.

The upshot:
             y[k] -> 0 as k gets large.

See it happen:

Clear[k];
\(	ColumnForm[
      Table[{y[k]\ , \((MatrixPower[A, k] . {y[0], y[1]})\)[\([1]\)]}, {k, 
          10, 100, 10}]]\)

Withering away to 0.

T.3.b.i) What is a linear homogeneous difference equation?

What is a linear homogeneous difference equation?

A  linear homogeneous difference equation is any difference equation that can be put in matrix form as above.

T.3.b.ii) The message

What's the message?

Once you have put a linear homogeneous difference equation in matrix form, you can analyze its behavior in by knowing how the matrix power function plots out. This is what was done above. It you want to see some more, go on.

T.3.c.i) Another one

Use the matrix form of
		y[k + 2] - 1.4 y[k + 1] + 0.3  y[k] = 0
to study the behavior of solutions y[k] as k->infinity

Look at the difference equation:
        y[k + 2] = 1.4 y[k + 1] - 0.3 y[k];
Make the coefficient matrix:

A = (({{0, 1}, {1.4, -0.3}}));
MatrixForm[A]

The matrix form of
        y[k + 2] - 1.4 y[k + 1] + 0.3  y[k] = 0
is 
       {y[k], y[k + 1]} = (A^k) . {y[0], y[1]}.
Check the eigenvalues of A

Clear[eigval, eigvect];
\(	{eigval[1], eigval[2]}\  = \ Eigenvalues[A]\)

Now check the eigenvectors:

{eigvect[1], eigvect[2]}  =  Eigenvectors[A]

The juicy info:
            | eigval[1] |  > 1 and 
	| eigval[1] |  > | eigval[2] | 
This tells you that unless {y[0], y[1]} is a multiple of eigenvector[2], you are guaranteed that as k gets large,
         
       || {y[k], y[k + 1]} ||  = 
        (( || A)^k) . {y[0], y[1]} ||   ->  infinity.
And this tells that unless {y[0], y[1]} is a multiple of eigvect[2], you are guaranteed that as k gets large,
            {y[k], y[k + 1]}
gravitates to the line through {0,0} defined by eigvect[1]; so that as k gets large
            y[k + 1]/y[k] -> (0.782245/(-0.622971)):

-(0.782245/0.622971)

On the other hand, because | eigval[2] |  < 1:

eigval[2]

This tells you that if {y[0], y[1]} is a multiple of eigvect[2], you are guaranteed that as k gets large,
         
       || {y[k], y[k + 1]} ||  = 
        (( || A)^k) . {y[0], y[1]} ||   ->  {0, 0}.
And this tells you that k gets large, 
            {y[k], y[k + 1]}
gravitates to the line through {0,0} defined by eigvect2].:

eigvect[2]

This tells you that as k gets large
            y[k + 1]/y[k] -> (0.6909/0.72295):

(0.6909/0.72295)

T.3.c.ii) Yet another one

Use the matrix form of
		y[k + 2] - 1.2 y[k + 1] + 0.3  y[k] = 0
to study the behavior of solutions y[k] as k->infinity

Look at the difference equation:
        y[k + 2] = 1.2 y[k + 1] - 0.3 y[k];
Make the coefficient matrix:

A = (({{0, 1}, {1.2, -0.3}}));
MatrixForm[A]

The matrix form of
        y[k + 2] - 1.4 y[k + 1] + 0.3  y[k] = 0
is 
       {y[k], y[k + 1]} = (A^k) . {y[0], y[1]}.
Check the eigenvalues of A

Clear[eigval, eigvect];
\(	{eigval[1], eigval[2]}\  = \ Eigenvalues[A]\)

Now check the eigenvectors:

{eigvect[1], eigvect[2]}  =  Eigenvectors[A]

The juicy info:
            | eigval[1] |  > 1 and 
	| eigval[1] |  > | eigval[2] | 
This tells you that unless {y[0], y[1]} is a multiple of eigenvector[2], you are guaranteed that as k gets large,
         
       || {y[k], y[k + 1]} ||  = 
        (( || A)^k) . {y[0], y[1]} ||   ->  infinity.
And this tells that unless {y[0], y[1]} is a multiple of eigvect[2], you are guaranteed that as k gets large,
            {y[k], y[k + 1]}
gravitates to the line through {0,0} defined by eigvect[1]; so that as k gets large
            y[k + 1]/y[k] -> (0.782245/(-0.622971)):

-(0.782245/0.622971)

On the other hand, because | eigval[2] |  < 1:

eigval[2]

This tells you that if {y[0], y[1]} is a multiple of eigvect[2], you are guaranteed that as k gets large,
         
       || {y[k], y[k + 1]} ||  = 
        (( || A)^k) . {y[0], y[1]} ||   ->  {0, 0}.
And this tells you that k gets large, 
            {y[k], y[k + 1]}
gravitates to the line through {0,0} defined by eigvect2].:

eigvect[2]

This tells you that as k gets large
            y[k + 1]/y[k] -> (0.6909/0.72295):

(0.6909/0.72295)

T.4)  How do you pronounce the words "eigenvector" and "eigenvalue"?  
What about the algebra behind calculating them?

T.4.a) How do you pronounce the words "eigenvector" and "eigenvalue"?

How do you pronounce the words "eigenvector" and "eigenvalue"?

The words "eigenvector" and "eigenvalue" are bastardized words. "Eigen" comes from German, "vector"  and "value" come from English. You pronounce "eigen" as if it were German and you pronounce "vector" and "value" as usual.
"Eigen" is pronounced I - gen with a hard g as in "garden" or "geezer," not the soft g as in "general" or "genetics." 
The German word "eigen" corresponds to the English word "own."

T.4.b.i)  Using roots of the characteristic polynomial Det[A - λ Identity] to calculate eigenvalues of a given square matrix A

Go with this sample matrix A:

dim  =  4;
A  =  
        Table[Random[Real, {-5, 5}], {i, 1, dim}, {j, 1, 
            dim}];
MatrixForm[A]

Put characteristicpoly[λ] = Det[A - λ IdentityMatrix]:

Clear[
        characteristicpoly, λ];
characteristicpoly[λ_] =Det[A - λ IdentityMatrix[dim]]

Find when characteristicpoly[λ] is 0:

Solve[characteristicpoly[λ]==0,λ]

Compare:

Eigenvalues[A]

Bingo.
The solutions λ of
                  characteristicpoly[λ] = 0 
 are the eigenvalues of A.
Explain why the same thing will work for any square matrix A.

Saying that λ is an eigenvalue for a given matrix A is the same as saying that there is a vector X ≠ {0,0,. . . ,0} with
              A.X = λ X.
              
This is the same as saying that there is a vector X ≠ {0,0,. . . ,0} with
              A.X - λ X. = {0,0,. . . ,0}.
              
And this is the same as saying that there is a vector X ≠ {0,0,. . . ,0} with 
          (A - λ Identity).X = {0,0,. . . ,0}. 
. 
This is the same as saying that the matrix  (A - λ Identity) hits a non-zero vector into {0,0,...,0}.

This is the same as saying that the  matrix (A - λ Identity) has at least one zero SVD stretch factor.

This is the same as saying that 
              characteristicpoly[λ] = Det[ (A - λ Identity)] = 0. 
              
Try it out on more random square matrices A:

dim  =  Random[Integer, {2, 10}];
A  =  
        Table[Random[Real, {-5, 5}], {i, 1, dim}, {j, 1, 
            dim}];
MatrixForm[A]

Clear[characteristicpoly, λ];
characteristicpoly[λ_]  = 
      Det[A  -  λ IdentityMatrix[dim]]

Solve[characteristicpoly[λ] == 0, λ]

Eigenvalues[A]

Rerun all four cells until you are satisfied.

T.4.b.ii)  Using roots of the characteristic polynomial Det[A - λ Identity] to calculate eigenvectors of a given square matrix A

Go with this sample matrix A:

dim  =  4;
A  =  
        Table[Random[Real, {-5, 5}], {i, 1, dim}, {j, 1, 
            dim}];
MatrixForm[A]

Put characteristicpoly[λ] = Det[A - λ IdentityMatrix]:

Clear[
        characteristicpoly, λ];
characteristicpoly[λ_] =Det[A - λ IdentityMatrix[dim]]

The roots of the characteristic polynomial are the eigenvalues of A:

Clear[eigenvalue, k];
eigenvalue[k_] :=Eigenvalues[A][[k]];
Table[eigenvalue[k],{k, 1, dim}]

How can you use the calculated eigenvalues to come up with the corresponding eigenvectors?

If you have your hands on eigenvalue[k], then the corresponding eigenvectors X are solutions of 
               A.X = eigenvalue[k] X.
               
This is the same as saying that
          A.X - eigenvalue[k] X = {0,0,.. . ,0}.
 
This is the same as saying that
           (A - eigenvalue[k] Identity).X = {0,0,. . .,0}.
 
This is the same as saying that
          X is in the null space 
  of the matrix   (A - eigenvalue[k] Identity) 
          .
Here is the unit eigenvector corresponding to eigenvalue[1]:

eigenvector1 =Flatten[NullSpace[A-eigenvalue[1] IdentityMatrix[
                dim]]]

Check to be sure that  A.eigenvector1 = eigenvalue[1] eigenvector1:

A.eigenvector1

eigenvalue[1] eigenvector1

Here are the unit eigenvectors corresponding to each of the eigenvalues:

Clear[eigenvector, k];
eigenvector[k_] :=Flatten[NullSpace[A-eigenvalue[k] IdentityMatrix[dim]]];
ColumnForm[Table[eigenvector[k],{k, 1, dim}]]

Check them out:

Table[A.eigenvector[k] == eigenvalue[k] eigenvector[k],{k, 1, dim}]

Try it for more random square matrices A:

dim  = 6;
A  =  
        Table[Random[Real, {-5, 5}], {i, 1, dim}, {j, 1, 
            dim}];
MatrixForm[A]

Clear[eigenvalue, k];
eigenvalue[k_] :=Eigenvalues[A][[k]];

Clear[eigenvector, k];
eigenvector[k_] :=Flatten[NullSpace[A-eigenvalue[k] IdentityMatrix[dim]]];
ColumnForm[Table[eigenvector[k],{k, 1, dim}]]

ColumnForm[Eigenvectors[A]]

T.4.b.iii) Is  the route of characteristic polynomials and null spaces a recommended way of calculating eigenvalues and eigenvectors?

Is the route of characteristic polynomials and null spaces a recommended way of calculating eigenvalues and eigenvectors in daily life?

Heavens to Betsy - NO!

All good computer based mathematics processors such as Mathematica  or MatLab have professionally written built-in instructions that are optimized for calculation of eigenvalues and eigenvectors. 
Using these instructions is the recommended way of calculating eigenvalues and eigenvectors in daily life.

The route of characteristic polynomials and null spaces is the path of necessity for those unfortunate folks condemned to a life of pencil and paper hand calculation. 

You should know about route of characteristic polynomials and null spaces because you might meet one of these unfortunate folks out on the street or in a coffee house sometime.

Take pity on them. These folks are forced to spend hour after hour in pure useless calculation.

T.4.b.iv) Should you worry about calculation of eigenvectors and eigenvalues by hand?

Should you worry about calculation of eigenvectors and eigenvalues by hand?

Only if you enjoy self-inflicted pain. 
In 2D,the cheat sheet chart problem following this one gives as easy a hand path as there is.

T.5)  Traces, determinants and the 2D matrix exponential cheat sheet chart

T.5.a.i) Eigenvalues of  a 2D matrix A are
            (1/2) (trace + (Sqrt[(trace^2) - 4 det]))   and (1/2) (trace - (Sqrt[(trace^2) - 4 det]))

Here comes a cleared matrix:

Clear[a, b, c, d];
A=(({{a, b}, {c, d}}));
MatrixForm[A]

Now calculate the eigenvalues of A:

Eigenvalues[A]

Two quantities stand out:
           (a  +  d) and a d  -  b c
Folks have special names for these quantities.  They say
           trace of A  =  a  +  d
and 
           determinant of A  =  a d  -  b c.

{trace, det}  = {a + d, Det[A]}

You can easily express the eigenvalues of A in terms of the trace of A and the determinant of A:

ExpandAll[Eigenvalues[A]]==ExpandAll[{(1/2) (trace - (Sqrt[(trace^2) - 4 det])),(1/2) (trace + (Sqrt[(trace^2) - 4 det]))}]

So the two eigenvalues of A are given by the formulas:
        (1/2) (trace + (Sqrt[(trace^2) - 4 det])) 
and 
         (1/2) (trace - (Sqrt[(trace^2) - 4 det])).

These formulas make it possible for you to calculate the eigenvalues of a coefficient matrix with pencil and paper or with a cheap pocket calculator.

Lots of reference books have charts, of the type you will see below, that help you to see how the matrix exponential
           (E^(A t)) . starter 
 plots out once you know the trace and the determinant.
In this problem, you will participate in building one of these charts. To start building the chart, you go with trace and determinant axes and then you plot the curve
              (trace^2) - 4 determinant = 0:

cutoff = 
        Plot[((trace^2)/4), {trace, -5, 5}, 
          PlotStyle -> {{Red, Thickness[0.01]}}, 
          AxesLabel -> {"trace", "det"}, 
          PlotRange -> {-6, 6}, 
          DisplayFunction -> Identity];
\(\(cutofflabel = 
        Graphics[
          Text[FontForm["\!\(trace\^2\) - 4 det == 0", {"Times", 
                10}], {3.0, \(-4.0\)}]];\)
\)
\(\(pointer = 
        Arrow[{2, 1} - {3.0, \(-3.6\)}, Tail -> {3.0, \(-3.6\)}, 
          VectorColor -> Black];\)
\)
Show[cutoff, cutofflabel, pointer, 
        DisplayFunction -> $DisplayFunction];

And then you annotate this plot as follows:

twosuckercutoff = 
        Graphics[{Blue, Thickness[0.01], 
            Line[{{-5, 0}, {0, 0}}]}];
\(\(twosuckerlabel = 
        Graphics[
          Text[FontForm[" Both real and 
 both negative", {"Times", 
                10}], {\(-4\), 1}]];\)
\)
twopropcutoff = 
        Graphics[{Blue, Thickness[0.01], 
            Line[{{0, 0}, {5, 0}}]}];
\(\(twoproplabel = 
        Graphics[
          Text[FontForm[" Both real and 
 both positive", {"Times", 
                10}], {4, 1}]];\)
\)
propswirlcutoff = 
        Graphics[{CadmiumOrange, Thickness[0.01], 
            Line[{{0, 0}, {0, 6}}]}];
\(\(propswirllabel = 
        Graphics[
          Text[FontForm[" p + I q 
 and p - I q
 with p > 0 
 and q ≠ 0", {"Times", 10}], {2, 4}]];\)
\)
\(\(pureswirlerlabel = 
        Graphics[
          Text[FontForm[" p
 u
 r ", {"Times", 12}], {0, 
              7}]];\)
\)
suckswirlcutoff = 
        Graphics[{CadmiumOrange, Thickness[0.01], 
            Line[{{0, 0}, {0, 6}}]}];
\(\(suckswirllabel = 
        Graphics[
          Text[FontForm["p + I q 
 and p - I q
 with p < 0 
 and q ≠ 0", {"Times", 10}], {\(-2\), 4}]];\)
\)
negdetcutoff = 
        Graphics[{Blue, Thickness[0.01], 
            Line[{{-5, 0}, {5, 0}}]}];
\(\(negdetlabel = 
        Graphics[
          Text[FontForm["Both real and 
 one > 0 and one < 0", {"Times", 10}], {0, \(-1.5\)}]];\)
\)
\(\(pureswirlerlabel = 
        Graphics[
          Text[FontForm["0 + I q and 0 - I q on positive vertical axis", {"Times", 10}], {0, 7}]];\)
\)
chart = 
        Show[cutoff, cutofflabel, pointer, twosuckercutoff, twosuckerlabel, 
          twopropcutoff, twoproplabel, propswirlcutoff, propswirllabel, 
          suckswirlcutoff, suckswirllabel, negdetcutoff, negdetlabel, 
          pureswirlerlabel, 
          DisplayFunction -> $DisplayFunction];

Some folks even carry a copy of this chart in their wallets.
How do folks use this chart to predict how (E^(A t)) . starter plots out?

Here it is for a sample matrix A

a = -1.4;
b = 0.9;
c = 
        0.4;
d = -0.5;
A=(({{a, b}, {c, d}}));

MatrixForm[A]

You calculate 
        trace = a + d
and 
        det = a d - b c
by hand. And then you plot the resulting point 
       {trace, det}
on the chart:

trace = a + d;
det = a d - b c;
\(\(chartpointpoint = 
        Graphics[{Magenta, PointSize[0.04], 
            Point[{trace, det}]}];\)
\)
Show[chart, chartpointpoint];

Both eigenvalues are negative. This tells you that no matter what  {a,b} is, you can be certain that
(E^(A t)) . starter -> {0, 0} as t gets large.

Another:

a = -1.4;
b = -3.0;
c = 
        2.1;
d = 1.4;
A=(({{a, b}, {c, d}}));
MatrixForm[A]

trace = a + d;
det = a d - b c;
chartpointpoint = 
          Graphics[{HotPink, PointSize[0.04], Point[{trace, det}]}];

Show[chart, chartpointpoint];

This point {trace,det}:

{trace, det}

lands on the positive vertical axis. 
This tells you that as t advances, (E^(A t)) . starter oscillates on an ellipse centered at {0,0}.

Here is the same thing for some random matrices:

a = Random[Real, {-2, 2}];
b = 
        Random[Real, {-2, 2}];
c = 
        Random[Real, {-2, 2}];
d = 
        Random[Real, {-2, 2}];
A=(({{a, b}, {c, d}}));
MatrixForm[A]

trace = a + d;
det = a d - b c;
chartpointpoint = 
          Graphics[{HotPink, PointSize[0.04], Point[{trace, det}]}];

Show[chart, chartpointpoint];

Neat cheat sheet.
Print one up and put it in your wallet.

T.5.a.ii) Explanations

Explain why the chart is correct.

p + I q and p - I q with p < 0 and q ≠ 0 explanation

Here are some sample points {trace, det} that plot out in the two swirling sucker region of the chart:

Clear[trace, k];
\(	\(trace[k_]\  := \ \(trace[k] = \ 
          Random[Real, {\(-4\), \(-0.1\)}]\);\)\)
samplepoints  =  
        Table[{trace[k], (trace[k]^2)/4  +  
              Random[Real, {0, 5  -  (trace[k]^2)/4}]}, {k, 1, 
            100}];
samplepointplot  =  
        ListPlot[samplepoints, PlotStyle -> {Banana, PointSize[0.03]}, 
          DisplayFunction -> Identity];
\(\(Show[samplepointplot, chart, PlotRange -> All, 
        AxesLabel -> {"trace", "det"}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> $DisplayFunction];\)\)

The p + I q and p - I q with p < 0 and q ≠ 0 region corresponds to the case that
         trace of A < 0,
         determinant of A > 0, and
        (((trace of A))^2) < 4 determinant of A.
Explanation:
The eigenvalues of A are
          p + I q = (1/2) (trace + (Sqrt[(trace^2) - 4 det]))
and
          p - I q = (1/2) (trace - (Sqrt[(trace^2) - 4 det])).
Because 
          det > 0
and
          (trace^2) - 4 det < 0, 
you see q is not 0. This gives you the whirl. 
Because trace < 0, you see that p < 0. This gives you the suck.

0 + I q and 0 - I q explanation:

Here are some sample points {trace, det} that plot out in 0 + I q and 0 - I q  region of the chart:

samplepoints  =  
        Table[{0, Random[Real, {0, 6}]}, {k, 1, 10}];
samplepointplot  =  
        ListPlot[samplepoints, PlotStyle -> {CadmiumOrange, PointSize[0.04]}, 
          DisplayFunction -> Identity];
\(\(Show[samplepointplot, chart, PlotRange -> All, 
        AxesLabel -> {"trace", "det"}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> $DisplayFunction];\)\)

The  0 + I q and 0 - I q region corresponds to the case that
         trace of A = 0,
and 
         determinant of A > 0.
Explanation:
The eigenvalues of A are   
         p + I q = (1/2) (trace + (Sqrt[(trace^2) - 4 det])) 
and
         p - I q = (1/2) (trace - (Sqrt[(trace^2) - 4 det])).

Because 
        trace = 0,
        p + I q = (1/2) (0 + (Sqrt[0 - 4 det]))
and
        p - I q = (1/2) (0 - (Sqrt[0 - 4 det])).
And because det  >  0,
        p + I q = (1/2) (0 + 2 I (Sqrt[Abs])[det]) 
and
        p - I q = (1/2) (0 - 2 I (Sqrt[Abs])[det])

p + I q and p - I q with p > 0 and q ≠ 0  explanation:

Here are some sample points {trace, det} that plot out in the p + I q and p - I q with p > 0 and q ≠ 0 region of the chart:

Clear[trace, k];
\(	\(trace[k_]\  := \ \(trace[k] = \ 
          Random[Real, {0.1, 4}]\);\)\)
\(\(samplepoints = 
        Table[\ {trace[k], \((trace[k]^2)\)/4\  + \ 
              Random[Real, {0, 
                  5\  - \ \((trace[k]^2)\)/4}]}, \ \ 
\ \ \ \ \ \ \ {k, 1, 
            100}];\)\)
samplepointplot  =  
        ListPlot[samplepoints, PlotStyle -> {Yellow, PointSize[0.03]}, 
          DisplayFunction -> Identity];
\(\(Show[samplepointplot, chart, PlotRange -> All, 
        AxesLabel -> {"trace", "det"}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> $DisplayFunction];\)\)

The two  p + I q and p - I q with p > 0 and q ≠ 0 region signal the case that
          trace of A > 0,
          determinant of A > 0, and
         (((trace of A))^2) < 4 determinant of A.
Explanation:
The eigenvalues of A are
           p + I q = 
        (1/2) (trace + (Sqrt[(trace^2) - 4 det])) 
and
           p - I q = (1/2) (trace - (Sqrt[(trace^2) - 4 det])).

Because 
         (trace^2) - 4 det < 0, 
you see q is not 0.  
Because trace > 0, you see that p > 0.

Both real and both negative explanation:

Here are some sample points {trace, det} that plot out in the both real and both negative  region of the chart:

Clear[trace, k];
\(	\(trace[k_]\  := \ \(trace[k] = \ 
          Random[Real, {\(-5\), \(-0.3\)}]\);\)\)
samplepoints  =  
        Table[{trace[k],  Random[Real, {0, (trace[k]^2)/4}]}, {k, 1, 
            100}];
samplepointplot  =  
        ListPlot[samplepoints, PlotStyle -> {ManganeseBlue, PointSize[0.03]}, 
          DisplayFunction -> Identity];
\(\(Show[samplepointplot, chart, 
        PlotRange -> 
          All, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> $DisplayFunction];\)\)

The both real and both negative region corresponds to the case that
          trace of A < 0,
          determinant of A > 0, and
         (((trace of A))^2) > 4 determinant of A.
Explanation:
  
The eigenvalues of A are
         (1/2) (trace - (Sqrt[(trace^2) - 4 det])) 
and
         (1/2) (trace + (Sqrt[(trace^2) - 4 det])).

Because 
        (trace^2) - 4 det > 0, 
you see that neither eigenvalue involves I = (Sqrt[-1]). So there is no whirl.
Because det > 0, you are guaranteed that
       \(\(\(\ \)\(0 < \@\(trace\^2 - 4\ det\) < \@trace\^2 = \(-trace\)\)\)\).

This guarantees that both eigenvalues
       (1/2) (trace - (Sqrt[(trace^2) - 4 det])) < 0
       (1/2) (trace + (Sqrt[(trace^2) - 4 det])) < 0.
The upshot: Both eigenvalues are negative.

Both real and both positive explanation:

Here are some sample points {trace, det} that plot out in the two pure propeller region of the chart:

Clear[trace, k];
\(	\(trace[k_]\  := \ \(trace[k] = \ 
          Random[Real, {0.3, 5}]\);\)\)
samplepoints  =  
        Table[{trace[k],  Random[Real, {0, (trace[k]^2)/4}]}, {k, 1, 
            100}];
samplepointplot  =  
        ListPlot[samplepoints, PlotStyle -> {LightCoral, PointSize[0.03]}, 
          DisplayFunction -> Identity];
\(\(Show[samplepointplot, chart, 
        PlotRange -> 
          All, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> $DisplayFunction];\)\)

The both real and both positive  region corresponds to case that
          trace of A > 0,
          determinant of A > 0, and
         (((trace of A))^2) > 4 determinant of A.
Explanation:
  
The eigenvalues of A are
         (1/2) (trace - (Sqrt[(trace^2) - 4 det])) 
and
         (1/2) (trace + (Sqrt[(trace^2) - 4 det]))
        

Because 
        (trace^2) - 4 det > 0, 
you see that neither eigenvalue involves I = (Sqrt[-1]). So there is no whirl.
Because det > 0, you are guaranteed that
       0 < (Sqrt[(trace^2) - 4 det]) < (Sqrt[(trace^2)]) = trace.

This guarantees that both eigenvalues
       (1/2) (trace - (Sqrt[(trace^2) - 4 det])) > 0
       (1/2) (trace + (Sqrt[(trace^2) - 4 det])) > 0.
The upshot: Both eigenvalues are positive.

Both real with one positive and one negative explanation

Here are some sample points {trace, det} that plot out in the both real with one positive and one negative region of the chart:

Clear[trace, k];
\(	\(trace[k_]\  := \ \(trace[k] = \ 
          Random[Real, {\(-5\), 5}]\);\)\)
samplepoints  =  
        Table[{trace[k],  Random[Real, {-4, 0}]}, {k, 1, 
            100}];
samplepointplot  =  
        ListPlot[samplepoints, PlotStyle -> {Turquoise, PointSize[0.03]}, 
          DisplayFunction -> Identity];
\(\(Show[samplepointplot, chart, 
        PlotRange -> 
          All, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ DisplayFunction -> $DisplayFunction];\)\)

The both real with one positive and one negative  region corresponds to the case that
          determinant of A < 0.
Explanation:
  
The eigenvalues of A are
         (1/2) (trace - (Sqrt[(trace^2) - 4 det])) 
and
         (1/2) (trace + (Sqrt[(trace^2) - 4 det])).

Because det < 0, you are guaranteed that
         (trace^2) - 4 det > 0,
so that neither eigenvalue involves I = (Sqrt[-1]). So there is no whirl.
Also because det < 0, you are guaranteed that
        (Sqrt[(trace^2) - 4 det]) > (Sqrt[(trace^2)]) = Abs[trace].
This tells you that A has one negative eigenvalue - namely
        (1/2) (trace - (Sqrt[(trace^2) - 4 det])).
And the other eigenvalue
        (1/2) (trace + (Sqrt[(trace^2) - 4 det]))
 is automatically positive.