MGM.06 Beyond 3D
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].";

Needs["Statistics`ContinuousDistributions`"];

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

G.1) Linear systems ready for your analysis

G.1.a.i) A linear system

Here's a matrix:

hitdim = 3;
hangdim = 5;
A=(({{0.722, -2.855, 2.841}, {1.315, -0.229, -1.267}, {2.884, -1.006, 0.789}, {0.192, 2.081, -2.131}, {1.929, -1.171, 2.328}}));
MatrixForm[A]

And a linear system based on this matrix:

Clear[x, y, k];
X = 
        Table[x[k], {k, 1, hitdim}];
Y = {4.7308, -5.5899, -2.3665, -4.5566, 3.1896};

linearsystem=ColumnForm[Thread[A . X == Y]]

Is the coefficient matrix of full rank?
Determine whether this linear system has
-> no solution,
-> many solutions,
-> exactly one solution.
Say how you came to your conclusion.
If it has no solution, then give the best approximate solution

G.1.a.ii) A related linear system

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

MatrixForm[A]

Here's a new linear system based on the same matrix:

Y = {3.6, -2.9, 3.55, -2.037, 6.63};

linearsystem=ColumnForm[Thread[A . X == Y]]

Determine whether this linear system has
-> no solution,
-> many solutions,
-> exactly one solution.
Say how you came to your conclusion. 
If it has no solution, then give the best approximate solution

G.1.b) A new linear system

Here's a matrix:

hitdim = 6;
hangdim = 5;

Clear[x, y, k];
A=(({{0.722, -2.855, 2.841, 0.341, -2.395, 0.888}, {1.315, -0.229, -1.267, 0.582, 1.953, -0.011}, {2.884, -1.006, 0.789, 0, -2.372, 2.222}, {0.192, 2.081, -2.131, 5.216, -3.333, 0.799}, {1.929, -1.171, 2.328, 0, 4.998, -2.007}}));

MatrixForm[A]

And a linear system based on this matrix:

X = Table[x[k], {k, 1, hitdim}];
\(\(Y = {\(-12.474\), 12.6209, 0.1854, 4.2912, 13.2443};\)
\)
ColumnForm[Thread[  A . X == Y]]

Is the coefficient matrix of full rank?
Determine whether this linear system has
-> no solution,
-> many solutions,
-> exactly one solution.
Say how you came to your conclusion.
If it has no solution, then give the best approximate solution.

G.1.c) A new linear system

Here's a matrix:

hitdim = 5;
hangdim = 5;

Clear[x, y, k];
A=(({{1.0, 0, 2.0, 0, 0}, {0, 1.0, 2.0, 0, 0}, {0, -1.0, -2.0, 0, 0}, {0, 0, 0, 1, 0}, {0, -0, 0, 0, 1}}));

MatrixForm[A]

And a linear system based on this matrix:

Clear[x];
X = Table[x[k], {k, 1, hitdim}];
\(\(Y = {2.4, 3.6, \(-3.7\), 0.0, 0.6};\)
\)
ColumnForm[Thread[ A . X == Y]]

Is the coefficient matrix of full rank?
Determine whether this linear system has
-> no solution,
-> many solutions,
-> exactly one solution.
Say how you came to your conclusion.
If it has no solution, then give the best approximate solution.

G.2) Nuts and bolts

G.2.a.i) Hitdim and hangdim

Here is a matrix for your viewing pleasure:

A = (({{-0.48876834635327504`, -0.26932739192216704`, 1.779692518541546`, -0.46806326837622736`}, {1.6161806555615859`, -0.675514053657053`, -0.23639265281896038`, -1.551473351685659`}, {1.1117550505617517`, -0.7014076389160944`, 1.2666317714035773`, -0.991857835009845`}, {1.9303725741926079`, -1.511316721927285`, 0.5212590952404645`, 1.850685872694322`}, {1.136037801999222`, -1.2381966146621632`, 0.5323207293871524`, -0.41248948997274537`}, {-1.4622484207455675`, -1.1218413780792167`, -0.058440230889615474`, -1.6761964116079464`}}))

This matrix hits on hitdimD and hangs its hits in hangdimD.
Fill the blanks: For this matrix 
hitdim =............. 
hangdim = ...........
The SVD hanger frame vectors for this matrix come from  .......D.
The SVD aligner frame vectors for this matrix come from  .......D

Cols and Ytest

G.2.a.ii) Stretches

Here's a matrix A which hits on 6D and hangs in 6D:

\(\(A = {{\(-3.000\), 0.410, \(-3.390\), \(-0.105\), 3.721, 
            0.000}, 
\ \ \ \ \ \ \ \ \ {1.515, 
            2.948, \(-2.29\), \(-5.372\), \(-0.496\), 
            1.241}, 
\ \ \ \ \ \ \ \ \ \ {\(-1.224\), \(-1.130\), 3.648, 
            3.662, 2.022, \(-1.919\)}, 
\ \ \ \ \ \ \ \ \ \ \ {3.071, \(-3.157\), 0.876, 3.446, \(-1.178\), 
            3.070}, 
\ \ \ \ \ \ \ \ \ \ \ {\(-0.697\), \(-0.964\), \(-1.788\), \(-0.824\), 2.383, \(-0.685\)}, 
\ \ \ \ \ \ \ \ \ \ \ {0.696, 0.227, 
            0.673, \(-1.312\), 3.00, 2.985}};\)
\)
MatrixForm[A]

Here is a calculation of the nonzero SVD stretch factors of A:

SingularValues[A][[2]]

What is it about this setup that signals in advance that when you have this matrix A as the coefficient matrix of a linear system,A.X = Y, then you are guaranteed that there is exactly one solution vector X for each given 6D vector Y?

G.2.a.iii) Rank

Here's a matrix A which hits on 4D and hangs in 5D:

\(\(A = {{\(-2.1\), 3.8, \(-2.4\), 3.3}, 
	{6.7, \(-1.5\), 2.4, 
            3.0}, 
	{3.6, 3.4, \(-0.3\), \(-0.7\)}, 
	\ {2, 1.3, 
            3.9, \(-3.2\)}, 
	\ \ {\(-3.3\), 1.2, 
            3.5, \(-2.9\)}};\)
\)
MatrixForm[A]

Here is a calculation of the rank of A:

rank = 
      Length[SingularValues[
            A][[2]]]

What is it about this setup that signals in advance that when you have this matrix A as the coefficient matrix of a linear system, then the linear system has exactly one solution or has no solution at all.

G.2.b.i) Determinant of coefficient matrix is 0

Here is a matrix which hits on 5D and hangs in 5D:

\(\(A = {{2.686, \(-0.666\), \(-3.0464\), \(-0.1568\), \(-1.881\)}, 
\ \ \ \ \ \ 		{3.068, 2.667, 
            1.596, \(-3.594\), \(-2.986\)}, 
\ \ \ \ \ \ 		{\(-1.918\), 
            2.915, \(-3.093\), \(-0.0829\), \(-0.354\)}, 
\ \ \ \ \ \ 		{1.15, 5.582, \(-1.497\), \(-3.6769\), \(-3.34\)}, 
\ \ \ \ \ \ 		{\(-0.3820\), \(-3.333\), \(-4.6424\), 3.4372, 1.105}};\)
\)
MatrixForm[A]

And the corresponding linear system of equations:

hitdim = 5;
\(\(hangdim = 5;\)
\)
Clear[x, y, k];
X = Table[x[k], {k, 1, hitdim}];
\(\(Y = Table[y[k], {k, 1, hangdim}];\)
\)
linearsystem  =  ColumnForm[Thread[A . X  ==   Y]]

Look at this calculation:

Det[A]

Use the result to comment on the following statement:
If for a certain given
        Y = {y[1],y[2],y[3],y[4],y[5]}, 
this linear system has one solution then it has many more.

G.2.b.ii) Exhibiting a solution and exhibiting a non-solution

Stay with the same linear system as in part i) immediately above:

linearsystem

Come up with a specific example of a
         Y = {y[1], y[2], y[3], y[4], y[5]}
for which the linear system has a solution.
Put answer here.

A.{0, 0, 0, 0, 0}

Then come up with a specific example of a
         Y = {y[1], y[2], y[3], y[4], y[5]}
for which the linear system has no solution.

Try:

Y  =  Table[Random[Real, {-2, 2}], {j, 1, 5}]

G.2.c.i) More than one solution

Here is a matrix which hits on 5D and hangs in 3D:

\(\(A = {{2.686, \(-0.666\), \(-3.0464\), \(-0.1568\), \(-1.881\)}, 
\ \ \ \ \ \ 		{3.068, 2.667, 
            1.596, \(-3.594\), \(-2.986\)}, 
	\ \ \ \ \ \ 	{\(-1.918\), 
            2.915, \(-3.093\), \(-0.0829\), \(-0.354\)}};\)
\)
MatrixForm[A]

And the corresponding linear system of equations:

hitdim = 5;
\(\(hangdim = 3;\)
\)
Clear[x, y, k];
X = Table[x[k], {k, 1, hitdim}];
\(\(Y = Table[y[k], {k, 1, hangdim}];\)
\)
linearsystem  =  ColumnForm[Thread[A . X  ==   Y]]

The rank of the coefficient matrix is:

rank = 
      Length[SingularValues[
            A][[2]]]

You look at this and say that no matter what Y = {y[1],y2],y[3]}  this linear system has more than one solution. You are right. Explain why you are right.

G.2.c.ii) Coming up with the X  of minimum norm that makes A.X as close to Y as it can be

Here is yet another matrix:

A=(({{2.686, -0.666, -3.046}, {3.068, 2.667, 1.596}, {-0.382, -3.333, -4.647}, {1.035, -2.421, 0.892}, {3.817, 4.732, -3.987}}));
MatrixForm[
      A]

And a corresponding linear system of equations:

hitdim = 3; hangdim = 5;
Y = {-5.03, 3.023, -7.854, 7.780, -16.905};
Clear[x, k];
\(\(X = Table[x[k], {k, 1, hitdim}];\)
\)
linearsystem  =  ColumnForm[Thread[A . X  ==   Y]]

In effort to go after solutions of this system, you do an SVD analysis of the system and put:

rank = 
          Length[SingularValues[
                A][[2]]];

Clear[hangerframe, k];
hangerframe[k_]:=SingularValues[A][[1]][[k]];

Ytest  =  (∑_(k = 1)^rank)(Y . 
              hangerframe[k]) hangerframe[k]

And you compare:

Ytest
Y

Seeing that 
          Ytest ≠Y, 
you announce that the given linear system has no solution. Your job is to do the next best thing: 

Come up with the X = {x[1],x[2],x[3]} of minimum norm that makes A.X as close to Y as it can be.

G.2.c.ii) Plotting

Here's a matrix A (({{1.5`, -1.1`}, {-0.9`, 0}, {0.9`, 0.7`}})):

A=(({{1.5, -1.1}, {-0.9, 0}, {0.9, 0.7}}));
MatrixForm[A]

Here's a piece of the plane you get when you hit all of 2D with A:

A = {{1.5, -1.1}, {-0.9, 0}, {0.9, 0.7}};
\(\({h, k} = {6, 6};\)
\)
Clear[x, y];
planepiece = 
        ParametricPlot3D[A . {x, y}, {x, -h, h}, {y, -k, k}, 
          PlotPoints -> {2, 2}, 
          DisplayFunction -> Identity];
\(\(scaledAcolumns = {Arrow[A . {h, 0}, Tail -> {0, 0, 0}, 
            VectorColor -> Red], 
          Arrow[A . {0, k}, Tail -> {0, 0, 0}, 
            VectorColor -> Red]};\)
\)
Show[planepiece, scaledAcolumns, ThreeAxes[10], 
        PlotRange -> All, Axes -> True, 
        AxesLabel -> {"x", "y", "z"}, Boxed -> False, 
        ViewPoint -> CMView, 
        DisplayFunction -> $DisplayFunction];

Throw in the plot of the point {1.1, 6.5, 9.1}:

point = {1.1, 6.5, 9.1};
\(\(pointplot = 
        Graphics3D[{PointSize[0.02], Red, Point[point]}];\)
\)
Show[planepiece, scaledAcolumns, pointplot, ThreeAxes[10], 
        PlotRange -> All, Axes -> True, 
        AxesLabel -> {"x", "y", "z"}, Boxed -> False, 
        ViewPoint -> CMView, 
        DisplayFunction -> $DisplayFunction];

The plotted point is not on the plane.
Your job is to come up with the point on the plane closest to the point plotted above.

G.2.c.ii) More plotting

Stay with the same plane as in part i) immediately above and look at this circle which is parallel to the xy-plane:

Clear[x, y, t];
{x[t_], y[t_], z[t_]} = {5 Cos[t] + 4, 5 Sin[t] + 7, 
          9};
\(\({tlow, thigh} = {0, 2\ π};\)
\)
\(\(curveplot = 
        ParametricPlot3D[{x[t], y[t], z[t]}, {t, tlow, thigh}, 
          DisplayFunction -> Identity];\)
\)
Show[planepiece, scaledAcolumns, curveplot, ThreeAxes[10], 
        PlotRange -> All, Axes -> True, 
        AxesLabel -> {"x", "y", "z"}, Boxed -> False, 
        ViewPoint -> CMView, 
        DisplayFunction -> $DisplayFunction];

Your job is to plot the curve on the plane that is composed of the points on the plane closest to corresponding points {x[t],y[t],z[t]} on the circle. 
Do it.

G.3) Phenoms asking to be explained

G.3.a.i)    Xapprox and PseudoInverse

Look at this:

hitdim = Random[Integer, {3, 7}];
hangdim = 
        Random[Integer, {3, 
            7}];

\(\(A\  = \ 
        Table[Random[Real, {\(-4\), 4}], {i, 1, hangdim}, {j, 1, 
            hitdim}];\)
\)
MatrixForm[A]

Look at this:

Y = Table[Random[Real, {-4, 4}], {i, 1, hangdim}]

Make Xapprox

Clear[alignerframe, stretch, hangerframe, 
        k];
hangerframe[
        k_] := SingularValues[
            A][[1]][[k]]
stretch[
        k_] := SingularValues[
            A][[2]][[k]]
alignerframe[
            k_] := SingularValues[
                A][[3]][[k]];

rank = 
          Length[SingularValues[
                A][[2]]];

Xapprox=Expand[(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]]

PseudoInverse[A] . Y

Rerun all cells a couple of times and then describe what you see and explain why you see it.

G.3.a.ii)    Ytest and A.(PseudoInverse[A].Y)

Look at this matrix:

hitdim = Random[Integer, {3, 7}];
hangdim = 
        Random[Integer, {3, 
            7}];

\(\(A\  = \ 
        Table[Random[Real, {\(-4\), 4}], {i, 1, hangdim}, {j, 1, 
            hitdim}];\)
\)
MatrixForm[A]

Look at this:

Y = Table[Random[Real, {-4, 4}], {i, 1, hangdim}]

Make Ytest:

rank = Length[SingularValues[A][[2]]];
Clear[
        hangerframe];
hangerframe[k_] := SingularValues[A][[1]][[k]];
 Ytest=(∑_(k = 1)^rank)(Y.hangerframe[k]) hangerframe[k]

Compare:

A.(PseudoInverse[A].Y)

Rerun all cells a couple of times and then describe what you see and explain why you see it.

G.3.a.iii)    A.X ={0,0,. . .,0} has a guaranteed solution

All clued-in matrix folks know that when you go with any coefficient matrix A , the linear system
          A.X ={0,0,. . .,0} has a solution.
Write down one.

G.3.a.iv)   When you have a linear system with fewer linear equations than unknowns

All clued-in matrix folks know that when you have a linear system with fewer linear equations than unknown x's, then that linear system has a solution other than {0,0,....0}. 
Explain how they know this

G.3.a.v)   Calculus Cal on when you have a linear system with the same number of equations  as unknowns

Just as you are getting ready to leave your computer and go out and relax , you notice a  repulsive figure enter the computer lab.  He sits down next to you and your worst fears are confirmed.  He is the infamous pest Calculus Cal.  Cal is a regular nose-picker; he has foul breath and a strong hint of a garlic/ammonia body odor goes with him at all times.  

To make matters worse, Cal thinks he knows everything and is very vocal about it.  But his knowledge is almost always lacking an important consideration. 

You try to sneak out but Cal catches your eye and says to you through a waft of foul breath.

Cal saw you doing part ii) immediately above. He says to you "When you have a linear system with the same number of equations as unknowns, you can count on an exact solution."
What do you say to Cal?  And then tell him where to go.

G.4) Build an inverse, build a pseudo inverse

G.4.a.i) Build an inverse

Here's a square matrix A:

hitdim  =  Random[Integer, {4, 6}];
hangdim  =  hitdim 
A  = 
        Table[Random[Real, {-2, 2}], {i, 1, hangdim}, {j, 1, 
            hitdim}];
MatrixForm[A]

Here is the full SVD of this matrix:

rank=Length[SingularValues[A][[2]]];
Clear[alignerframe, k];
alignerframe[k_]:=SingularValues[A][[3,k]];
aligner = Table[alignerframe[k], {k, 1, rank}];
MatrixForm[
      aligner]

Clear[stretch];
stretch[k_]:=SingularValues[A][[2,k]];
stretcher = 
          DiagonalMatrix[Table[stretch[k], {k, 1, rank}]];
MatrixForm[stretcher]

Clear[hangerframe, k];
hangerframe[k_]:=SingularValues[A][[1,k]];
hanger = Transpose[Table[hangerframe[k], {k, 1, rank}]];
MatrixForm[
      hanger]

The inverse of A is:

MatrixForm[Inverse[A]]

Use the hanger matrix, the stretch factors , and the aligner calculated above to duplicate this calculation of the  inverse of A

G.4.b.i) Build a transpose

Here is a matrix A:

hitdim  =  Random[Integer, {4, 7}];
hangdim  =  
        hitdim  +  (-1)^Random[Integer, {0, 1}] Random[
              Integer, {1, 3}];
A  = 
        Table[Random[Real, {-2, 2}], {i, 1, hangdim}, {j, 1, 
            hitdim}];
MatrixForm[A]

Here is the full SVD of this matrix:

rank=Length[SingularValues[A][[2]]]

Clear[alignerframe, k];
alignerframe[k_]:=SingularValues[A][[3,k]];
aligner = Table[alignerframe[k], {k, 1, rank}];
MatrixForm[
      aligner]

Clear[stretch];
stretch[k_]:=SingularValues[A][[2,k]];
stretcher = 
          DiagonalMatrix[Table[stretch[k], {k, 1, rank}]];
MatrixForm[stretcher]

Clear[hangerframe, k];
hangerframe[k_]:=SingularValues[A][[1,k]];
hanger = Transpose[Table[hangerframe[k], {k, 1, rank}]];
MatrixForm[
      hanger]

The transpose of A is:

MatrixForm[Transpose[A]]

Use the hanger matrix, the stretch factors , and the aligner calculated above to duplicate this calculation of the transpose of A.

G.4.b.ii)  Rank of (A^t) = rank of  A

Go back and review what you did in part i). 
Then say why it is very vivid that  
            rank of (A^t) = rank of  A.
 for this matrix A or any other matrix A.

G.4.c.i) Build a pseudoinverse

Here is a matrix A:

hitdim  =  Random[Integer, {4, 5}];
hangdim  =  
        hitdim  +  (-1)^Random[Integer, {0, 1}] Random[
              Integer, {1, 2}];
A  = 
        Table[Random[Real, {-2, 2}], {i, 1, hangdim}, {j, 1, 
            hitdim}];
MatrixForm[A]

Here is the full SVD of this matrix:

rank=Length[SingularValues[A][[2]]]

Clear[alignerframe, k];
alignerframe[k_]:=SingularValues[A][[3,k]];
aligner = Table[alignerframe[k], {k, 1, rank}];
MatrixForm[
      aligner]

Clear[stretch];
stretch[k_]:=SingularValues[A][[2,k]];
stretcher = 
          DiagonalMatrix[Table[stretch[k], {k, 1, rank}]];
MatrixForm[stretcher]

Clear[hangerframe, k];
hangerframe[k_]:=SingularValues[A][[1,k]];
hanger = Transpose[Table[hangerframe[k], {k, 1, rank}]];
MatrixForm[
      hanger]

The pseudoinverse of A is:

MatrixForm[PseudoInverse[A]]

Use the hanger matrix, the stretch factors, and the aligner calculated above to duplicate this calculation of the pseudoinverse of A.

G.4.c.ii) Rank of PseudoInverse[A] = rank of  A

Here's a new random matrix:

hitdim  =  Random[Integer, {4, 5}];
hangdim  =  
        hitdim  +  Random[Integer, {1, 2}];
A  = 
        Table[Random[Real, {-2, 2}], {i, 1, hangdim}, {j, 1, 
            hitdim}];
MatrixForm[A]

The pseudo inverse of A is:

MatrixForm[PseudoInverse[A]]

The rank of A is:

rankA=Length[SingularValues[A][[2]]]

The rank of PseudoInverse[A] is:

rankpseudoinverseA=Length[SingularValues[A][[2]]]

Go back and review what you did in part i). 
Then say why it is very vivid that  
            rank of PseudoInverse[A] = rank of  A 
 for this matrix A and all other matrices A.

G.5) Agree or disagree

Here are the basic facts:

◼ The rank of a matrix is the number of non-zero SVD stretch factors for A. 

   This is also the number of hangerframe vectors for A.

◼ A given matrix A is of full rank if the rank of A = hitdim.

◼ Given a linear system A.X = Y,  put 
          Ytest = (∑_(k = 1)^rank)(Y . 
                  hangerframe[k]) hangerframe[k].

◼ If  Ytest = Y, and A is of full rank, then there is exactly one solution of A . X = Y and that solution is
          Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]       (exactly determined).

◼ If  Ytest = Y, and A is of not of full rank, then 
          Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k] 
 is just one of infinitely many solutions  A . X = Y (under determined).

◼ If  Ytest = Y, then
          Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k] 
 is the solution of  A . X = Y of minimum norm.

◼ If Ytest ≠ Y, then the linear system A . X = Y has no solutions at all (over determined).

◼ If A is invertible, then
            (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]=(A^(-1)) . Y.

◼ For any matrix A (invertible or not)
            (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]=PseudoInverse[A] . Y.

◼ For any matrix A (invertible or not)
            Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]  
is as close to  Y as any other possible hit with A.

◼  For any matrix A (invertible or not)
A. (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]   = Ytest 
is as close to  Y as any other possible hit with A.

Indicate whether you agree or disagree with each of the following statements. 
You are not asked to explain your responses, but feel free to throw in an explanation if you want to.

G.5.a.i) Using  Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]

You are given a matrix A and use it to make a linear system 
A . X = Y.  You check the rank of A and learn A is of full rank. Then you do your SVD analysis of A and calculate               
             Ytest = (∑_(k = 1)^rank)(Y . 
                  hangerframe[k]) hangerframe[k]
and then you hit learn that
              Ytest = Y.
You now have enough information to announce that there is exactly one solution of the linear system A.X = Y.

G.5.a.ii) Using  Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]

You are given a matrix A and use it to make a linear system 
A . X = Y.  You check the rank of A and learn A is not of full rank. Then you do your SVD analysis of A and calculate               
             Ytest = (∑_(k = 1)^rank)(Y . 
                  hangerframe[k]) hangerframe[k]
and then you hit learn that
              Ytest = Y.
You now have enough information to announce that there is exactly one solution of the linear system A.X = Y.

G.5.a.iii) Using Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]

You are given a matrix A and use it to make a linear system 
A . X = Y. . Then you do your SVD analysis of A and calculate               
            Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]
and learn that
             Ytest ≠ Y.
You now have enough information to announce that there are no solutions of  the linear system A.X = Y.

G.5.d.vi) Using Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]

You are a given square A (so that hitdim = hangdim) and use it to make a linear system A.X = Y. 
 You do your SVD analysis of A and calculate               
       Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k].
And then you find that
        Ytest  ≠  Y.
This tells you that
        Det[A] = 0.

G.5.b.i) Using rank

You are given a matrix A and use it to make a linear system A . X = Y.  You check the rank of A and learn A is of full rank. 
You now have enough information to announce that there is exactly one solution of the linear system A.X = Y.

G.5.b.ii) Xpprox and Ytest

You are given a matrix A and use it to make a linear system 
A . X = Y.  You do your SVD analysis of A and calculate  
            Xapprox = (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

and

Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]
At this stage, you are guaranteed that 
             A.Xapprox = Ytest.

G.5.b.iii) Using  Xapprox = (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

You are given a matrix A and use it to make a linear system 
A . X = Y.  You check the rank of A and learn A is not of full rank. Then you do your SVD analysis of A and calculate               
             \(Xapprox = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\)
and then you hit Xapprox with A and learn that
            A.Xapprox = Y.
You now have enough information to announce that the linear system A.X = Y has lots of solutions.

G.5.b.iv) Using  Xapprox = (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

You are given a matrix A and use it to make a linear system 
A . X = Y.  You check the rank of A and learn A is not of full rank. Then you do your SVD analysis of A and calculate               
             \(Xapprox = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\)
and then you hit Xapprox with A and learn that
             A.Xapprox ≠ Y.
You now have enough information to announce that the linear system A.X = Y no solution at all.

G.5.c.i) Full rank

You are given a linear system 
         A.X = Y.
You check the rank of A and find that A is of full rank. 
This guarantees that the given linear system has exactly one solution (and no more).

G.5.c.ii) Rank deficient

You are given a linear system 
         A.X = Y.
You check the rank of A and find that A is not of full rank. 
This guarantees that the given linear system has no solution at all.

G.5.c.iii) Rows and rank

Saying that a matrix A is of full rank is the same as saying that the rank of A is equal to the number of horizontal rows of A.

G.5.c.iv) Columns and rank

Saying that a matrix A is of full rank is the same as saying that the rank of A is equal to the number of vertical columns of A.

G.5.c.v)  hitdim > hangdim

You are given a linear system 
         A.X = Y.
with hitdim > hangdim
It is automatic that the coefficient matrix A is not of full rank.

G.5.c.vi)  Full rank

You are given a linear system 
         A.X = Y.
You check the rank of A and find that A is of full rank. This tells you that the system has exactly one solution or no solutions at all.

G.5.c.vii)  Same number of equations as variable x[i]'s

You are given a linear system 
         A.X = Y
with the same number of equations as variable x[i]'s
You check the rank of A and find that A is of full rank. This tells you that the system has exactly one solution.

G.5.c.viii) Rank deficient

You are given a linear system 
            A.X = Y. 
You check the rank of A and find that A is not of full rank. 
This tells you that the system either has no solution or has many solutions.

G.5.c.ix) Rank and hangdim

Given a matrix hitting on hitdimD and hanging in hangdimD, you can be sure that 
that the rank of A cannot be bigger than hangdim.

G.5.c.x) Rank and hitdim

Given a matrix hitting on hitdimD and hanging in hangdimD, you can be sure that 
that the rank of A cannot be bigger than hitdim.

G.5.d.i) Many solutions?

You are a given matrix A that hits on 8D and hangs in 5D.   If the rank of A is 5, then no matter what Y in 5D you go with, you are guaranteed that the corresponding linear system A.X = Y has many solutions.

G.5.d.ii) Rank of A is the same as the rank of (A^t)?

Given any matrix A, you are guaranteed that the rank of A is the same as the rank of (A^t).

G.5.d.iii)  Det[A] ≠ 0

Given any square matrix A (so that hitdim = hangdim), saying that A is invertible is the same as saying that Det[A] ≠ 0.

G.5.d.iv) Rank > hangdim?

It is possible to encounter a matrix A that hangs its hits in 4D and whose rank is 6.

G.5.d.v) Stretches

It is possible to encounter a matrix that hits on 6D and hangs its hits in 4D and has five non-zero SVD stretch factors.

G.5.d.vi) Hanger frame

It is possible to encounter a matrix that hangs its hits in 4D and whose SVD hanger frame consists of five mutually perpendicular unit vectors.

G.5.d.v) Stretches

It is possible to encounter a matrix that hits on 6D and hangs its hits in 4D and has exactly four non-zero SVD stretch factors.

G.5.d.vi) Hanger frame

It is possible to encounter a matrix that hangs its hits in 7D whose SVD hanger frame consists of six mutually perpendicular unit vectors.

G.5.d.vii) hitdim < hangdim

You are given a linear system 
         A.X = Y.
with hitdim < hangdim
It is automatic that the coefficient matrix A is of full rank.

G.5.d.viii) Hitting on 5D and hanging in 8D

You are a given matrix A that hits on 5D and hangs in 8D.   Regardless of the rank of A, there are guaranteed to be many Y's in 8D for which the linear system A.X = Y has no solution.

G.5.e.i) Stretches and determinants

Here is a 3D matrix shown with its nonzero SVD stretch factors:

A = (({{0.5, 0, 1.1}, {1., -0.2, 1.5}, {3.5, -0.8, 4.9}}));
SingularValues[A][[2]]

This is enough to tell you that the rank of A is 2 and that the determinant of A is 0

G.5.e.ii)  Stretches and determinants

Here is a 4D matrix shown with its SVD stretch factors:

A = (({{1., -0.4, 0.8, -1. }, {2., -0.4, 3., 0}, {-0.5, 0, -1.1, -0.5}, {3.5, -0.8, 4.9, -0.5}}));
SingularValues[A][[2]]

This is enough to tell you that the rank of A is 2 and that the determinant of A is:

7.30167   1.22704

G.5.e.iii) Four linear equations in four variables

Here is a linear system  consisting of four equations involving four variables x[1],x[2],x[3] and x[4]:

(Subscript[a, 1]) x[1] + (Subscript[b, 1]) x[2] + (Subscript[c, 1]) x[3] + (Subscript[d, 1]) x[4]  =  y[1]
(Subscript[a, 2]) x[1] + (Subscript[b, 2]) x[2] + (Subscript[c, 2]) x[3] + (Subscript[d, 2]) x[4]  =  y[2]
(Subscript[a, 3]) x[1] + (Subscript[b, 3]) x[2] + (Subscript[c, 3]) x[3] + (Subscript[d, 3]) x[4]  =  y[3]
(Subscript[a, 4]) x[1] + (Subscript[b, 4]) x[2] + (Subscript[c, 4]) x[3] + (Subscript[d, 4]) x[4]  =  y[4]

Given specific values for all the constants

(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3]), (Subscript[a, 4]),  (Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4]),(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3]),(Subscript[c, 4]),(Subscript[d, 1]),(Subscript[d, 2]),(Subscript[d, 3]),(Subscript[d, 4]) and y[1], y[2], y[3], y[4],

you are guaranteed that the resulting linear system has exactly one solution for each of
        x[1], x[2], x[3], x[4]
because you have four equations in four unknowns.

G.5.e.iv) Det[A] ≠ 0

You are given a square matrix A (so that hitdim = hangdim) A and use it to make a linear system A.X = Y.  You calculate Det[A] and find that Det[A] ≠ 0. You know have enough information to announce that the linear system has exactly one solution.

G.5.f.i) Approximate solution

You are given a matrix A and use it to make a linear system 
A . X = Y.  You do your SVD analysis of A and calculate               
             \(Xapprox = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\)
and then you announce that Xapprox the best approximate solution you can get.

G.5.f.ii)  PseudoInverse

You are given a matrix A that hits on 10D and hangs in 13D.   Take a Y in 13D , do your SVD analysis of A and you  calculate               
             \(∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\)
and then you announce that  this is the same as
             PseudoInverse[A] . Y.

G.5.f.iii) PseudoInverse

You are given an invertible matrix A that hits on 9D and hangs in 9D.   You take a Y in 9D , do your SVD analysis of A and you  calculate               
             \(∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\).
And then you announce that  this is the same as
             PseudoInverse[A] . Y.

G.5.f.iv) Inverse and PseudoInverse

If A is an invertible matrix, then
             (A^(-1))  =  PseudoInverse[A].

G.5.f.v) PseudoInverse

A square matrix has in inverse if its determinant is not 0, but all other matrices have a pseudo inverse.

G.5.f.vi) Ytest and invertibility

You are a given square A (so that hitdim = hangdim) and use it to make a linear system A.X = Y. 
 You do your SVD analysis of A and calculate               
        Ytest =   (∑_(k = 1)^rank)(Y . 
                hangerframe[k]) hangerframe[k].
And then you find that
        Ytest  ≠  Y.
This tells you that A is not invertible.

G.5.f.vii) Solution of smallest norm

You are given a matrix A and use it to make a linear system 
A . X = Y.  You check the rank of A and learn A is not of full rank. Then you do your SVD analysis of A and calculate               
             \(Xapprox = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\)
and then you hit Xapprox with A and learn that
             A.Xapprox= Y.
You now have enough information to announce that the linear system A.X = Y is underdetermined (lots of solutions) and to announce that if Xother is any solution, then
       (Sqrt[Xtest . Xtest]) < (Sqrt[Xother . Xother])

G.6) Measurement errors lead to the need for approximate solutions

G.6.a) Measurement errors lead to the need for approximate solutions

You have a summer job with C&M Engineering Co in Latrobe, Pennsylvania, home of Rolling Rock.
Today you and your gang are out in the field measuring the heights of four nearby mountains in the Allegheny range. 
The gang sets up the transit instrument and goes to work. 
You say: "I'll call x[k] the height of the kth mountain."  
The gang makes the following measurements (in feet) from afar:
        x[1] = 1268.1,
        x[2] = 1317.4,
        x[3] = 1723.7,
        x[4] = 1003.1. 
Just to check these measurements, the gang claims to the top of mountain 4 and measures the how much higher than mountain 4 the other mountains are. This leads to the following new measurements:
        x[1] - x[4] = 264.1,
        x[2] - x[4] = 315.0,
        x[3] - x[4] = 719.8.
You put these equations in matrix form:

A=(({{1.0, 0, 0, 0}, {0, 1.0, 0, 0}, {0, 0, 1.0, 0}, {0, 0, 0, 1.0}, {1.0, 0, 0, -1.0}, {0, 1.0, 0, -1.0}, {0, 0, 1.0, -1.0}}));
MatrixForm[
      A]

\(
\(Y = {1268.1, 1317.4, 1723.7, 1003.1, 264.1, 315.0, 
          719.8};\)
\)
hitdim = 4; hangdim = 7;
Clear[x, k];
\(\(X = Table[x[k], {k, 1, hitdim}];\)
\)
ColumnForm[Thread[A . X == Y]]

You can see that this system is has no solution because the fifth equation says
       x[1] - x[4] = 264.1, 
But the first equation combined with the fourth equation says
        x[1] - x[4] is:

1268.1 - 1003.1

This tells you that there is no hope using your measurements to find exact solutions for the four heights 
      {x[1], x[2], x[3], x[4]}. 
But you have no way of knowing which of the measurements are most accurate or which are least accurate. 
So you decide to use all the measurements and you look at the linear system again:

linearsystem

And you say I'll find the x[1], x[2], x[3] and x[4] that make the vector defined by the left sides as close as it can be to the vector defined by the right hand sides.
Do it.

G.6.a.ii) On a different day

On a different day, you and your gang are out in the field measuring the heights of six nearby mountains in the Allegheny range.  Again the gang sets up the transit instrument and goes to work. You say: "I'll call x[k] the height of the kth mountain."  
The gang makes the following measurements (in feet) from afar:
        x[1] = 968.3
        x[2] = 1110.7
        x[3] = 1346.8
        x[4] = 1452.1
        x[5] = 1623.2
        x[6] = 1856.8
Then the gang claims to the top of mountain 1 and measures the how much higher than mountain 1 the other mountains are. This leads to the following new measurements:
        x[2] - x[1] = 145.7
        x[3] - x[1] = 380.1
        x[4] - x[1] = 483.9
        x[5] - x[1] = 658.4
       x[6] - x[1] = 882.5
Using all the measurements made above, give your best shot at the values for 
        x[1], x[2], x[3], x[4], x[5] and x[6] 
that you are going to report. 
Give a short justification for your answer.

G.7) Linear systems  in least square line fit.
 Linear systems in Fourier fit with sine and cosine waves: Trig polynomials.
 Linear systems in Polynomial fit: Interpolating polynomials

G.7..a) Linear systems  in least square line fit

Here is a plot of some points:

\(\(points = {{\(-0.069\), \(-0.148\)}, {0.410, 0.142}, {0.733, 
            0.2377}, 
	\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 	{1.186, 
            0.584}, {1.618, 0.649}, {1.949, 0.969}};\)
\)
pointplot = 
        ListPlot[points, PlotStyle -> {PointSize[0.02], GreenishUmber}, 
          AxesOrigin -> {0, 0}, 
          AxesLabel -> {"t", "y"}];

As you can see, you can't put a straight line through these points, but you can do the next best thing:

Clear[fitter, t];
fitter[t_] = Fit[points, {1, t}, t]

\(\(fitplot = 
        Plot[fitter[t], {t, \(-0.2\), 2.2}, 
          PlotStyle -> {{Thickness[0.01], Black}}, 
          DisplayFunction -> Identity];\)
\)
Show[pointplot, fitplot, 
        DisplayFunction -> $DisplayFunction];

That's nice. 
If you want to know how to come up with the formula for that line, matrix methods can tell you a lot.
Go with this cleared line function:

Clear[x, t, line];
line[t_] = x[1]  +  t  x[2]

Look at the points again:

points

To try to fit a line through all these points you take the first slot of each point, substitute it into the cleared line function and set the result equal to the second slot of the same point:

ColumnForm[
      Table[line[points[[k, 1]]]  ==  points[[k, 2]], {k, 1, 
          Length[points]}]]

This linear system has no solution. 
Why? Because the points do not line up in a perfect straight line. 

So you do the next best thing: You come up with the x[1] and x[2] that make the 6D vector defined by the left hand side as close to the 6D vector defined by the right hand side as it can be.
Use matrix methods to do this and compare your resulting x[1] and x[2]  with:

fitter[t]

G.7.b.i)  Linear systems in Fourier fit with sine and cosine waves

Trig polynomials are combinations of Sine and Cosine waves.
Here are some samples of trig polynomials and their plots:

\(\(high = Random[Integer, {1, 5}];\)
\)
Clear[trigpoly];
\(trigpoly[t_] = 
      Random[Real, {\(-1\), 
            1}] + ∑\+\(k = 1\)\%high Random[Real, {\(-1\), 1}]\ Cos[
              k\ π\ t] + ∑\+\(k = 1\)\%high\(\( Random[
              Real, {\(-1\), 1}]\)\(\ \)\(Sin[
              k\ π\ t]\)\(
\)\)\)
Plot[trigpoly[t], {t, 0, 4}, 
        PlotStyle -> {{Navy, Thickness[0.01]}}, 
        AspectRatio -> (1/GoldenRatio), 
        AxesLabel -> {"t", "trigpoly[t]"}];

Historical Marker:
At the December, 1807 meeting of the French Academy of Sciences, mathematician and engineer Joseph Fourier announced that any periodic function can be approximated by  sums of sine and cosine waves (trig polynomials). Those in attendance thought he was slightly off his rocker and laughed at him, but eventually they were all forced to change their minds because time has confirmed that Fourier's idea was one of the great scientific break-throughs of all time.  In this problem, you get the opportunity to get your hands on Fourier's idea in action.

Here's how it goes.
You start with any periodic function f[t] and identify its period L:

Clear[f, t];
f[t_] = 
        t - 2 Floor[(t/2)] - 4 Abs[0.5 t - Round[0.5 t]];
L = 2;
\(\(cycles = 3;\)
\)
Plot[f[t], {t, 0, cycles L}, 
        PlotStyle -> {{Thickness[0.02], Navy}}, 
        AxesLabel -> {"t", "f[t]"}, 
        PlotLabel -> "cycles" cycles, 
        Epilog -> {{Red, Thickness[0.02], 
              Line[{{0, 0}, {L, 0}}]}, {Text["One Period", {(L/2), 
                  0.1}]}}, AspectRatio -> Automatic];

Plot one period of f[t] and throw in an odd number of equally spaced points on the plot of f[t]:

periodplot = 
          Plot[f[t], {t, 0, L}, PlotStyle -> {{Thickness[0.02], Navy}}, 
            AxesLabel -> {"t", "f[t]"}, 
            DisplayFunction -> Identity];

jump = (L/8);
fpoints = Table[N[{t, f[t]}], {t, jump, L - jump, jump}];

pointplot = 
          ListPlot[fpoints, PlotStyle -> {Goldenrod, PointSize[0.06]}, 
            DisplayFunction -> Identity];

setup = 
        Show[periodplot, pointplot, 
          DisplayFunction -> $DisplayFunction];

Define a trig polynomial whose constituent Cosine and Sine waves are fixed to match the period of f[t] and has one undetermined coefficient x[k] for each point

b = (1/2) (Length[fpoints] - 1);

Clear[trigpoly, x, t, k];
trigpoly[t_]=x[1]+(∑_(k = 2)^(b + 1))x[k] Cos[(((k - 1) (2 π) t)/L)]+(∑_(k =  b + 2)^(Length[fpoints]))x[k] Sin[(((k - b - 1) (2 π) t)/L)]

To get the approximation, you go after the x[k]'s that make the plot trigpoly[t] go right through the plotted points. This results in the following linear system:

\(linearsystem\  = \ 

      ColumnForm[
        Table[trigpoly[fpoints[\([k, 1]\)]]\  == \ fpoints[\([k, 2]\)], 
{k, 
            1, Length[fpoints]}]]\)

Pluck off the coefficient matrix A:

\(\(linearformulas = 
        Table[trigpoly[
            fpoints[[k, 1]]], {k, 1, 
            Length[fpoints]}];\)
\)
\(\(variables = Table[x[k], {k, 1, Length[fpoints]}];\)
\)
A = Outer[D, {linearformulas}, 
            variables][[1]];
MatrixForm[A]

Check the rank of A:

rank = 
      Length[SingularValues[
            A][[2]]]

Good. 
A is of full rank. This tells you that Xsol is the only solution.
Get the x[k]'s:

Y = 
          Table[fpoints[[k, 2]], {k, 1, 
              Length[fpoints]}];

Clear[alignerframe, stretch, hangerframe, k];
alignerframe[k_]:=SingularValues[A][[3,k]];

stretch[k_]:=SingularValues[A][[2,k]];

hangerframe[k_]:=SingularValues[A][[1,k]];

 Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

The trig polynomial you are after is:

Clear[trigfitter];
trigfitter[t_]=trigpoly[t]/.Table[x[k] -> 
                Xsol[[k]],{k, 0, Length[fpoints]}]

See it:

\(\(fitplot = 
        Plot[trigfitter[t], {t, 0, L}, 
          PlotStyle -> {{CadmiumOrange, Thickness[0.015]}}, 
          DisplayFunction -> Identity];\)
\)
Show[setup, fitplot, 
        DisplayFunction -> $DisplayFunction];

Nailed all the points.
See it over some extra cycles:

cycles = 3;
Plot[{f[t], trigfitter[t]}, {t, 0, cycles L}, 
        PlotStyle -> {{Thickness[0.02], Navy}, {Thickness[0.015], 
              CadmiumOrange}}, AxesLabel -> {"t", ""}, 
        PlotLabel -> "cycles" cycles, 
        Epilog -> {{Red, Thickness[0.02], 
              Line[{{0, 0}, {L, 0}}]}, {Text["One Period", {(L/2), 
                  0.1}]}}, AspectRatio -> Automatic];

But the approximation is not all that great.

Copy paste and edit reducing jump from (L/8) to (L/16) and use the resulting plot to assess how well the resulting new trigfitter[t] approximates the original function f[t].

G.7.b.ii) More Fourier fit

Go with this square wave function f[t]:

Clear[f, t];
f[t_] = Sign[Sin[π t]]  +  1;
L = 2;
\(\(cycles = 4;\)
\)
Plot[f[t], {t, 0, cycles L}, 
        PlotStyle -> {{Thickness[0.02], Navy}}, 
        AxesLabel -> {"t", "f[t]"}, 
        PlotLabel -> "cycles" cycles, 
        Epilog -> {{Red, Thickness[0.02], 
              Line[{{0, 0}, {L, 0}}]}, {Text["One Period", {(L/2), 
                  0.1}]}}];

Come up with a trig polynomial that approximates this function fairly well.
Show off your work with winning plots.

G.7.c.i) How you make an interpolating polynomial

Here are some points:

tlow = Random[Real, {-1, 1}];
\(\(thigh = Random[Real, {3, 4}];\)
\)
jump = ((thigh - tlow)/6);
\(\(points = 
        Table[{t, Random[Real, {0, 2}]}, {t, tlow, thigh, 
            jump}];\)
\)
pointplot = 
        ListPlot[points, PlotStyle -> {Coral, PointSize[0.03]}, 
          AxesLabel -> {"t", "y"}];

There is an old adage that says there is a unique polynomial of least degree whose plot runs through all these plotted points.
Mathematica can calculate this polynomial it is:

Clear[leastdegree, t];
leastdegree[t_] = 
      Expand[InterpolatingPolynomial[points, t]]

Here are the same points shown together with a plot of this function.

\(\(curveplot = 
        Plot[leastdegree[t], {t, tlow, thigh}, 
          PlotStyle -> {{Thickness[0.01], Indigo}}, 
          DisplayFunction -> Identity];\)
\)
Show[pointplot, curveplot, pointplot, PlotRange -> All, 
        DisplayFunction -> $DisplayFunction];

The polynomial curve runs right through the points.
One way of seeing where this interpolating polynomial comes from involves linear systems.
First you look at the points:

points

Next you make a cleared polynomial function of t with one undetermined coefficient x[k] for each given point {t[k],y[k]}:

Clear[clearpoly, x, k];
clearpoly[t_]=(∑_(k = 1)^(Length[points]))x[k] (t^(k - 1))

To nail down the 
        x[1], x[2], x[3], x[4], x[6], x[7] 
that make this into the interpolating polynomial running through all the points, you take the first slot of each point, substitute it into the cleared polynomial function and set the result equal to the second slot of the same point:

\(linearsystem\  = \ 

      ColumnForm[
        Table[clearpoly[points[\([k, 1]\)]]\  == \ points[\([k, 2]\)], {k, 1, 
            Length[points]}]]\)

You can see that this is a linear system with the same number of variable x[k]'s as equations.
To save a lot of typing, you can get the coefficient matrix A this way:

\(\(linearformulas = 
        Table[clearpoly[
            points[[k, 1]]], {k, 1, 
            Length[points]}];\)
\)
\(\(xvariables = Table[x[k], {k, 1, Length[points]}];\)
\)
A = Outer[D, {linearformulas}, 
            xvariables][[1]];
MatrixForm[A]

Now look at:

Det[A]

How does this signal that the linear system has exactly one solution?
Put answer here.

Use matrix methods to solve the linear system in part i) for 
           x[1],x[2],x[3],x[4],x[5],x[6] and x[7], 
substitute your x-values into the cleared polynomial and compare the resulting polynomial to:

Expand[InterpolatingPolynomial[points, t]]

Describe what you get.

G.7.c.i) Theory of interpolating polynomials

To make an interpolating polynomial, you start with given points in the form {t,y}  all with different t- slots.

Clear[t, y, i];
points = Table[{t[i], y[i]}, {i, 1, 6}]

Next you make a cleared polynomial function of t with one undetermined coefficient x[i] for each given point {t[i],y[i]}:

Clear[clearpoly];
clearpoly[t_]=(∑_(k = 1)^(Length[points]))x[k] (t^(k - 1))

To nail down the 
        x[1], x[2], x[3], x[4], x[6] 
that make this into the interpolating polynomial running through all the points, you take the first slot of each point, substitute it into the cleared polynomial function and set the result equal to the second slot of the same point:

\(linearsystem\  = \ 

      ColumnForm[
        Table[clearpoly[t[k]]\  == \ y[k], {k, 1, 
            Length[points]}]]\)

You can see that this is a linear system with the same number of variable x[k]'s as equations.
To save a lot of typing, you can get the coefficient matrix A this way:

linearformulas = 
        Table[clearpoly[t[k]], {k, 1, Length[points]}];
\(\(xvariables = Table[x[k], {k, 1, Length[points]}];\)
\)
A = Outer[D, {linearformulas}, 
            xvariables][[1]];
MatrixForm[A]

That's the right coefficient matrix A and now that you have it, let Mathematica calculate Det[A] in factored form:

Factor[Det[A]]

How does this signal that if the given points all have different first slots (t[i]'s), then you are guaranteed that for each given {y[1],y[2],y[3],y[4],y[5],y[6]} this linear system has exactly one solution for {x[1],x[2],x[3],x[4],x[5],x[6]}?
How  does signal that there is precisely one polynomial of the form
       poly[t] = a[1] + a[2] t + a[3] (t^2)+ a[4] (t^3) + a[5] (t^4)+ a[6] (t^5)
with the property that
       poly[t[k]] = y[k] for k =1,2,3,4,5,6

G.8) Linear systems in calculus and diffeq

G.8.a) The method of undetermined coefficients: Systems of linear equations in integration

When you go with a function f[x] and try to calculate 
           (∫_0^t) f[x] dx, 
you look for a function g[x] with 
          g'[x]=f[x]. 
If you find g[x], then you can announce that
         (∫_0^t) f[x] dx = g[t] - g[0]

In the case 
          f[x] = (x^3) (e^(-x)), 
it's a good idea to look for g[x] like this:

Clear[x, a, k];
g[x_] = (∑_(k = 0)^3) a[
            k + 1] (x^k) (E^(-x))

Look at g'[x] with an eye toward setting the coefficients a[1],a[2],a[3] and a[4] so that 
g'[x] = f[x] = (x^3) (e^(-x)):

g'[x]

Collect the terms and factor:

Collect[g'[x], {1, x, (x^2), (x^3)}]

So g'[x] = -(e^(-x)) a[1] + (e^(-x)) a[2] + 
        x (-(e^(-x)) a[2] + 
              2 (e^(-x)) a[
                  3]) - (e^(-x)) (x^3) a[4] + 
        (x^2) (-(e^(-x)) a[3] + 
              3 (e^(-x)) a[4])

You want this to be f[x] = (x^3) (E^(-x)), so you write
                         -a[1] + a[2] = 0

- a[2] + 2 a[3] = 0

-a[3] + 3 a[4] = 0

-a[4] = 1

Solve these equations by hand and then use matrix methods to confirm your hand calculation. Use your solutions for 
a[1], a[2], a[3], a[4] to calculate
        (∫_0^t) (x^3) (E^(-x)) dx

Solving the equations
                         -a[1] + a[2] = 0

- a[2] + 2 a[3] = 0

-a[3] + 3 a[4] = 0

-a[4] = 1
by hand is really easy.

You use the last equation to nail down 
          a[4] = -1.
 Put   a[4] = -1. into the equation
              -a[3] + 3 a[4] = 0 
to get
              -a[3] - 3 = 0.
This gives a[3] = -1.
Now you take over.

To do this by matrix methods, use an inverse matrix.

G.8.b) Systems of linear equations for coming up with formulas for solutions of linear differential equations

You have the job of coming up with a formula for the solution of the fourth order differential equation
 y''''[t] + 2.9 y'''[t] + 3.04 y''[t] + 1.356 y'[t] + 0.216 y[t] = 0
with 
   y[0] = 2.3, y'[0] = -2.7, y''[0] = -3.7 and y'''[0] = 0.8 
  Remembering the characteristic equation, you type:

Clear[z];
charequation = 
      (z^4) + 2.9 (z^3) + 3.04 (z^2) + 1.356 z + 0.216 == 0

And solve it:

Solve[charequation, z]

Next you put:

Clear[trialy, t, x];
trialy[t_] = 
      x[1] (E^(-t)) + x[2] (E^(-0.9 t)) + x[3] (E^(-0.6 t)) + 
        x[4] (E^(-0.4 t))

Next,your incorporate the starter data
  y[0] = 2.3, y'[0] = -2.7, y''[0] = -3.7 and y'''[0] = 0.8 :

\(ColumnForm[
      Thread[{trialy[0]\  == \ 
            2.3, 
\(trialy'\)[0]\  == \ \(-2.7\), 
\(\(trialy'\)'\)[
              0]\ \  == \ \(-3.7\), 
\(\(\(trialy'\)'\)'\)[0]\  == \ 
            0.8}]]\)

Hey dude!  That's a linear system.
To nail down a formula the solution of  
  y''''[t] + 13.4 y'''[t] + 48.41 y''[t] + 51.484 y'[t] + 8.466 y[t] = 0
with 
   y[0] = 2.3, y'[0] = -2.7, y''[0] = -3.7 and y'''[0] = 0.8, 
   you solve the linear system above for x[1],x[2],x[3] and x[4] and insert the results in trialy[t]. 
Use matrix methods to do this.

G.8.c) Systems of linear  equations for the method of undetermined coefficients in diffeq

Advanced folks know that all solutions of the forced, damped oscillator 
  y''[t] + 0.4 y'[t] + 1.3 y[t] = 
              1.2 Sin[1.2 t] + 0.5 Cos[2.3 t] + 0.6 Sin[4.1 t]
have the same long term behavior regardless of what the given initial data on y[0] and y'[0] are. So to see what the long term behavior any solution y[t], all you need it to come up with any one particular solution of it and plot it.
Glancing at the right side, you get the idea you can make a solution out of:

Clear[undeterminedy, t, x];
\(undeterminedy[t_] = 
	
      x[1]\ Sin[1.2\ t] + x[2]\ Cos[1.2\ t] + 
	x[3]\ Sin[2.3\ t] + 
        x[4]\ Cos[2.3\ t] + 
	x[5]\ Sin[4.1\ t] + 
        x[6]\ Cos[4.1\ t]\)

To try to determine the x[1], x[2], x[3], x[4] ,x[5] and x[6] that give a solution, you replace y[t] by undeterminedy[t] in
                y''[t] + 0.4 y'[t] + 1.3 y[t]:

Clear[y];
left=Expand[y''[t]+0.4 y'[t]+1.3 y[t]/.{y[t] -> undeterminedy[t],y'[t]->undeterminedy'[t],y''[t]->undeterminedy''[t]}]

You want this to be  1.2 Sin[1.2 t] + 0.5 Cos[2.3 t] + 0.6 Sin[4.1 t], so you type:

\(\(right = 
        1.2\ Sin[1.2\ t] + 0.5\ Cos[2.3\ t] + 0.6\ Sin[4.1\ t];\)
\)
\(ColumnForm[
      Thread[{Coefficient[left, Sin[1.2\ t]]\  == \ 
            Coefficient[right, Sin[1.2\ t]], 

          Coefficient[left, Cos[1.2\ t]]\  == \ 
            Coefficient[right, Cos[1.2\ t]], 

          Coefficient[left, Sin[2.3\ t]]\  == \ 
            Coefficient[right, Sin[2.3\ t]], 

          Coefficient[left, Cos[2.3\ t]]\  == \ 
            Coefficient[right, Cos[2.3\ t]], 

          Coefficient[left, Sin[4.1\ t]]\  == \ 
            Coefficient[right, Sin[4.1\ t]], \ 

          Coefficient[left, Cos[4.1\ t]]\  == \ 
            Coefficient[right, Cos[4.1\ t]]}]]\)

That's a 6D linear system.
Solve by hand or use matrix methods to solve it (your choice). 
Then substitute your solutions for x[1],x[2],x[3],x[4],x[5] and x[6] into trialy[t], thereby obtaining a formula for one solution y[t] of
  y''[t] + 0.4 y'[t] + 1.3 y[t] = 
              1.2 Sin[1.2 t] + 0.5 Cos[2.3 t] + 0.6 Sin[4.1 t].
Exhibit the long term behavior of this oscillator by plotting your formula over a healthy range of t's.

G.9) Null space and rank

G.9.a.i) Null frames

Here's a random matrix that hits on 5D and hangs in 3D

hitdim = 6; hangdim = 3;
Clear[a, i, j];
\(\(a[i_, j_] := Random[Real, {\(-2\), 2}];\)
\)
A = Table[a[i, j], {i, 1, hangdim}, {j, 1, hitdim}];
MatrixForm[A]

Clear[nullframe];
nullframe[
          k_] := nullframe[k] = NullSpace[A][[
            k]];
zeroes = Table[nullframe[k], {k, 1, Length[NullSpace[A]]}]

You can think of these as SVD aligner frames for A that correspond to zero stretch factors.

Table[A . nullframe[k], {k, 1, Length[NullSpace[A]]}]

Here are the SVD aligner frames for A that correspond to positive stretch factors:

\(\(rank = 
        Length[\(SingularValues[
              A]\)[[2]]];\)
\)
Clear[alignerframe, k];
alignerframe[k_] := SingularValues[A][[3, 
          k]];
nonzeroes = Table[alignerframe[k], {k, 1, rank}]

Notice that:

Length[zeroes] + Length[nonzeroes] == hitdim

Is this just an isolated accident or will it happen in every set up?
Explain.

G.9.a.ii) No null frames

Here's a random matrix of full rank:

hitdim = 4; hangdim = 6;
Clear[a, i, j];
\(\(a[i_, j_] := Random[Real, {\(-2\), 2}];\)
\)
A = Table[a[i, j], {i, 1, hangdim}, {j, 1, hitdim}];
MatrixForm[A]

rank = 
      Length[SingularValues[
            A][[2]]]

Now look at this:

NullSpace[A]

Nothing there.
Why was this guaranteed to happen?

G.9.b.i) No null frames

Here's a linear system:

\(\(A = {{2.686, \(-0.666\), \(-3.0464\), \(-2.1568\), \(-1.881\)}, {3.068, 2.667, 4.596, \(-3.594\), \(-2.986\)}, {\(-1.918\), 
            2.915, \(-3.093\), \(-0.0829\), \(-0.354\)}, {1.15, 
            5.582, \(-1.497\), \(-6.6769\), \(-3.34\)}, {\(-0.3820\), \(-5.333\), \(-4.6424\), 3.4372, 1.105}};\)
\)
hitdim = 5; hangdim = 5;
Clear[x, k];
X = Table[x[k], {k, 1, hitdim}];
\(\(Y = {\(-2.34\), 1.84, \(-2.98\), 7.04, \(-6.70\)};\)
\)
ColumnForm[ Thread[A . X == Y]]

Intent on solving it quickly, you look at:

Xapprox = PseudoInverse[A] . Y;
A . Xapprox

Y

NullSpace[A]

And announce that this system has exactly one solution. And its solution is:

Xapprox

Are you right?

G.9.b.ii) Some null frames

Here's another linear system:

\(\(A = {{2.686, \(-0.666\), \(-3.0464\), \(-0.1568\), \(-1.881\)}, {3.068, 2.667, 1.596, \(-3.594\), \(-2.986\)}, {\(-1.918\), 
            2.915, \(-3.093\), \(-0.0829\), \(-0.354\)}, {1.15, 
            5.582, \(-1.497\), \(-3.6769\), \(-3.34\)}, {\(-0.3820\), \(-3.333\), \(-4.6424\), 3.4372, 1.105}};\)
\)
hitdim = 5; hangdim = 5;
Clear[x, k];
X = Table[x[k], {k, 1, hitdim}];
\(\(Y = {1.71986, 14.9293, 23.1721, 38.1014, \(-13.2094\)};\)
\)
linearsystem  =  ColumnForm[Thread[A . X  ==   Y]]

Intent on solving it quickly, you look at:

Xapprox = PseudoInverse[A] . Y;
A . Xapprox

Y

NullSpace[A]

And announce that this system is underdetermined (lots of solutions) and the solution of minimum norm is:

Xsol

Are you right?

G.9.b.iii) Some null frames

Here's another linear system:

\(\(A = {{2.686, \(-0.666\), \(-3.0464\), \(-0.1568\), \(-1.881\)}, {3.068, 2.667, 1.596, \(-3.594\), \(-2.986\)}, {\(-1.918\), 
            2.915, \(-3.093\), \(-0.0829\), \(-0.354\)}, {1.15, 
            5.582, \(-1.497\), \(-3.6769\), \(-3.34\)}, {\(-0.3820\), \(-3.333\), \(-4.6424\), 3.4372, 1.105}};\)
\)
hitdim = 5; hangdim = 5;
Clear[x, k];
X = Table[x[k], {k, 1, hitdim}];
\(\(Y = {0.12, 3.14, 2.72, \(-5.06\), \(-2.43\)};\)
\)
linearsystem  =  ColumnForm[Thread[A . X  ==   Y]]

Intent on solving it quickly, you look at:

Xapprox = PseudoInverse[A] . Y;
A . Xapprox

Y

And announce that this system is overdetermined.
 Are you right?

G.10) Columns and rank, rows and rank

G.10.a) Columns and rank

Here is a matrix entered via its vertical columns:

Clear[column, j];
column[1] = {0.96, 1.58, 0.91, -0.22, 0.97};
column[2] = {0.57, 
          1.55, -1.42, -1.42, -1.45};
column[
          3] = {2.07, -2.91, -2.74, -1.40, -2.87};
column[4] = {-0.76, 0.74, -2.56, -0.61, 
          0.78};
\(\(column[5] = {1.36, \(-0.93\), 2.73, 2.38, \(-2.49\)};\)
\)
A = Transpose[Table[column[j], {j, 1, 5}]];
MatrixForm[A]

This matrix is of full rank:

rank = 
      Length[SingularValues[
            A][[2]]]

Your job is to come up with a nonzero 5D vector
        Z = {z[1],z[2],z[3],z[4],z[5]} 
so that when you reset  
        column[3] = Z, 
then the resulting matrix is not of full rank.
Check yourself by calculating the rank of the resulting matrix.

G.10.b) Rows and rank

Here is a matrix entered via its horizontal rows:

Clear[row, i];
row[1] = {1.96, 1.78, 2.91, -3.22};
row[2] = {3.69, 2.55, -0.42, -2.27};
row[3] = {2.17, -3.01, 4.74, -2.30};
\(\(row[4] = {\(-1.76\), 1.74, \(-0.96\), 0.97};\)
\)
A = Table[row[i], {i, 1, 4}];
MatrixForm[A]

This matrix is of full rank:

rank = 
      Length[SingularValues[
            A][[2]]]

Your job is to come up with a nonzero 4D vector
        Z = {z[1],z[2],z[3],z[4]} 
so that when you reset  
        row[2] = Z, 
then the resulting matrix is not of full rank.
Check yourself by calculating the rank of the resulting matrix.

G.11) Inverses and eigenvalues via row reduction

G.11..a.i) Inverses via row reduction

Here's random 4D matrix A

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

Augment A with this right column:

Clear[y];
Y  =  Table[y[i], {i, 1, dim}];
 augmentedA  = AppendRows[A, Transpose[{Y}]];
MatrixForm[augmentedA]

Look at this:

rowechelon  =  
        Expand[RowReduce[augmentedA]];
MatrixForm[rowechelon]

And look at this:

Inverse[A] . Y

Describe what you see and explain why you see it.

G.11.a.ii) Eigenvalues via row reduction

▪ Here's random 2D matrix A :

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

The eigenvalues of A are the λ's for which the matrix 
            Aλ = A -λ IdentityMatrix[dim]
 is not invertible:

Aλ  =  
        A  - λ IdentityMatrix[dim];
MatrixForm[Aλ ]

Augment Aλ with this right column:

Clear[y];
Y  =  Table[y[i], {i, 1, dim}];
 augmentedAλ  = 
        AppendRows[Aλ, Transpose[{Y}]];
MatrixForm[augmentedAλ]

Look at this:

rowechelon  = 
        Together[
          RowReduce[augmentedAλ]];
MatrixForm[rowechelon]

Remembering that the eigenvalues of A are the λ's for which Aλ is not invertible, give a quadratic equation you can solve to come up with the eigenvalues of A. Explain why your equation is correct.
Put answer here.

▪ Compare your equation with:

Det[Aλ]  ==  0

Say why your are not surprised.
Put answer here.

▪ Here's random 4D matrix A :

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

The eigenvalues of A are the λ's for which the matrix 
            Aλ = A -λ IdentityMatrix[dim]
 is not invertible:

Aλ  =  
        A  - λ IdentityMatrix[dim];
MatrixForm[Aλ ]

Augment Aλ with this right column:

Clear[y];
Y  =  Table[y[i], {i, 1, dim}];
 augmentedAλ  = 
        AppendRows[Aλ, Transpose[{Y}]];
MatrixForm[augmentedAλ]

Look at this:

rowechelon  = 
        Together[
          RowReduce[augmentedAλ]];
MatrixForm[rowechelon]

Remembering that the eigenvalues of A are the λ's for which  Aλ is not invertible, give a fourth degree equation you can solve to come up with the eigenvalues of A. Explain why your equation is correct.
Put answer here.

▪ Compare your equation with:

Det[Aλ]  ==  0

Say why your are not surprised.
Put answer here.