MGM.09 Eigensense: Diagonalizable Matrices,
Matrix Exponential, Matrix Powers and Dynamical Systems
GIVE IT A TRY!

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].";

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

G.1) Eigensense*

A kD vector X is an eigenvector of a  kD matrix A if
               A . X 
points in the same or the opposite direction as X.
In this case there is a number s (called corresponding eigenvalue)  so that 
            A . X = s X.

G.1.a.i) Eigenvectors with positive eigenvalues

Here's a matrix A:

A=(({{0.9, -0.3}, {-0.4, 1.3}}));
MatrixForm[A]

Here are its eigenvectors shown with a random unit vector:

Clear[eigenvector];
\(\({eigenvector[1], eigenvector[2]} = Eigenvectors[A];\)
\)
s = Random[Real, {0, N[2 π]}];
randomvector = {Cos[s], Sin[s]};
\(\(ranger = 1.5;\)
\)
before = 
        Show[Arrow[eigenvector[1], Tail -> {0, 0}, 
            VectorColor -> Indigo], 
          Arrow[eigenvector[2], Tail -> {0, 0}, 
            VectorColor -> Indigo], 
          Arrow[randomvector, Tail -> {0, 0}, 
            VectorColor -> Indigo], Axes -> True, 
          PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
          AxesLabel -> {"x", "y"}];

Here's what you get when you hit all three plotted vectors by the matrix A:

after = 
        Show[Arrow[A . eigenvector[1], Tail -> {0, 0}, 
            VectorColor -> Red], 
          Arrow[A . eigenvector[2], Tail -> {0, 0}, 
            VectorColor -> Red], 
          Arrow[A . randomvector, Tail -> {0, 0}, 
            VectorColor -> Red], 
          PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
          Axes -> True, AxesLabel -> {"x", "y"}, 
          PlotLabel -> "After"];

Here are both plots:

Show[before, after, PlotLabel -> "Both"];

Describe how this plot reveals which two of the original vectors are the eigenvectors of A.
Describe how this plot reveals that both eigenvalues of A are positive with one eigenvalue smaller than 1 and the other eigenvalue bigger than 1.

G.1.a.ii) One eigenvector with a positive eigenvalue; one eigenvector with a negative eigenvalue

Here's a new matrix A:

A=(({{1.1, -0.7}, {0.7, -0.8}}));
MatrixForm[A]

Here are its eigenvectors shown with a random unit vector:

Clear[eigenvector];
\(\({eigenvector[1], eigenvector[2]} = Eigenvectors[A];\)
\)
s = Random[Real, {0, N[2 π]}];
randomvector = {Cos[s], Sin[s]};
\(\(ranger = 1.7;\)
\)
before = 
        Show[Arrow[eigenvector[1], Tail -> {0, 0}, 
            VectorColor -> Indigo], 
          Arrow[eigenvector[2], Tail -> {0, 0}, 
            VectorColor -> Indigo], 
          Arrow[randomvector, Tail -> {0, 0}, 
            VectorColor -> Indigo], Axes -> True, 
          PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
          AxesLabel -> {"x", "y"}, 
          PlotLabel -> "Before"];

Here's what you get when you multiply all three plotted vectors by the matrix A:

after = 
        Show[Arrow[A . eigenvector[1], Tail -> {0, 0}, 
            VectorColor -> Red], 
          Arrow[A . eigenvector[2], Tail -> {0, 0}, 
            VectorColor -> Red], 
          Arrow[A . randomvector, Tail -> {0, 0}, 
            VectorColor -> Red], 
          PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
          Axes -> True, AxesLabel -> {"x", "y"}, 
          PlotLabel -> "After"];

Here are both plots:

Show[before, after, PlotLabel -> "Both"];

Describe how this plot reveals which two of the original vectors are the eigenvectors of A.
Describe how this plot reveals that one of the eigenvalues is positive and that the other is negative.

G.1.b.i) Pick out the eigenvectors

Here's a new matrix A:

A=(({{1.425, -0.129904}, {-0.129904, 1.275}}));
MatrixForm[A]

Here are a bunch of color-coded unit vectors:

Clear[vectorcolor, t];
\(\(vectorcolor[t_] = 
        RGBColor[0.5\ \((Cos[t] + 1)\), 0.5\ \((Sin[t] + 1)\), 
          0.5];\)
\)
\(\(unitvector[t_] = {Cos[t], Sin[t]};\)
\)
before = 
        Show[Table[
            Arrow[unitvector[t], Tail -> {0, 0}, 
              VectorColor -> vectorcolor[t]], {t, 0, 
              2 π, (π/6)}], Axes -> True, 
          AxesLabel -> {"x", "y"}, 
          PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}}, 
          PlotLabel -> "Before"];

Here are the same vectors after they have been multiplied by A.

after = 
        Show[Table[
            Arrow[A . unitvector[t], Tail -> {0, 0}, 
              VectorColor -> vectorcolor[t]], {t, 0, 
              2 π, (π/6)}], Axes -> True, 
          AxesLabel -> {"x", "y"}, 
          PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}}, 
          PlotLabel -> "After"];

Here are both plots together:

Show[before, after, PlotLabel -> "Both"];

In this plot, you see enough information to pick out the two eigenvectors of A and their negatives. Say how you made the identification.

G.1.b.ii) Rotation matrices can't have real eigenvectors

Here is a random number s  with 0 < s < (π/2) together with the matrix A whose hits rotate everything about {0,0} by s counterclockwise radians :

s  =  Random[Real, {0.1, (π/2)  -  0.1}];
A  =  
        Transpose[{{Cos[s], Sin[s]}, {Cos[s + (π/2)], 
              Sin[s + (π/2)]}}];
MatrixForm[A]

Here's what happens when you diagonalize A:

Clear[eigenvector];
spanners  =   Eigenvectors[A];
SpannerMatrix  =  Transpose[spanners];
MatrixForm[
      SpannerMatrix]

The diagonal matrix comes from the eigenvalues of A:

diagonalmatrix  =  DiagonalMatrix[Eigenvalues[A]];
MatrixForm[
      diagonalmatrix]

Here's a look at
          SpannerMatrix.diagonalmatrix.(SpannerMatrix^(-1)):

MatrixForm[SpannerMatrix.diagonalmatrix.Inverse[SpannerMatrix]]

Here's a look at A:

MatrixForm[A]

So:
         A  =  SpannerMatrix . diagonalmatrix . (SpannerMatrix^(-1))
and A is diagonalizable.
Notice that complex numbers were needed to do the diagonalization. 
Explain why you could have been predicted this in advance.

G.2) Diagonalizations, matrix exponentials and matrix powers

G.2.a.i) Diagonalizing  a 2D matrix

Here's a 2D matrix A:

A = (({{2.1972, -0.9859}, {0.7394, -1.7804}}));
MatrixForm[A]

Here are the eigenvalues of A:

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

And the corresponding unit eigenvectors of A:

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

See them:

\(\(eigenplot\  = \ 
\ \ \ \ \ \ \ \ \ {Arrow[eigenvector[1], 
            Tail -> {0, 0}, VectorColor -> Red], 
\ \ \ \ \ \ \ \ \ \ Arrow[
            eigenvector[2], Tail -> {0, 0}, 
            VectorColor -> Indigo]};\)
\)
\(\(labels\  = \ {Graphics[{Black, 
              Text["eigenvector[1]", 
                0.5\ \ eigenvector[
                    1]]}], 
\ \ \ \ \ \ \ \ \ \ Graphics[{Black, 
              Text["eigenvector[2]", 
                0.5\ \ eigenvector[2]]}]};\)\)
\(\(ranger\  = \ 1;\)
\)
\(\(Show[eigenplot, labels, Axes -> True, 
        AxesLabel -> {"x", "y"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotRange -> {{\(-ranger\), ranger}, {\(-ranger\), 
              ranger}}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotLabel -> "Two eigenvectors"];\)\)

The eigenvectors span all of 2D and this tells you that A is diagonalizable.
Diagonalize A by coming up with SpannerMatrix and diagonalmatrix 
so that
          SpannerMatrix.diagonalmatrix.(SpannerMatrix^(-1))
duplicates A.

G.2.a.ii) Building the matrix exponential  (E^(A t))

Stay with the same matrix A as in part i):

A = (({{2.1972, -0.9859}, {0.7394, -1.7804}}));
MatrixForm[A]

In part i), you diagonalized A by coming up with SpannerMatrix and diagonalmatrix 
so that
          A = SpannerMatrix.diagonalmatrix.(SpannerMatrix^(-1)).
          
Stay with the same SpannerMatrix, but come up with newdiagonalmatrix so that
          (E^(A t))= SpannerMatrix.newdiagonalmatrix.(SpannerMatrix^(-1)).
Check yourself by activating:

Clear[t];
MatrixForm[Simplify[MatrixExp[A t]]]

G.2.b.i) Reading the eigenvalues of A from the formula for  (E^(A t))

Here's new 2D matrix A:

A = (({{-1.48, -1.43}, {-2.55, 3.82}}));
MatrixForm[A]

Look at this calculation of   (E^(A t)):

Clear[t];
MatrixForm[Simplify[MatrixExp[A t]]]

If you are clued in, you can read the eigenvalues of A from the output directly above. 
Type them here:
           eigenvalue[1] = .......................... 
           eigenvalue[2] = .......................... .

G.2.b.ii) Another one

Here's new 2D matrix A:

A =(({{-1.48, -1.43}, {2.55, 0.82}}));
MatrixForm[A]

Look at this calculation of   (E^(A t)):

MatrixForm[Simplify[ComplexExpand[MatrixExp[A t]]]]

If you are clued in, you can read the eigenvalues p + I q and p - I q of A from the output directly above. 
Type them here:
           eigenvalue[1] = .......................... 
           eigenvalue[2] = .......................... .

G.2.c) The square root of a matrix

Here's new 2D matrix A:

A =(({{2.48, 1.03}, {1.13, 0.82}}));
MatrixForm[A]

Apply the square root function to this matrix to make a new matrix B with
                   B.B = A

Eigenvalues[A]

G.2.d) The derivative of  (E^(A t)) with respect to t

Here's new 2D matrix A:

A = (({{Random[Real, {-3, 3}], Random[
                    Real, {-3, 3}]}, {Random[Real, {-3, 3}], Random[
                    Real, {-3, 3}]}})) ;
MatrixForm[A]

Look at this calculation of D[(E^(A t)),t]  (the derivative of (E^(A t)) with respect to t):

MatrixForm[Simplify[ComplexExpand[D[MatrixExp[A t], t]]]]

Look at this calculation of  A . (E^(A t)):

MatrixForm[Simplify[ComplexExpand[A . MatrixExp[A t]]]]

Rerun all three cells several times and then say why you are not surprised with what you see.

G.2.e.i) Matrix powers (A^k) for rotation matrices A

Here is a random right hand perpendicular frame:

s = Random[Real, {0.3, 0.7}];
\(\(ranger = 1;\)
\)
Clear[perpframe];
\(\({perpframe[1], 
          perpframe[2]} = {{Cos[s], Sin[s]}, {Cos[s\  + \ π\/2], 
            Sin[s\  + \ π\/2]}};\)
\)
\(\(frameplot = {Table[
            Arrow[perpframe[k], Tail -> {0, 0}, 
              VectorColor -> Indigo, HeadSize -> 0.2], {k, 1, 2}], 
          Graphics[Text["perpframe[1]", 0.6\ perpframe[1]]], 
          Graphics[
            Text["perpframe[2]", 0.6\ perpframe[2]]]};\)
\)
before = 
        Show[frameplot, 
          PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
          Axes -> True, AxesLabel -> {"x", "y"}, 
          PlotLabel -> s "=s"];

When you want to rotate everything by s counterclockwise radians you hit with this matrix:

A  =  Transpose[{perpframe[1], perpframe[2]}];
MatrixForm[A]

See what happens when you take two random starting points and hit them with increasing matrix powers of A:

\(
\(Clear[iterationplot, interationplotter, scaler, point, pointcolor, 
        starter, khigh, k, j, iterationpoints];\)\)
point[k_, starter_]  := MatrixPower[A, k]  . starter ;
\(pointcolor[k_, khigh_]\  := \ 
      RGBColor[Sin[\((Pi/2)\)\ k/\((khigh + 0.1)\)], 
        Cos[\((Pi/2)\)\ k/\((khigh\  + \ 0.1)\)], 0]\ 
\)
\(\(iterationpoints[starter_, khigh_]\  := 

        Table[Graphics[\ \ {pointcolor[k, khigh], PointSize[0.03], 
              Point[MatrixPower[A, k]\  . starter]}], {k, 0, 
            khigh}];\)
\)
\(\(iterationplotter[khigh_] := 
        Show[Table[iterationpoints[starter[j], khigh], {j, 1, 2}], 
          Axes -> True, AxesLabel -> {"x", "y"}, 
          PlotLabel -> \*""Hits with A,(A^2),..,(A^khigh)"", 
		PlotRange -> {{\(-3\), 3}, {\(-3\), 3}}, AspectRatio -> Automatic, 
          DisplayFunction -> Identity];\)

\)
ss  =  Random[Real, {0,  2 Pi}];
sss  =  Random[Real, {0,  2 Pi}];
starter[1] =  
        Random[Real, {0.5, 1}] {Cos[ss], Sin[ss]} ;
\(\(starter[
          2] = \ \ Random[Real, {1.5, 3}] {Cos[sss], 
            Sin[sss]}\ ;\)

\)
khigh  =  0;
\(\(Show[iterationplotter[khigh], PlotLabel -> "Two starting points", 
        DisplayFunction -> $DisplayFunction];\)
\)

Here are the same two starter points together with the points you get when you hit each of the original starting points with A:

khigh  =  1;
Show[iterationplotter[khigh], 
        DisplayFunction -> $DisplayFunction];

Here are the same two starter points together with the points you get when you hit each of the original starting points with A and A.A:

khigh  =  2;
Show[iterationplotter[khigh], 
        DisplayFunction -> $DisplayFunction];

Here are the same two starter points together with the points you get when you hit each of the original starting points with A, A . A = (A^2), A . A . A = (A^3), A . A . A . A = (A^4), 
      and A . A . A . A . A = (A^5):

khigh  =  5;
Show[iterationplotter[khigh], 
        DisplayFunction -> $DisplayFunction];

Here are the same two starter points together with the points you get when you hit each of the original starting points with A, A . A = (A^2), A . A . A = (A^3),  .   .   . , (A^9), 
      and (A^10):

khigh  =  10;
Show[iterationplotter[khigh], 
        DisplayFunction -> $DisplayFunction];

Here are the same four starter points together with the points you get when you hit each of the original starting points with A, A . A = (A^2), A . A . A = (A^3),  .   .   . , (A^49), 
      and (A^50):

khigh  =  50;
Show[iterationplotter[khigh], 
        DisplayFunction -> $DisplayFunction];

Describe what you see and explain why you see it.
Is  every integral power of A  (such as (A^4)) also a rotation matrix?

G.2.e.ii) If A is rotation matrix, does  (E^(A t))oscillate on closed curves as t advances

Here's a new right hand perpendicular frame and the rotation matrix based on it:

s = Random[Real, {0.3, 0.7}];
\(\(ranger = 1;\)
\)
Clear[perpframe];
\(\({perpframe[1], 
          perpframe[2]} = {{Cos[s], Sin[s]}, {Cos[s\  + \ π\/2], 
            Sin[s\  + \ π\/2]}};\)
\)
\(\(frameplot = {Table[
            Arrow[perpframe[k], Tail -> {0, 0}, 
              VectorColor -> Indigo, HeadSize -> 0.2], {k, 1, 2}], 
          Graphics[Text["perpframe[1]", 0.6\ perpframe[1]]], 
          Graphics[
            Text["perpframe[2]", 0.6\ perpframe[2]]]};\)
\)
before = 
        Show[frameplot, 
          PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
          Axes -> True, AxesLabel -> {"x", "y"}, 
          PlotLabel -> s "=s"];
rotationmatrix  =  
        Transpose[{perpframe[1], perpframe[2]}];
MatrixForm[rotationmatrix]

When you want to use this perpendicular frame to rotate, you hit with this matrix:

A  =  rotationmatrix;
MatrixForm[A]

Look at the eigenvalues of A:

Eigenvalues[A]

Use what you see to answer this question:
If A is a rotation matrix, does (E^(A t))oscillate on closed curves as t advances?

G.2.e.iii) For which rotation matrices A does  (E^(A t))oscillate on closed curves as t advances

Here is a cleared rotation matrix:

Clear[perpframe, s];
\(\({perpframe[1], 
          perpframe[2]}\  = \ 
{{Cos[s], Sin[s]}, {Cos[s + π\/2], 
            Sin[s + π\/2]}};\)\)
\(\(A\  = \ Transpose[{perpframe[1], perpframe[2]}];\)
\)
MatrixForm[A]

Here is Mathematica's calculation of the eigenvalues of A:

Eigenvalues[A]

The question is this:
What values of s (with 0 ≤ s < 2π)  guarantee that (E^(A t))oscillates on closed curves as t advances?
Show off your answer with some plots.

Fill in your pick for s and run this code.

s = ??????;
				
Clear[trajectoryplot, eigenplot, scaler, eigenvector, 
              starter, thigh, k, j];

trajectoryplot[thigh_, starter_]  :=  
            ParametricPlot[MatrixExp[A t] . starter, {t, 0, thigh}, 
              PlotStyle -> {{Thickness[0.01], 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[{NavyBlue, PointSize[0.03], 
                  Point[starter[j]]}], {j, 1, 6}]; 

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}]};
starter[5] =  {Random[Real, {-10, -4}], 0};
starter[6] =  {Random[Real, {4, 10}], 0};

thigh  =  0.01;
Show[Table[trajectoryplot[thigh, starter[j]], {j, 1, 6}] , 
            starterplots, PlotLabel -> "Starting points", 
            DisplayFunction -> $DisplayFunction];

Here is how the curves
 \(E\^\(A\ t\) . starter[1], E\^\(A\ t\) . starter[2], 
      E\^\(A\ t\) . starter[3], 
E\^\(A\ t\) . starter[4], 
      E\^\(A\ t\) . starter[5], E\^\(A\ t\) . starter[6]\)
plot out as t advances from 0 to 8

\(\(thigh\  = \ 2;\)
\)
\(\(Show[Table[trajectoryplot[thigh, starter[k]], {k, 1, 6}]\ , 
	
        Table[pointer[thigh, starter[k]], {k, 1, 6}], starterplots, 
        PlotRange -> All, 
        PlotLabel -> \*""Plots of \!\(E\^\(A\\\\\\\\ t\)\).starter[k] for 0 < t"" < thigh, 
        DisplayFunction -> $DisplayFunction];\)\)

G.3) Using eigenvalues to analyze behavior of solutions of systems of linear differential equations

G.3.a.i) A system of linear diffeq's

To get a 2D system of linear differential equations, you go with a 2D matrix coefficient matrix A and put
          {x'[t],y'[t]}=A . {x[t], y[t]} . 
Here's one:

A =(({{-1.4, 0.8}, {0.5, -0.7}}));
MatrixForm[
      A]
Clear[x, y, 
        t];
linearsystem  =  ({x'[t], y'[t]}  ==  
            A . {x[t], y[t]});
ColumnForm[
      Thread[linearsystem]]

Look at this plot of six trajectories coming from solutions 
             {x[t], y[t]} = (E^(A t)) . starter 
 of this system of differential equations starting at six random points as t advances from 0 to 1.5:

Clear[trajectoryplot, starterplot, starter, thigh, 
        matrixexpA];
matrixexpA[t_]  =  MatrixExp[A t];
\(	\(trajectoryplot[starter_, 
          thigh_]\  := 
	\ {ParametricPlot[\ \ matrixexpA[t]\  . 
              starter, {t, 0, thigh}, 
            PlotStyle -> {{Thickness[0.01], 
                  CadmiumOrange}}, \ \ \ \ 
			\ \ \ PlotLabel -> "Solutions of {x'[t],y'[t]} = A.{x[t],y[t]}", \ \ \ 
			\ \ \ PlotRange -> All, \ AxesLabel -> {"x", "y"}, AxesOrigin -> {0, 0}, 
            DisplayFunction -> Identity], 
		
		ArrowHead[
            matrixexpA[0.85\ \ thigh]\  . 
              starter, \(matrixexpA'\)[0.85\ thigh]\  . 
              starter, 
			\ \ \ \ \ \ \ \ HeadSize -> 0.2, 
            VectorColor -> Black, 
            Aperture -> 0.4], 
		
\ \ \ \ \ Graphics[{NavyBlue, 
              PointSize[0.03], Point[starter]}], 
          Graphics[{NavyBlue, 
              Text[thigh\ "= thigh", {\(-0.7\), 1}]}]};\)
\)
starter[1] =  {Random[Real, {1, 3}], 
          Random[Real, {1, 3}]};
starter[2] =  {Random[Real, {-3, -1}], 
          Random[Real, {-3, -1}]};
starter[3] =  {Random[Real, {-3, -1}], 
          Random[Real, {1, 3}]};
starter[4] =  {Random[Real, {1, 3}], 
          Random[Real, {-3, -1}]};
starter[5] =  {Random[Real, {-3, 3}], 
          Random[Real, {-3, 3}]};
\(\(starter[6] = \ {Random[Real, {\(-3\), 3}], 
          Random[Real, {\(-3\), 3}]};\)
\)
thigh  =  1.5;
\(\(Show[Table[trajectoryplot[starter[k], thigh], {k, 1, 6}], 

        DisplayFunction -> $DisplayFunction];\)\)

Here's what happens as t advances from 0 to 7.5:

thigh  =  7.5;
\(\(Show[Table[trajectoryplot[starter[k], thigh], {k, 1, 6}], 

        DisplayFunction -> $DisplayFunction];\)\)

Here are the eigenvalues of the coefficient matrix A:

Eigenvalues[A]

Use the eigenvalue information to explain why the plots come out the way they do. 
How will trajectories coming from other solutions look?
Are they all attracted to {0,0} as t gets large?

G.3.a.ii) Swirling around {0,0}

To get a 2D system of linear differential equations, you go with a 2D matrix coefficient matrix A and put
          {x'[t],y'[t]}=A . {x[t], y[t]} . 
Here's another one:

A = (({{-1.4, 1.6}, {-1.5, -0.7}}));
MatrixForm[
      A]
Clear[x, y, 
        t];
linearsystem  =  ({x'[t], y'[t]}  ==  
            A . {x[t], y[t]});
ColumnForm[
      Thread[linearsystem]]

Look at this plot of six trajectories coming from solutions  
         {x[t], y[t]} = (E^(A t)) . starter 
starting at six random points as t advances from 0 to 1.5:

Clear[trajectoryplot, starterplot, starter, thigh, 
        matrixexpA];
matrixexpA[t_]  =  MatrixExp[A t];
\(	\(trajectoryplot[starter_, 
          thigh_]\  := 
	\ {ParametricPlot[\ \ matrixexpA[t]\  . 
              starter, {t, 0, thigh}, 
            PlotStyle -> {{Thickness[0.01], 
                  CadmiumOrange}}, \ \ \ \ 
			\ \ \ PlotLabel -> "Solutions of {x'[t],y'[t]} = A.{x[t],y[t]}", \ \ \ 
			\ \ \ PlotRange -> All, \ AxesLabel -> {"x", "y"}, AxesOrigin -> {0, 0}, 
            DisplayFunction -> Identity], 
		
		ArrowHead[
            matrixexpA[0.85\ \ thigh]\  . 
              starter, \(matrixexpA'\)[0.85\ thigh]\  . 
              starter, 
			\ \ \ \ \ \ \ \ HeadSize -> 0.2, 
            VectorColor -> Black, 
            Aperture -> 0.4], 
		
\ \ \ \ \ Graphics[{NavyBlue, 
              PointSize[0.03], Point[starter]}], 
          Graphics[{NavyBlue, 
              Text[thigh\ "= thigh", {\(-1\), 1}]}]};\)
\)
starter[1] =  {Random[Real, {1, 3}], 
          Random[Real, {1, 3}]};
starter[2] =  {Random[Real, {-3, -1}], 
          Random[Real, {-3, -1}]};
starter[3] =  {Random[Real, {-3, -1}], 
          Random[Real, {1, 3}]};
starter[4] =  {Random[Real, {1, 3}], 
          Random[Real, {-3, -1}]};
starter[5] =  {Random[Real, {-3, 3}], 
          Random[Real, {-3, 3}]};
\(\(starter[6] = \ {Random[Real, {\(-3\), 3}], 
          Random[Real, {\(-3\), 3}]};\)
\)
thigh  =  1.5;
\(\(Show[Table[trajectoryplot[starter[k], thigh], {k, 1, 6}], 

        DisplayFunction -> $DisplayFunction];\)\)

Here's what happens as t advances from 0 to 4.5:

thigh  =  4.5;
\(\(Show[Table[trajectoryplot[starter[k], thigh], {k, 1, 6}], 

        DisplayFunction -> $DisplayFunction];\)\)

Here are the eigenvalues of the coefficient matrix A:

Eigenvalues[A]

Use the eigenvalue information to explain why the plots come out the way they do. 
How will trajectories coming from other solutions look?
Are they all attracted to {0,0} as t gets large?

G.3.a.iii) Moving away from {0,0}

To get a 2D system of linear differential equations, you go with a 2D matrix coefficient matrix A and put
          {x'[t],y'[t]}=A . {x[t], y[t]} . 
Here's another yet another one:

A = (({{1.4, -0.4}, {-0.5, 0.6}}));
MatrixForm[
      A]
Clear[x, y, 
        t];
linearsystem  =  ({x'[t], y'[t]}  ==  
            A . {x[t], y[t]});
ColumnForm[
      Thread[linearsystem]]

Look at this plot of six trajectories coming from solutions  
         {x[t], y[t]} = (E^(A t)) . starter 
starting at six random points as t advances from 0 to 1.5:

Clear[trajectoryplot, starterplot, starter, thigh, 
        matrixexpA];
matrixexpA[t_]  =  MatrixExp[A t];
\(	\(trajectoryplot[starter_, 
          thigh_]\  := 
	\ {ParametricPlot[\ \ matrixexpA[t]\  . 
              starter, {t, 0, thigh}, 
            PlotStyle -> {{Thickness[0.01], 
                  CadmiumOrange}}, \ \ \ \ 
			\ \ \ PlotLabel -> "Solutions of {x'[t],y'[t]} = A.{x[t],y[t]}", \ \ \ 
			\ \ \ PlotRange -> All, \ AxesLabel -> {"x", "y"}, AxesOrigin -> {0, 0}, 
            DisplayFunction -> Identity], 
		
		ArrowHead[
            matrixexpA[0.85\ \ thigh]\  . 
              starter, \(matrixexpA'\)[0.85\ thigh]\  . 
              starter, 
			\ \ \ \ \ \ \ \ HeadSize -> 0.2, 
            VectorColor -> Black, 
            Aperture -> 0.4], 
		
\ \ \ \ \ Graphics[{NavyBlue, 
              PointSize[0.03], Point[starter]}], 
          Graphics[{NavyBlue, 
              Text[thigh\ "= thigh", {\(-1\), 1}]}]};\)
\)
starter[1] =  {Random[Real, {1, 3}], 
          Random[Real, {1, 3}]};
starter[2] =  {Random[Real, {-3, -1}], 
          Random[Real, {-3, -1}]};
starter[3] =  {Random[Real, {-3, -1}], 
          Random[Real, {1, 3}]};
starter[4] =  {Random[Real, {1, 3}], 
          Random[Real, {-3, -1}]};
starter[5] =  {Random[Real, {-3, 3}], 
          Random[Real, {-3, 3}]};
\(\(starter[6] = \ {Random[Real, {\(-3\), 3}], 
          Random[Real, {\(-3\), 3}]};\)
\)
thigh  =  1.5;
\(\(Show[Table[trajectoryplot[starter[k], thigh], {k, 1, 6}], 

        DisplayFunction -> $DisplayFunction];\)\)

Here's what happens as t advances from 0 to 7.5:

thigh  =  7.5;
\(\(Show[Table[trajectoryplot[starter[k], thigh], {k, 1, 6}], 

        DisplayFunction -> $DisplayFunction];\)\)

Here are the eigenvalues of the coefficient matrix A:

Eigenvalues[A]

Use this information to help to identify a unit vector pointing in the same direction as the line you see above.

G.3.a.iv) Swirling away from {0,0}

To get a 2D system of linear differential equations, you go with a 2D matrix coefficient matrix A and put
          {x'[t],y'[t]}=A . {x[t], y[t]} . 
Here's another one:

A = (({{0.4, -0.8}, {0.5, 0.3}}));
MatrixForm[
      A]
Clear[x, y, 
        t];
linearsystem  =  ({x'[t], y'[t]}  ==  
            A . {x[t], y[t]});
ColumnForm[
      Thread[linearsystem]]

Look at this plot of six trajectories coming from solutions  
         {x[t], y[t]} = (E^(A t)) . starter 
starting at six random points as t advances from 0 to 1.5:

Clear[trajectoryplot, starterplot, starter, thigh, 
        matrixexpA];
matrixexpA[t_]  =  MatrixExp[A t];
\(	\(trajectoryplot[starter_, 
          thigh_]\  := 
	\ {ParametricPlot[\ \ matrixexpA[t]\  . 
              starter, {t, 0, thigh}, 
            PlotStyle -> {{Thickness[0.01], 
                  CadmiumOrange}}, \ \ \ \ 
			\ \ \ PlotLabel -> "Solutions of {x'[t],y'[t]} = A.{x[t],y[t]}", \ \ \ 
			\ \ \ PlotRange -> All, \ AxesLabel -> {"x", "y"}, AxesOrigin -> {0, 0}, 
            DisplayFunction -> Identity], 
		
		ArrowHead[
            matrixexpA[0.95\ \ thigh]\  . 
              starter, \(matrixexpA'\)[0.95\ thigh]\  . 
              starter, 
			\ \ \ \ \ \ \ \ HeadSize -> 0.6, 
            VectorColor -> Black, 
            Aperture -> 0.4], 
		
\ \ \ \ \ Graphics[{Blue, 
              PointSize[0.03], Point[starter]}], 
          Graphics[{Blue, 
              Text[thigh\ "= thigh", {\(-1\), 1}]}]};\)
\)
starter[1] =  {Random[Real, {1, 3}], 
          Random[Real, {1, 3}]};
starter[2] =  {Random[Real, {-3, -1}], 
          Random[Real, {-3, -1}]};
starter[3] =  {Random[Real, {-3, -1}], 
          Random[Real, {1, 3}]};
starter[4] =  {Random[Real, {1, 3}], 
          Random[Real, {-3, -1}]};
starter[5] =  {Random[Real, {-3, 3}], 
          Random[Real, {-3, 3}]};
\(\(starter[6] = \ {Random[Real, {\(-3\), 3}], 
          Random[Real, {\(-3\), 3}]};\)
\)
thigh  =  1.5;
\(\(Show[Table[trajectoryplot[starter[k], thigh], {k, 1, 6}], 

        DisplayFunction -> $DisplayFunction];\)\)

Here's what happens as t advances from 0 to 4.5:

thigh  =  4.5;
\(\(Show[Table[trajectoryplot[starter[k], thigh], {k, 1, 6}], 

        DisplayFunction -> $DisplayFunction];\)\)

Here are the eigenvalues of the coefficient matrix A:

Eigenvalues[A]

Use the eigenvalue information to explain why the plots come out the way they do. 
How will trajectories coming from other solutions look?

G.3.a.v) Merging onto a line

To get a 2D system of linear differential equations, you go with a 2D matrix coefficient matrix A and put
          {x'[t],y'[t]}=A . {x[t], y[t]} . 
Here's another one:

A = (({{1.4, 0.8}, {0.5, -0.9}}));
MatrixForm[
      A]
Clear[x, y, 
        t];
linearsystem  =  ({x'[t], y'[t]}  ==  
            A . {x[t], y[t]});
ColumnForm[
      Thread[linearsystem]]

Look at this plot of six trajectories coming from solutions  
         {x[t], y[t]} = (E^(A t)) . starter 
starting at six random points as t advances from 0 to 1.5:

Clear[trajectoryplot, starterplot, starter, thigh, 
        matrixexpA];
matrixexpA[t_]  =  MatrixExp[A t];
\(	\(trajectoryplot[starter_, 
          thigh_]\  := 
	\ {ParametricPlot[\ \ matrixexpA[t]\  . 
              starter, {t, 0, thigh}, 
            PlotStyle -> {{Thickness[0.01], 
                  CadmiumOrange}}, \ \ \ \ 
			\ \ \ PlotLabel -> "Solutions of {x'[t],y'[t]} = A.{x[t],y[t]}", \ \ \ 
			\ \ \ PlotRange -> All, \ AxesLabel -> {"x", "y"}, AxesOrigin -> {0, 0}, 
            DisplayFunction -> Identity], 
		
		ArrowHead[
            matrixexpA[0.85\ \ thigh]\  . 
              starter, \(matrixexpA'\)[0.85\ thigh]\  . 
              starter, 
			\ \ \ \ \ \ \ \ HeadSize -> 0.2, 
            VectorColor -> Black, 
            Aperture -> 0.4], 
		
\ \ \ \ \ Graphics[{NavyBlue, 
              PointSize[0.03], Point[starter]}], 
          Graphics[{NavyBlue, 
              Text[thigh\ "= thigh", {\(-1\), 1}]}]};\)
\)
starter[1] =  
        0.5 {Random[Real, {1, 3}], Random[Real, {1, 3}]};
starter[2] =  
        0.5 {Random[Real, {-3, -1}], 
            Random[Real, {-3, -1}]};
starter[3] =  
        0.5 {Random[Real, {-3, -1}], 
            Random[Real, {1, 3}]};
starter[4] =  
        0.5 {Random[Real, {1, 3}], 
            Random[Real, {-3, -1}]};
starter[5] =  
        0.5 {Random[Real, {-3, 3}], 
            Random[Real, {-3, 3}]};
\(\(starter[6] = \ 
        0.5\ {Random[Real, {\(-3\), 3}], 
            Random[Real, {\(-3\), 3}]};\)
\)
thigh  =  0.5;
\(\(Show[Table[trajectoryplot[starter[k], thigh], {k, 1, 6}], 
        AspectRatio -> 1, 

        DisplayFunction -> $DisplayFunction];\)\)

Here's what happens as t advances from 0 to 1.5:

thigh  =  1.5;
\(\(Show[
        Table[trajectoryplot[starter[k], thigh], {k, 1, 
            6}], 
AspectRatio -> 1, 

        DisplayFunction -> $DisplayFunction];\)\)

Here's what happens as t advances from 0 to 10:

thigh  = 10;
\(\(Show[
        Table[trajectoryplot[starter[k], thigh], {k, 1, 
            6}], 
AspectRatio -> 1, 

        DisplayFunction -> $DisplayFunction];\)\)

Here are the eigenvalues of the coefficient matrix A:

Eigenvalues[A]

Use this information to help to identify a unit vector pointing in the same direction as the line you see above.

G.4) The determinant of any diagonalizable matrix is the product of its eigenvalues

G.4.a)  The determinant of any diagonalizable matrix is the product of its eigenvalues

One day, one of your friends comes by, holds out a dusty book and says, "This book I checked out of the library says that the determinant of any diagonalizable matrix is the product of its eigenvalues." 
You decide to look into this and type:

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

Then you calculate the product of the eigenvalues of A:

product  =  Times @@ Eigenvalues[A]

And then you calculate the determinant of A:

determinant  =  Det[A]

Bingo.
You try another:

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

product  =  Times @@ Eigenvalues[A]

determinant  =  Det[A]

You say: "Let me see that book." 
In the book it says, "That the determinant of any diagonalizable matrix is the product of its eigenvalues is a quick consequence of these facts: 
▪ Fact 1: If A,B and C are square matrices each with k rows and k columns, then
                  Det[A.B.C] = Det[A] Det[B] Det[C]
▪ Fact 2: The determinant of any diagonal matrix is the product of its diagonal entries."

You say: "Taking these facts for granted, I see how to use these two facts to explain why the determinant of any diagonalizable matrix is the product of its eigenvalues." 
Your friend says, "Explain it to me."
Do it.

G.5) Continuous dynamical systems: Physical and Chemical Sciences
Analyzing the effect of resistance in simple parallel electrical circuits.
Chemical reactions

G.5.a.i)  Analyzing the effect of resistance in simple parallel electrical circuits

You are fortunate enough to have your own computer set-up right in your dorm room. One night, Brian, the EE student who lives across the hall comes in and says, "Part of my EE 250 homework for tomorrow is to analyze the effect of varying the size of the resistance on the current in a simple parallel electrical circuit.  I don't have a clue about where to start."  

At first, you think, "I'm scared because I'm a life science major and I don't know diddly squat about electrical circuits."  But you don't let on. Instead, you say, "Let me see that EE 250 textbook." It turns out to be thicker and heavier than you had expected.  In the book you see that in a simple parallel electrical circuit, the main measurements are 
       x[t]  =  voltage drop across the capacitor and
       y[t]  =  current through the inductor
Reading on, you spot that these measurements are related through the linear system
       x'[t]  =    y[t]/L
       y'[t]  = -x[t]/c - y[t]/(c r) 
where L,  c,  r are given positive numbers with
        L  =  inductance
        c  =  capacitance, and
        r  =  resistance.
You say, "Brian I don't even know what those words mean, but I think I can help you." 
You say that the coefficient matrix for this linear system is:

\(Clear[A, L, r, c]; 

    A[L_, r_, c_] = {{0, 1\/L}, {\(-\(1\/c\)\), \(-\(1\/\(c\ r\)\)\)}}; 

    MatrixForm[A[L, r, c]]\)

Now you ask Brian, "Is it OK to go with c  =  0.8,  L  =  5.0?
Eagerly, he says, "Sure!"
You enter c  =  0.8,  L  =  5.0 and then look at the eigenvalues of A[L, r, c] for r running from 0.1 to 3.7

c = 0.8; L = 5.0; \(Clear[r]; 


    ColumnForm[
      Table[{r "= r", \((Eigenvalues[A[L, r, c]])\)}, {r, 0.1, 3.7, 0.4}]]
      \)

You say that the current function y[t] will plot out a lot differently for the smaller positive resistance measurements r than for larger resistance measurements.
What is it about the output immediately above that gave you this idea?
What do you mean?
How are the shapes of the plots related to the eigenvalues above?

G.5.a.ii) When the resistance measurement r is very big

Now you ask Brian whether he is interested in what happens when the resistance measurement r is very big?  He replies," I never thought about that.  What do you think?"
You say, "Let's see."

c = 0.8; L = 5.0; \(Clear[r]; 
Table[{r, Eigenvalues[A[L, r, c]]}, {r, 10, 510, 50}]\)

You say that when you go with very large r's, the plots of the corresponding current functions y[t] don't change much as you change r.  
What is it about the output immediately above that gives you this idea?
What do you mean?

G.5.b.i) An eigenvalue = 0

When you come across this linear system
        x'[t]  =  -0.75 x[t]  +  0.25 y[t]
        y'[t]  =    0.75 x[t]  -  0.25 y[t],
you immediately write down this coefficient matrix:

A=(({{-0.75, 0.25}, {0.75, -0.25}}));
MatrixForm[A]

You check it out:

Clear[x, y, t];
linearsystem={x'[t],y'[t]}==A . {x[t], y[t]};

ColumnForm[Thread[linearsystem]]

You ask Mathematica for the eigenvectors and eigenvalues of A:

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

Clear[eigenvalue];
{eigenvalue[1], eigenvalue[2]} = Chop[Eigenvalues[A]]

As explained in the tutorials, solution pairs {x[t], y[t]} are given by
               {x[t], y[t]} = (E^(A t)) . starter.
Explain this: 
If the starter is not a multiple of eigenvector[1]  and  
                {x[t], y[t]} = (E^(A t)) . starter 
is the corresponding solution pair, then as t gets large  y[t]/x[t] tends to:

((-0.948683)/(-0.316228))

G.5.b.ii) Chemical reactions

Given quantities of two reactants X and Y are mixed together. Agree that x[t] measures amount of X at time t units after the reaction begins.  And agree that y[t]measures the amount of Y at time t units after the reaction begins. In the reaction under study here,
        x'[t]  =  - k1 x[t]  +  k2 y[t]
        y'[t]  =     k1 x[t]  -  k2 y[t].  
where k1 and k2 are positive reaction rates.

This system of linear differential equations reflects the fact that X reacts with Y to produce more Y and Y reacts with X to produce more X. Both reactions are reversible. 
Look at the eigenvalues and eigenvectors of the cleared coefficient matrix:

Clear[k1, k2];
A =(({{-k1, k2}, {k1, -k2}}));
MatrixForm[
      A]

Eigenvalues[A]

Eigenvectors[A]

The question here is: 
What is the ultimate ratio Y to X as t goes to infinity.  
Does it depend on the starting values  x[0] and y[0]?

You start with x[0] and y[0] both positive; so it is impossible to start on a multiple of {-1,1}.

G.6) Continuous dynamical systems: Life Sciences
Reservoir Models for drug metabolization
Population models and the ultimate sex ratio

G.6.a.i) Reservoir Models for drug metabolization

Lots of folks deal with the hot issue of models for drug metabolization.  The simplest models for these systems are two compartment models such as this one:

Clear[box, color];
\(box[a_, color_] := 
      Graphics[{color, 
          Line[{{a[[1]] - 1, 
                a[[2]] - 
                  1}, {a[[1]] - 1, 
                a[[2]] + 
                  1}, {a[[1]] + 1, 
                a[[2]] + 
                  1}, {a[[1]] + 1, 
                a[[2]] - 
                  1}, {a[[1]] - 1, 
                a[[2]] - 
                  1}}]}]
\)
\(\(generalmodel = 
        Show[box[{0, 0}, GreenDark], 
          box[{4, 0}, Red], 
		Arrow[{2, 0}, Tail -> {1, 1\/2}], 
          Arrow[{\(-2\), 0}, 
            Tail -> {3, \(-\(1\/2\)\)}], 
		Arrow[{0, \(-1\)}, 
            Tail -> {4, \(-1\)}, 
            VectorColor -> CadmiumYellow], 
		Arrow[{0, \(-1\)}, 
            Tail -> {0, \(-1\)}, 
            VectorColor -> Brown], 
		AspectRatio -> Automatic, 
          Epilog -> {Text["k1", {2, 3\/4}], 
              Text["k2", {2, \(-\(3\/4\)\)}], 
              Text["k3", {15\/4, \(-\(3\/2\)\)}], 
              Text["k4", {1\/4, \(-\(3\/2\)\)}], 
              Text["Gastro-
 Intestinal", {0, 0}], 
              Text["Blood", {4, 0}]}];\)\)

Agree that 
→ x[t] measures the amount of the substance in the gastro-intestinal tract 
→ y[t] measures the amount of the substance in the blood t time units after the substance has been ingested into the gastro-intestinal track.
 
In the case that k3  =  0 and k4  =  0, you get this simplified model:

Show[box[{0, 0}, GreenDark], box[{4, 0}, Red], 
        Arrow[{2, 0}, Tail -> {1, (1/2)}], 
        Arrow[{-2, 0}, Tail -> {3, -(1/2)}], 
        AspectRatio -> Automatic, 
        Epilog -> {Text["k1", {2, (3/4)}], 
            Text["k2", {2, -(3/4)}], 
            Text["Gastro-\n Intestinal", {0, 0}], 
            Text["Blood", {4, 0}]}];

The corresponding linear system is
     x'[t]  =  -k1 x[t]  +  k2 y[t]
     y'[t]  =    k1 x[t]  -  k2 y[t].

This simplified model isn't realistic because it doesn't allow for elimination through urination or defecation. This following diagram sets up the model for elimination through urination.

urinemodel = 
        Show[box[{0, 0}, GreenDark], box[{4, 0}, Red], 
          Arrow[{2, 0}, Tail -> {1, (1/2)}, 
            VectorColor -> NavyBlue], 
          Arrow[{-2, 0}, Tail -> {3, -(1/2)}, 
            VectorColor -> NavyBlue], 
          Arrow[{0, -1}, Tail -> {4, -1}, 
            VectorColor -> CadmiumYellow], 
          AspectRatio -> Automatic, 
          Epilog -> {Text["k1", {2, (3/4)}], 
              Text["k2", {2, -(3/4)}], 
              Text["k3", {(15/4), -(3/2)}], 
              Text["Gastro-\n Intestinal", {0, 0}], 
              Text["Blood", {4, 0}]}];

Write down the linear system resulting from this diagram.

G.6.a.ii) Trying it out

Stay with the model:

Show[urinemodel, AspectRatio -> Automatic];

Go with specific the specific numbers 
         k1  =  0.5,  k2  =  0.5, and k3 = 0.6

Find formulas for x[t] and y[t] given starting data 
        x[0]  =  5 (milligrams) and y[0]  =  0.

Plot {x[t], y[t]} and give the limiting values of x[t] and y[t] as t goes to infinity.

G.6.b.i) Linear population  models

Agree that female[t] measures the population of females and that male[t] measures the population of males t years after a given reference year. One linear model for this setup is
         male'[t] = - dm male[t] + bm (male[t]  +  female[t])
        female'[t] = - df female[t] + bf (male[t]  +  female[t]).

The coefficient matrix for this linear population model is:

Clear[dm, bf, dm, df];
A=(({{bm - dm, bm}, {bf, bf - df}}));
MatrixForm[
      A]

Check:

Clear[male, female, t];
popmodel={male'[t],female'[t]}==A.{male[t], female[t]};

ColumnForm[
      Thread[popmodel]]

This checks.
When you go with bm  =  bf and dm  =  df, so that males and females are equally likely to be born and equally likely to die, then you get

ColumnForm[Thread[popmodel /. {bm -> bf, dm -> df}]]

And:

Eigenvalues[A /. {bm -> bf, dm -> df}]

Eigenvectors[A/.{bm -> bf, dm -> df}]

When you go with
          bm  =  bf and dm  =  df, 
as above, and you go with 
          2 bf  >  df,

which eigenvector is the ultimate eigenvector in the sense that as t gets large
{male[t],female[t]} gravitates to the line through {0,0} determined by this eigenvector?

Remember that bm,  bf,  dm,  and df are all positive.
And because populations cannot be negative, it is impossible to start on a multiple of {-1,1}.

G.6.b.ii)  Ultimate sex ratio = 1 if  2 bf  > df

Explain: 
When you go with
          bm  =  bf and dm  =  df, 
and you go with 
           2 bf  >  df,
then this model predicts that both sexes prosper, growing approximately exponentially, with ultimate sex ratio
         male[t]/female[t] -> 1  as  t -> ∞.
So that any huge initial excess of males or females disappears over the course of time regardless of starting data.

G.6.b.iii) Extinction if  2 bf  < df

Explain:
 When you go with
          bm  =  bf and dm  =  df, 
and you go with 
          2 bf  <  df,
then this model predicts that both sexes are headed towards extinction in the sense that
    male[t] -> 0   and   female[t] -> 0  as  t -> ∞.

G.6.c.i) When females are population dominant

In certain populations, the females are population dominant in the sense that the growth of the populations of both sexes depends only on the number of females. Examples are the wild horse population and some pure-bred dog populations in which very few males are allowed to breed while females are bred almost indiscriminately. 
Discuss why an appropriate linear population model for a female population dominant society is          
          male'[t]  =  - dm male[t]  +  bm female[t]
          female'[t]  =  - df female[t]  +  bf  female[t].

G.6.c.ii)  Female population dominant society with bf  >  df

Go with the linear population model for a female population dominant society:
          
          male'[t]  =  - dm male[t]  +  bm female[t]
          female'[t]  =  - df female[t]  +  bf  female[t].

Set up the matrix:

Clear[dm, bf, dm, df];
A= (({{-dm, bm}, {0, bf - df}}));
MatrixForm[
      A]

Check:

Clear[male, female, t];
linpopmodel={male'[t],female'[t]}==A . {male[t], female[t]};

ColumnForm[Thread[linpopmodel]]

Now do as you are compelled:

Clear[eigenvalue];
{eigenvalue[1], eigenvalue[2]} = Eigenvalues[A]

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

If bf  >  df, then which eigenvector is the ultimate eigenvector ?
How do you know that if bf  >  df, then this model predicts that the both male and female populations eventually grow approximately exponentially and the ultimate sex ratio
         male[t]/female[t] -> (((dm  +  bf  -  df))/bm) as t -> ∞ ?
Does this allow for the possibility that the male population tends to 0 because very few of them are needed to sustain the overall populations?

Remember that bm,  bf,  dm, and df are all positive.

G.6.c.iii) Controlling the ultimate sex ratio in a fascist female population dominant society

Imagine that you are a big shot in a fascist female population dominant society and imagine that you can control the size of dm (percentage death rate of males). How would you set dm in terms of bf, df and bm to guarantee that the ultimate sex ratio 
         male[t]/female[t] -> 1 as t -> ∞?

G.6.c.iv) What the model predicts if   bf  <  df

Look again at the linear population model for a female population dominant society:
          male'[t]  =  - dm male[t]  +  bm female[t]
          
      female'[t]  =  - df female[t]  +  bf  female[t].

Set up the matrix:

Clear[dm, bf, dm, df];
A=(({{-dm, bm}, {0, bf - df}}));
MatrixForm[A]

Check:

Clear[male, female, t];
linpopmodel={male'[t],female'[t]}==A . {male[t], female[t]};

ColumnForm[Thread[linpopmodel]]

Now do as you are compelled:

Clear[eigenvalue];
{eigenvalue[1], eigenvalue[2]} = Eigenvalues[A]

What does the model predict if bf  <  df?

G.7) Discrete dynamical systems: Life Sciences
Population dynamics and Mendelian genetic dynamics

G.7.a) Population dynamics

Agree that
 x[n] = the fraction of the human population living in urban areas  in year 10 n
 and

y[n]  is the fraction of the population living in rural areas in year 10 n;

so x[n + 1] is the fraction of the population living in urban areas and y[n + 1] is the fraction of the population living in rural areas in year 10 years later. 

The relationship between
      {x[n], y[n]}  and {x[n + 1], y[n + 1]} 
 is expressed in this discrete dynamical system:
       x[n + 1] = 0.98 x[n]  +  0.12 y[n]
       y[n + 1] = 0.02 x[n]  +  0.88 y[n] .
       
Interpretation: 
▪ Over a ten year period, 98% of the urban dwellers stay in urban areas and 12% of the rural dwellers move to urban areas. And
▪Over a ten year period, 2% of the urban dwellers move to rural areas and 88% of the rural dwellers continue to live in rural areas.

Set up a good matrix and calculate its eigenvalues and eigenvectors. 
Then use the results to calculate the ultimate ratio  x[n]/y[n] as n gets large.
In the long run, what percentage of the population lives in urban areas?

G.7.b.i) Discrete dynamical systems: Mendelian genetic dynamics:
Continued continued crossings with a hybrid

Gregor Mendel (1822-84)  laid the base for the modern theory of genetics by putting it on a mathematical basis.
In the Mendelian system, a gene pair is:
        ▪  dominant (D) if both of its slots are dominant,
        ▪  recessive (R) if both of its slots are recessive or 
        ▪  hybrid (H) if it has one  dominant and one recessive slot.

Here are the predicted outcomes according to Mendelian theory:
                Parents ------> Offspring
                D - D-------->  D 100%
                D - H --------> D 50% and H 50%
	     D - R---------> H 100%
	     H - H---------> D 25% , H 50%  and R 25% 
	     H - R---------> H 50%  and R 50 % 
	     R - R ---------> R 100%
Seedgrowers often practice the procedure known as continued crossings with a hybrid (H).
This means each new plant has one hybrid parent  and one parent from the general population.
Generally the hybrid parent has a certain desirable quality.
Start the process of continued crossings with a hybrid and put
Put 
      D[n] = fraction of dominants in generation n of the process,
      H[n] = fraction of hybrids in generation n of the process, and
      R[n] = fraction of recessives in generation n of the process.
This gives:
              D[n + 1]=0.5 D[n]+0.25 H[n] + 0.0 R[n]

H[n + 1]=0.5 D[n]+0.5 H[n]+0.5 R[n]

R[n + 1]=0.0 D[n]+0.25 H[n]+0.5 R[n]

Summary:
Go with  continued crossings with a hybrid, and put 
      D[n] = fraction of dominants in generation n.
      H[n] = fraction of hybrids in generation n.
      R[n] = fraction of recessives in generation n.
 This leads to the discrete dynamical system:
             D[n + 1] = 0.5 D[n] + 0.25 H[n] + 0.0 R[n],
              H[n + 1] = 0.5 D[n] + 0.5 H[n] + 0.5 R[n],
             R[n + 1] = 0.0 D[n] + 0.25 H[n] + 0.5 R[n].
Set up a good matrix A and calculate its eigenvalues and eigenvectors. 
Then use the results to predict the limiting values of  D[n], H[n] , and R[n] as n gets large.

G.7.b.ii) Mendelian genetic dynamics: Starting with a hybrid-dominant mating

Look at these calculations which show (A^k) . {1, 0, 0} for generations k = 0,1,2,...,18.

A = (({{0.5, 0.25, 0}, {0.5, 0.5, 0.5}, {0, 0.25, 0.5}}));
starter  =  {1, 0, 0};
ColumnForm[
      Table[{k , MatrixPower[A, k] . starter}, {k, 0, 18}]]

That 
         starter = {1, 0, 0}
says that in generation 0,the hybrids were always mated to a  dominant D.
Are there any recessives in the first generation?
What percentage of recessives are expected in generation 3?

G.7.b.iii) Mendelian genetic dynamics: Continued crossings with a dominant

Go back to:
               Parents ------> Offspring
                D - D-------->  D 100%
                D - H --------> D 50% and H 50%
	     D - R---------> H 100%
	     H - H---------> D 25% , H 50%  and R 25% 
	     H - R---------> H 50%  and R 50 % 
	     R - R ---------> R 100%
This time investigate the procedure known as continued crossings with a dominant (D).
This means each new plant has one D parent  and one parent from the general population.
     D[n] = fraction of dominants in generation n of the process,
     H[n] = fraction of hybrids in generation n of the process, and
     R[n] = fraction of recessives in generation n of the process.
Write down the resulting dynamical system. Set up a good matrix and calculate its eigenvalues and eigenvectors. 
Then use the results to calculate to predict the limiting values of  D[n], H[n] , and R[n] as n gets large.

Look at:
                D - D-------->  D 100%
                D - H --------> D 50% and H 50%
	     D - R---------> H 100%
	     H - H---------> D 25% , H 50%  and R 25% 
	     H - R---------> H 50%  and R 50 % ,
	     R - R ---------> R 100%
 For one thing,
            D[n + 1]=1.0 D[n]+0.50 H[n]+0.0 R[n]. 
Also
             H[n+1] = ????????? (You fill this in.)
And
              R[n + 1] = 0.0 D[n] + 0.0 H[n] + 0.0 R[n].

G.8) The Fibonacci difference equation

G.8.a.i) Fibonacci numbers come  from the difference equation
		y[k + 2] = y[k + 1] +  y[k]

You get the Fibonacci numbers from the difference equation
		y[k + 2] = y[k + 1] +  y[k], starting with y[0] = 0 and y[1] = 1.

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, . . ..

G.8.a.ii) Slamming out some Fibonacci numbers

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[49] and y[50]
in terms of y[0]  =  0 and y[1]  =  1.

G.8.a.iii)  The golden ratio

Check the eigenvalues and eigenvectors of your coefficient matrix A and use the result to predict the limiting value
of y[k + 1]/y[k]as k gets large.

N[GoldenRatio]

G.8.a.iv) Different values of y[0] and y[1]

When you go with the Fibonacci difference equation
y[k + 2] = y[k + 1] +  y[k] and start with y[0] > 0 and y[1] > 0 
(so that y[0] is not necessarily 0 and y[1] is not necessarily 1), 
then what is the limiting value of y[k + 1]/y[k]as k gets large?

G.9) Using the cheat sheet to analyze the behavior of  (E^(A t)) . starter

G.9.a.i)

Here is the cheat sheet chart from the Tutorials:

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];\)
\)
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];

To analyze a matrix A = (({{a, b}, {c, d}})),
you put
           trace   =  a  +  d
and 
           det  =  a d  -  b c 
and then plot the point {trace,det} on the chart.

G.9.a.i) Going to {0,0}

Take another look at the chart.

Show[chart];

Agree or disagree:
If  matrix A = (({{a, b}, {c, d}})) is a given matrix and
       trace   =  a  +  d < 0,
and
       det  =  a d  -  b c > 0,
then 
       (E^(A t)) . starter -> {0, 0}
no matter where the starter point is.

G.9.a.ii) Going to infinity

Take another look at the chart.

Show[chart];

Agree or disagree:
If  matrix A = (({{a, b}, {c, d}})) is a given matrix and
           trace   =  a  +  d > 0,
and
            det  =  a d  -  b c > 0,
then 
            ||(E^(A t)) . starter|| -> infinity
no matter where the starter point is (as long as starter ≠ {0,0}).

G.9.a.iii) Oscillating on closed curves if trace = 0

Take another look at the chart.

Show[chart];

Agree or disagree:
If  matrix A = (({{a, b}, {c, d}})) is a given matrix and
           trace   =  a  +  d = 0,
then as t advances
            (E^(A t)) . starter
oscillates on closed curves centered at {0,0}.

G.9.a.iv) Set a constant

Take another look at the chart.

Show[chart];

Here is a matrix A with one cleared entry:

Clear[a];
A = (({{a, -2}, {1, 1}}));
MatrixForm[A]

This time you have three jobs:
Go with starter  ≠ {0,0}.
▪Job 1: Set a so that (E^(A t)) . starter -> {0, 0} as t gets large.
Put answer here:

▪Job 2: Set a so that ||(E^(A t)) . starter|| -> ∞ as t gets large.
Put answer here:

▪Job 3: Set a so that as t advances (E^(A t)) . starter oscillates on closed curves centered at {0,0}.

G.9.b.i) Round off instability

Here's a new matrix A:

A=(({{-1.2, 1.5}, {-1.3, 1.2}}));
MatrixForm[A]

Check the eigenvalues of A:

Clear[eigenvalue];
{eigenvalue[1], eigenvalue[2]} = Eigenvalues[A]

The form is p  +  I q with 
        p  =  0 and q ≠ 0. 
This means as t advances (E^(A t)) . starter is guaranteed to advance on an ellipse centered at {0,0}.

ranger = 10;
\(\(starter = {Random[Real, {\(-\(ranger\/4\)\), ranger\/4}], 
          Random[Real, {\(-\(ranger\/4\)\), ranger\/4}]};\)
\)
MatrixExp[A t] . starter

\(\(endtime = 24;\)
\)
\(\(trajectoryplot = 
        ParametricPlot[MatrixExp[A\ t] . starter, {t, 0, endtime}, 
          PlotStyle -> {{CadmiumOrange, Thickness[0.01]}}, 
          DisplayFunction -> Identity];\)
\)
starterpoint = starter;
\(\(starterplot = 
        Graphics[{Red, PointSize[0.06], Point[starterpoint]}];\)
\)
Show[starterplot, trajectoryplot, 
        PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
        Axes -> True, AxesLabel -> {"x", "y"}, 
        DisplayFunction -> $DisplayFunction];

Now take the coefficient matrix A and build in some random noise like this and compare:

h = 0.049;
radomnoise = {{Random[Real, {-h, h}], 
            Random[Real, {-h, h}]}, {Random[Real, {-h, h}], 
            Random[Real, {-h, h}]}};
\(\(noisyA = A + radomnoise;\)
\)
MatrixForm[noisyA]
MatrixForm[A]

NoisyA is guaranteed to round off nicely to A.
Now look at a plot of 
         (E^(noisyA t)) . starter 
 together with a plot of 
           (E^(A t)) . starter:

h = 0.049;
randomnoise = {{Random[Real, {-h, h}], 
            Random[Real, {-h, h}]}, {Random[Real, {-h, h}], 
            Random[Real, {-h, h}]}};
\(\(noisyA = A + randomnoise;\)
\)
\(\(noisytrajectoryplot = 
        ParametricPlot[MatrixExp[noisyA\ t] . starter, {t, 0, endtime}, 
          PlotStyle -> {{NavyBlue, Thickness[0.01]}}, 
          DisplayFunction -> Identity];\)
\)
\(\(Show[starterplot, noisytrajectoryplot, trajectoryplot, 
        PlotRange -> {{\(-ranger\), ranger}, {\(-ranger\), ranger}}, 
        Axes -> True, AxesLabel -> {"x", "y"}, 
        DisplayFunction -> $DisplayFunction];\)
\)
MatrixForm[A]
MatrixForm[noisyA]

Describe what you see.
Comment on the statement: 
→ Even though 
1) the eigenvalues of A are of the form p + I q and p - I q with p = 0  and  q ≠ 0   and 
2) noisyA rounds off to A, it is quite unlikely that the eigenvalues of noisyA are of the form  are of the form p + I q and p - I q with p = 0  and  q ≠ 0.

G.9.b.ii) Explanation

One of the goals of mathematics is to explain why things work out the way they do.
To this end, take anther look at look at the cheat sheet chart from one of the Tutorials.

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];\)
\)
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];

Given a coefficient matrix:

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

The part of the chart corresponding eigenvalues of the form  p + I q and p - I q with p = 0  and  q ≠ 0 is 
        trace  =  a  +  d  =  0
and 
        det   =  a d  -  bc  >  0.

Try to use this information to explain:
If you start with any coefficient matrix A whose eigenvalues are of the  form  p + I q and p - I q with p = 0  and  q ≠ 0, and you make the slightest random error in entering A, then the resulting matrix is all but guaranteed to have eigenvalues not of the form that  p + I q and p - I q with p = 0  and  q ≠ 0.

G.9.b.iii) Other possibilities

Take another look at the chart:

Show[chart];

If you are given a matrix with 
            trace  <  0 and det  =  0,
then use the chart to predict what wildly different things can happen when very slight random errors are made when you enter the matrix.

G.10) Matrix exponential (E^(A t)) for matrices A with
           eigenvalue[1] = 0 and 0 > eigenvalue[2].
 Matrix powers (A^k) for matrices A with
          eigenvalue[1] = 1 and 1> | eigenvalue[2] |

G.10.a) Matrix exponential (E^(A t)) for matrices A with eigenvalue[1] = 0 and 0 > eigenvalue[2].

Here are facts from the Basics:

If A is any diagonalizable matrix and the absolute values all eigenvalues of A are negative , then no matter what X is, you are guaranteed that 
          (E^(A t)) . X -> {0, 0}as t gets large.

If A is any diagonalizable matrix and the absolute values of the of both eigenvalues of A are both positive, then no matter what as long as X ≠ {0, 0}, you are guaranteed that 
          (((||E))^(A t)) . 
            X||->infinity as t gets large.

If eigenvalue[1] > eigenvalue[2], then the ultimate direction of (E^(A t)) . X as t gets large is in the direction of eigenvector[1] (unless X is a multiple of eigenvector[2])

This does not give info about the case in which
           eigenvalue[1] =  0 and 0  > eigenvalue[2].
It is possible to predict exactly what will happen. If this intrigues you, go on.

Here's a  2D matrix A  with eigenvalue[1] = 0 and 0  > eigenvalue[2]

A = (({{-0.50, 0.40}, {0.25, -0.20}}));	

MatrixForm[
      A]

Check the eigenvalues and eigenvectors:

Eigenvalues[A]

Enter them so that  eigenvalue[1] = 0 and 0 >eigenvalue[2]

Clear[eigenvalue];
{eigenvalue[2], eigenvalue[1]}  =  Eigenvalues[A]

Enter the eigenvectors the corresponding way so that
         A . eigenvector[1] = eigenvalue[1] eigenvector[1]
and 
         A . eigenvector[2] = eigenvalue[2] eigenvector[2] :

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

Here are four random points
        {starter[1], starter[2], starter[3], starter[4]} 
in 2D shown with scaled versions of the eigenvectors of A:

\(\(Clear[trajectoryplot, eigenplot, scaler, starter, thigh, k, 
        j];\)
\)
\(\(trajectoryplot[thigh_, starter_]\  := \ 
        ParametricPlot[MatrixExp[A\ t] . starter, {t, 0, thigh}, 
          PlotStyle -> {{Thickness[0.01], CadmiumOrange}}, 
          DisplayFunction -> Identity];\)
\)
\(\(starterplots\  = \ 
        Table[Graphics[{NavyBlue, PointSize[0.03], Point[starter[j]]}], {j, 
            1, 4}];\)\ 
\)
\(\(scaler\  = \ 12;\)\ 

\)
\(\(eigenlabels\  = \ {Graphics[{Red, 
              Text["eigenvector[1]", 
                0.6\ scaler\ eigenvector[1]]}], 
		Graphics[{GrayLevel[
                0.5], Text["eigenvector[2]", 
                0.6\ scaler\ \ eigenvector[2]]}]};\)
\)
\(\(eigenplot[
          j_]\  := \ 
\ \ \ \ \ \ \ \ \ {Arrow[scaler\ eigenvector[j], 
            Tail -> {0, 0}, 
            VectorColor -> 
              Black], 
\ \ \ \ \ \ \ \ \ \ \ Arrow[\(-scaler\)\ \ eigenvector[j], Tail -> {0, 0}, 
            VectorColor -> Black]};\)
\)
starter[1] =  
        Random[Real, {4, 10}] eigenvector[1]  +  
          Random[Real, {4, 10}] eigenvector[2];
starter[2] =  
        Random[Real, {-10, -4}] eigenvector[1]  + 
          Random[Real, {-10, -4}] eigenvector[2];
starter[3] =  
        Random[Real, {-10, -4}] eigenvector[1]  + 
          Random[Real, {4, 10}] eigenvector[2];
\(\(starter[4] = \ 
        Random[Real, {4, 10}]\ eigenvector[1]\  + \ 
          Random[Real, {\(-10\), \(-4\)}] eigenvector[2];\)
\)
thigh  =  0.01;
\(\(Show[Table[trajectoryplot[thigh, starter[j]], {j, 1, 4}]\ , 
        starterplots, eigenplot[1], eigenplot[2], eigenlabels, 
        PlotLabel -> "Starting points with scaled eigenvectors", 
        DisplayFunction -> $DisplayFunction];\)
\)

Here is how the curves
 (E^(A t)) . starter[1], (E^(A t)) . starter[2], 
      (E^(A t)) . starter[3], (E^(A t)) . starter[4]
plot out as t advances from 0 to 0.5 shown with scaled eigenvectors of A:

thigh  =  0.5;
Show[Table[trajectoryplot[thigh, starter[j]], {j, 1, 4}] ,starterplots,eigenplot[1],eigenplot[2],eigenlabels,PlotLabel->thigh "\!\(\*
StyleBox[\"is\",\nFontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\" \",\nFontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\"thigh\",\nFontColor->RGBColor[0, 0, 1]]\)",DisplayFunction -> $DisplayFunction];

Here's what happens as t advances from 0 to 1.5 shown with scaled eigenvectors of A:

thigh  =  1.5;
Show[Table[trajectoryplot[thigh, starter[j]], {j, 1, 4}] ,starterplots,eigenplot[1],eigenplot[2],eigenlabels,PlotLabel->thigh "\!\(\*
StyleBox[\"is\",\nFontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\" \",\nFontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\"thigh\",\nFontColor->RGBColor[0, 0, 1]]\)",DisplayFunction -> $DisplayFunction];

Here's what happens as t advances from 0 to 3 shown with scaled eigenvectors of A:

thigh  =  3;
Show[Table[trajectoryplot[thigh, starter[j]], {j, 1, 4}] ,starterplots,eigenplot[1],eigenplot[2],eigenlabels,PlotLabel->thigh "\!\(\*
StyleBox[\"is\",\nFontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\" \",\nFontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\"thigh\",\nFontColor->RGBColor[0, 0, 1]]\)",DisplayFunction -> $DisplayFunction];

Here's what happens as t advances from 0 to 5 shown with scaled eigenvectors of A:

thigh  =  5;
Show[Table[trajectoryplot[thigh, starter[j]], {j, 1, 4}] ,starterplots,eigenplot[1],eigenplot[2],eigenlabels,PlotLabel->thigh "\!\(\*
StyleBox[\"is\",\nFontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\" \",\nFontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\"thigh\",\nFontColor->RGBColor[0, 0, 1]]\)",DisplayFunction -> $DisplayFunction];

Here's what happens as t advances from 0 to 20:

thigh  = 20;
Show[Table[trajectoryplot[thigh, starter[j]], {j, 1, 4}] ,starterplots,eigenplot[1],eigenplot[2],eigenlabels,PlotLabel->thigh "\!\(\*
StyleBox[\"is\",\nFontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\" \",\nFontColor->RGBColor[0, 0, 1]]\)\!\(\*
StyleBox[\"thigh\",\nFontColor->RGBColor[0, 0, 1]]\)",DisplayFunction -> $DisplayFunction];

As k gets large, the points (E^(A t)) . starter stall on the line through {0,0} defined by eigenvector[1].
You could have predicted this in advance because  eigenvalue[1] = 0 and 0 > eigenvalue[2]:.
What is different is: 
-> as k gets large, (E^(A t)) . starter does not go to {0, 0} and
  -> as k gets large,  (||(E^(A t)) . starter|| does not go to infinity. 
    
Instead (E^(A t)) starter gravitates to a certain specific point on the line through {0,0} defined by eigenvector[1].
The fact of the matter is that if 
       eigenvalue[1] = 0 and 0 > eigenvalue[2]:. 
and 
        starter = s eigenvector[1] + t eigenvector[2] 
then you are guaranteed that as t gets large
          (E^(A t)) . starter -> s eigenvector[1].
Your job is to explain why this is guaranteed.

In one of the Basics, this fact is explained:
 \(E\^\(A\ t\) . \((s\ eigenvect[1]\  + \ 
              t\ eigenvect[
                  2])\)

\ \ \ \ \ \  = \(s\ E\^\(eigenval[1]\ t\)\ eigvect[
                1] + t\ E\^\(eigenval[2]\ t\)\ eigvect[
                2]
	
		\  = \ \(E\^\(eigenval[1]\ t\)\) \((s\ \ eigvect[
                    1] + t\ E\^\(\((eigenval[2]\  - \ eigenval[1])\)\ t\)\ eigvect[2]\)\)\)
In this setup eigval[1] = 0 and 0 > eigval[2]. 
This gives
 \(E\^\(A\ t\) . \((s\ eigenvect[1]\  + \ 
              t\ eigenvect[2])\)
		
			\  = \ 
        s\ \ eigvect[1] + t\ E\^\(\((eigenval[2])\)\ t\)\ eigvect[2]\)
 Use the the fact that  eigval[2] < 0  to make your explanation.

G.10.b.i) Matrix powers (A^k) for matrices A with eigenvalue[1] = 1 and 1> | eigenvalue[2] | :

Here are facts freom the Basics:

If A is any diagonalizable matrix and the absolute values all eigenvalues of A are less than 1, then no matter what X is, you are guaranteed that 
          (A^k) . X -> {0, 0}as k gets large.

If A is any diagonalizable matrix and the absolute values of the of both eigenvalues of A are both bigger than 1, then no matter what as long as X ≠ {0, 0}, you are guaranteed that 
          (((||A))^k) . 
            X||->infinity as k gets large.

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])

This does not give info about the case in which
           eigenvalue[1] = 1 and 1> | eigenvalue[2] | :
It is possible to predict exactly what will happen. If this intrigiues you, go on.

Here's a random 2D matrix A  with eigenvalue[1] = 1 and 1> | eigenvalue[2] |

a  =  Random[Real, {0.1, 0.99}];
b  = 
        Random[Real, {0.1, 0.99}];
c  =  
        1 -  a;
d  =  1  -  b ;
A = (({{a, b}, {c, d}}));

MatrixForm[
      A]

Check the eigenvalues:

Clear[eigenvalue];
{eigenvalue[1], eigenvalue[2]}  =  Eigenvalues[A]

This gives  eigenvalue[1] = 1 and 1> | eigenvalue[2] | :
Here are two random points
        {starter[1], starter[2]} 
in 2D shown with scaled versions of the eigenvectors of A:

\(\(Clear[iterationplot, interationplotter, eigenplot, scaler, point, 
        pointcolor, eigenvector, eigenlabels, starter, khigh, k, j, 
        iterationpoints];\)
\)
point[k_, starter_]  := MatrixPower[A, k]  . starter ;
\(\(\(pointcolor[k_, khigh_]\)\(\ \)\(:=\)\(\ \)\(RGBColor[
        Sin[\((Pi/2)\)\ k/\((khigh + 0.1)\)], 
        Cos[\((Pi/2)\)\ k/\((khigh\  + \ 0.1)\)], 
        E^\((\(-k\)/\((khigh + 0.1)\))\)]\)\(\ \)\(
\)\)\)
\(\(iterationpoints[starter_, khigh_]\  := 

        Table[Graphics[\ \ {pointcolor[k, khigh], PointSize[0.03], 
              Point[MatrixPower[A, k]\  . starter]}], {k, 0, 
            khigh}];\)
\)
\(\(iterationplotter[khigh_] := 
        Show[Table[iterationpoints[starter[j], khigh], {j, 1, 2}], 
          Axes -> True, AxesLabel -> {"x", "y"}, 
          PlotLabel -> \*""Hits with A,(A^2),..,(A^khigh)"", 
		AspectRatio -> 1/GoldenRatio, 
          DisplayFunction -> Identity];\)
\)
\(\(scaler := \ 
        0.9\ Sqrt[
            point[1, starter[1]] . point[1, starter[1]]\ ];\)
\)
{eigenvector[1], eigenvector[2]}  =  Eigenvectors[A];
\(\(eigenplot[
          j_] := \ 
\ \ \ \ \ \ \ \ \ {Arrow[scaler\ eigenvector[j], 
            Tail -> {0, 0}, 
            VectorColor -> 
              Black], 
\ \ \ \ \ \ \ \ \ \ \ Arrow[\(-scaler\)\ \ eigenvector[j], Tail -> {0, 0}, 
            VectorColor -> Black]};\)
\)
\(\(eigenlabels\  = \ {Graphics[{Red, 
              Text["eigenvector[1]", 
                0.6\ scaler\ eigenvector[1]]}], 
		Graphics[{GrayLevel[
                0.5], Text["eigenvector[2]", 
                0.6\ scaler\ \ eigenvector[2]]}]};\)
\)
starter[1] =  
        Random[Real, {4, 10}] eigenvector[1]  +  
          Random[Real, {4, 10}] eigenvector[2]; 
\(\(starter[2] = \ 
        Random[Real, {\(-10\), \(-4\)}]\ eigenvector[1]\  + \ 
          Random[Real, {\(-10\), \(-4\)}]\ eigenvector[2];\)\ 
\)
khigh  =  0;
\(\(Show[iterationplotter[khigh], eigenplot[1], eigenplot[2], 
        eigenlabels, 
        PlotLabel -> "Four starting points with scaled eigenvectors", 
        DisplayFunction -> $DisplayFunction];\)
\)

Here are the same two starter points together with the points you get when you hit each of the original starting points with A:

khigh  =  1;
Show[iterationplotter[khigh], eigenplot[1], eigenplot[2], 
        eigenlabels, DisplayFunction -> $DisplayFunction];

Here are the same two starter points together with the points you get when you hit each of the original starting points with A and A.A:

khigh  =  2;
Show[iterationplotter[khigh], eigenplot[1], eigenplot[2], 
        eigenlabels, DisplayFunction -> $DisplayFunction];

Here are the same two starter points together with the points you get when you hit each of the original starting points with A, A.A and A.A.A:

khigh  =  3;
Show[iterationplotter[khigh], eigenplot[1], 
        eigenplot[2], eigenlabels, 
        DisplayFunction -> $DisplayFunction];

Here are the same two starter points together with the points you get when you hit each of the original starting points with A, A . A, (A^3), (A^4), and (A^5):

khigh  =  5;
Show[iterationplotter[khigh], eigenplot[1], eigenplot[2], 
        eigenlabels, 
        PlotLabel -> "Two starting points with scaled eigenvectors", 
        DisplayFunction -> $DisplayFunction];

Here are the same two starter points together with the points you get when you hit each of the original starting points with  A, A . A, (A^3), (A^4), . . ., (A^9) , 
      (A^10):

khigh  =  10;
Show[iterationplotter[khigh], eigenplot[1], 
        eigenplot[2], eigenlabels, 
        DisplayFunction -> $DisplayFunction];

As k gets large, the points (A^k) .starter gravitate to the line through {0,0} defined by eigenvector[1].
You could have predicted this in advance because  eigenvalue[1] = 1 and 1> | eigenvalue[2] | :.
What is different is: 
-> as k gets large, (A^k) . starter does not go to {0, 0} and
 -> as k gets large,  ||(A^k) . starter|| does not go to infinity. 
    
Instead (A^k) starter gravitates to a certain specific point on the line through {0,0} defined by eigenvector[1].
The fact of the matter is that if 
       eigenvalue[1] = 1 and 1> | eigenvalue[2] | :. 
and 
        starter = s eigenvector[1] + t eigenvector[2] 
then you are guaranteed that
          (A^k) . starter -> s eigenvector[1].
Your job is to explain why this is guaranteed.

In one of the Basics, this fact is explained:
 \(A\^k . \((s\ eigenvect[1]\  + \ 
              t\ eigenvect[2])\)
		
		 = \ \ \(eigval[
                1]\^k\) \((s\ \ eigvect[1] + \ 
              t\ \((eigval[2]\/eigval[1])\)\^k\ eigvect[2])\)\)
In this setup eigval[1] = 1 and1> | eigval[2] | . 
This gives
 \(A\^k . \((s\ eigenvect[1]\  + \ 
              t\ eigenvect[2])\)
		
		 = \ \ \((s\ \ eigvect[1] + \ 
            t\ eigval[2]\^k\ eigvect[2])\)\)
 Use the behavior of  eigval[2]^(k ) as k gets large to make your explanation.

G.10.b.ii) Identify the straight line

Stay with the same setup as in part i) so that eigenvalues of A are {eigenvalue[1], eigenvalue[2]} and 
-> eigenvalue[1] = 1 and 1> | eigenvalue[2] | . 
-> starter = s eigenvector[1] + t eigenvector[2] 
then you are guaranteed that as k gets large
           (A^k) . starter -> s eigenvector[1].
Explain this: As k gets large (A^k) . starter moves to s eigenvector[1] on a straight line. 
Identify this straight line.

G.11) More on the matrix exponential and related matters
A . (E^(A t)) = (E^(A t)) . A is guaranteed.
When you can be sure that A . B = B . A.
When you can be sure that (E^(A + B)) = (E^A) . (E^B)

G.11.a.i) Eigenvalues of  (E^A)

A 2D diagonalizable matrix is any matrix A of the form

  A =SpannerMatrix.(({{eigenvalue[1], 0}, {0, eigenvalue[2]}})).(SpannerMatrix^(-1)).

You calculate (E^A) by putting
    (E^A) =SpannerMatrix.(({{(E^eigenvalue)[1], 0}, {0, (E^eigenvalue)[2]}})).(SpannerMatrix^(-1))
   Seeing this, answer this question:
Given a diagonalizable matrix A, how are the eigenvalues of (E^A) related to the eigenvalues of A?

G.11.a.ii) Eigenvectors of  (E^A)

A 2D diagonalizable matrix is any matrix A of the form

  A =SpannerMatrix.(({{eigenvalue[1], 0}, {0, eigenvalue[2]}})).(SpannerMatrix^(-1)).

You calculate (E^A) by putting
    (E^A) =SpannerMatrix.(({{(E^eigenvalue)[1], 0}, {0, (E^eigenvalue)[2]}})).(SpannerMatrix^(-1))
Seeing this, answer this question:
Given a diagonalizable matrix A, how are the eigenvectors of (E^A) related to the eigenvectors of A?

G.11.b.i)    A . (E^(A t)) = (E^(A t)) . A

Here are two random 2D matrices A and B:

A = (({{Random[Real, {-2, 2}], Random[Real, {-2, 2}]}, {Random[Real, {-2, 2}], Random[Real, {-2, 2}]}})) ;
MatrixForm[A]

B = (({{Random[Real, {-2, 2}], Random[Real, {-2, 2}]}, {Random[Real, {-2, 2}], Random[Real, {-2, 2}]}})) ;
MatrixForm[B]

Here are A.B and B.A:

MatrixForm[A . B]

MatrixForm[B . A]

As you can see, the outcome 
          A . B = B . A
is pretty darn rare.
In spite of this, if A is diagonalizable, you can always count on having
       A . (E^(A t)) = (E^(A t)) . A.
Here's a start on the idea in 2D: The same idea (with more writing) works in all dimensions.
A 2D diagonalizable matrix is any matrix A of the form

  A =SpannerMatrix.(({{eigenvalue[1], 0}, {0, eigenvalue[2]}})).(SpannerMatrix^(-1))
 
 Just to keep everything from bouncing off the screen, use this shorthand:
            SM = SpannerMatrix
            e[1] = eigenvalue[1]
            e[2] = eigenvalue[2].
            
In this shorthand, a  2D diagonalizable matrix A is of the form
             A = SM.(({{e[1], 0}, {0, e[2]}})).(SM^(-1)).
This gives:
            A t = t SM.(({{e[1], 0}, {0, e[2] }})).(SM^(-1))=  SM.(({{t e[1], 0}, {0, t e[2] }})).(SM^(-1))
         
So.
             (E^(A t))=SM.(({{(E^(t e[1])), 0}, {0, (E^(t e[2]))}})).(SM^(-1)).
And so
  A . (E^(A t)) = SM.(({{e[1], 0}, {0, e[2]}})).(SM^(-1)).SM.(({{(E^(t e[1])), 0}, {0, (E^(t e[2]))}})).(SM^(-1))
 while
  (E^(A t)) . A = SM.(({{(E^(t e[1])), 0}, {0, (E^(t e[2]))}})).(SM^(-1)).SM.(({{e[1], 0}, {0, e[2]}})).(SM^(-1)).
  
  Your job is to take over and finish this explanation of why
               A . (E^(A t)) = (E^(A t)) . A.

G.11.b.ii)  Simultaneously diagonalizable

Given a diagonalizable 2D matrix A, then you are guaranteed a basis of 2D 
           {eigenvectA[1], eigenvectA[2]}
of unit eigenvectors of A.
Given another diagonalizable 2D matrix B, then you are guaranteed a basis of 2D 
            {eigenvectB[1], eigenvectB[2]} 
of unit eigenvectors of B.

Old matrix folks know that if A and B are simultaneously diagonalizable in the sense that
           eigenvectA[1] = eigenvectB[1]
and 
          eigenvectA[2] = eigenvectB[2],
then it is certain that
          A . B = B . A.
Your job is to explain why this is certain.

A =SpannerMatrixA.(({{eigenvalA[1], 0}, {0, eigenvalA[2]}})).(SpannerMatrixA^(-1))

B =SpannerMatrixB.(({{eigenvalB[1], 0}, {0, eigenvalB[2]}})).(SpannerMatrixB^(-1))

Because
           eigenvectA[1] = eigenvectB[1]
and 
           eigenvectA[2] = eigenvectB[2],
you are certain that
           SpannerMatrixA = SpannerMatrixB. 
Just for shorthand put 
           SM = SpannerMatrixA = SpannerMatrixB. 
In this shorthand,
           A =SM.(({{eigenvalA[1], 0}, {0, eigenvalA[2]}})).(SM^(-1))
and
           B =SM.(({{eigenvalB[1], 0}, {0, eigenvalB[2]}})).(SM^(-1)).
           
 Now you take over and finish the explanation of why A . B = B . A.

G.11.b.iii)   If A is diagonalizable, then A and  (E^(A t)) are simultaneously diagonalizable

Given a diagonalizable 2D matrix A, explain why   A and  (E^(A t)) are simultaneously diagonalizable.

G.11.b.iv) Calculus Cal says  (E^(A  +  B)) = (E^A) . (E^B)

One day while you are sitting at your favorite machine, that repulsive Calculus Cal comes in and says:"Given two 2D matrices A and B, you can be sure that (E^(A  +  B)) = (E^A) . (E^B)."
 You say:"Where did you get that idea?" 
 Cal says:"For numbers a and b, everyone knows that (E^(a  +  b)) = (E^a) (E^b) is a basic identity and because matrices are made of numbers, it follows that (E^(A  +  B)) = (E^A) . (E^B) for any 2D matrices A and B." 
 You decide to check this out:

A = (({{Random[Real, {-2, 2}], Random[Real, {-2, 2}]}, {Random[Real, {-2, 2}], Random[Real, {-2, 2}]}}));
MatrixForm[A]

B =  (({{Random[Real, {-2, 2}], Random[Real, {-2, 2}]}, {Random[Real, {-2, 2}], Random[Real, {-2, 2}]}}));
MatrixForm[B]

Here is (E^(A  +  B)):

MatrixForm[MatrixExp[A  +  B]]

Here is   (E^A) . (E^B):

MatrixForm[MatrixExp[A].MatrixExp[B]]

Rerun a couple of times and then tell Cal what you think of his theory.

G.11.b.v) If A and B are simultaneously diagonalizable, then    (E^(A + B)) = (E^A) . (E^B).

Given a diagonalizable 2D matrix A, then you are guaranteed a basis of 2D 
           {eigenvectA[1], eigenvectA[2]}
of unit eigenvectors of A.
Given another diagonalizable 2D matrix B, then you are guaranteed a basis of 2D 
            {eigenvectB[1], eigenvectB[2]} 
of unit eigenvectors of B.

Old matrix folks know that A and B are simultaneously diagonalizable in the sense that
           eigenvectA[1] = eigenvectB[1]
and 
          eigenvectA[2] = eigenvectB[2],
then it is certain that
          (E^(A + B)) = (E^A) . (E^B).
Your job is to explain why this is certain.

A =SpannerMatrixA.(({{eigenvalA[1], 0}, {0, eigenvalA[2]}})).(SpannerMatrixA^(-1))

B =SpannerMatrixB.(({{eigenvalB[1], 0}, {0, eigenvalB[2]}})).(SpannerMatrixB^(-1))

Because
           eigenvectA[1] = eigenvectB[1]
and 
           eigenvectA[2] = eigenvectB[2],
you are certain that
           SpannerMatrixA = SpannerMatrixB. 
Just for shorthand put 
           SM = SpannerMatrixA = SpannerMatrixB. 
In this shorthand,
           A =SM.(({{eigenvalA[1], 0}, {0, eigenvalA[2]}})).(SM^(-1))
and
           B =SM.(({{eigenvalB[1], 0}, {0, eigenvalB[2]}})).(SM^(-1)).
           In this shorthand,
           (E^A) =SM.(({{(E^eigenvalA)[1], 0}, {0, (E^eigenvalA)[2]}})).(SM^(-1))
and
           (E^B) =SM.(({{(E^eigenvalB)[1], 0}, {0, (E^eigenvalB)[2]}})).(SM^(-1)).
           
So
     (E^A) . (E^B)
= SM.(({{(E^eigenvalA)[1], 0}, {0, (E^eigenvalA)[2]}})).(SM^(-1)).SM.(({{(E^eigenvalB)[1], 0}, {0, (E^eigenvalB)[2]}})).(SM^(-1))
= SM.(({{(E^eigenvalA)[1], 0}, {0, (E^eigenvalA)[2]}})).(({{(E^eigenvalB)[1], 0}, {0, (E^eigenvalB)[2]}})).(SM^(-1))
=SM.(({{(E^eigenvalA)[1] (E^eigenvalB)[1], 0}, {0, (E^eigenvalA)[1] (E^eigenvalA)[2]}})).(SM^(-1)) = SM.(({{(E^(eigenvalA[1] + eigenvalB[1])), 0}, {0, (E^(eigenvalA[2] + eigenvalB[2]))}})).(SM^(-1)).
   Once you know that 
    A  +  B =  SM.(({{eigenvalA[1] + eigenvalB[1], 0}, {0, eigenvalA[2] + eigenvalB[2]}})).(SM^(-1)). 
 
 Then you will know that   
  (E^(A + B)) = SM.(({{(E^(eigenvalA[1] + eigenvalB[1])), 0}, {0, (E^(eigenvalA[2] + eigenvalB[2]))}})).(SM^(-1))
		=(E^A) (E^B) .    

To know that
 A  +  B =  SM.(({{eigenvalA[1] + eigenvalB[1], 0}, {0, eigenvalA[2] + eigenvalB[2]}})).(SM^(-1)),
 all you need to explain is:
 1) Why
           X = eigenvectA[1] = eigenvectB[1]
is an eigenvector of A + B and why
          (A + B) . X = (eigenvalA[1] + eigenvalB[1]) X
and
2) Why
           Y = eigenvectA[2] = eigenvectB[2]
is an eigenvector of A + B and why
          (A + B) . X = (eigenvalA[2] + eigenvalB[2]) X

Now you take over and explain why these two.