Physics 2301 -- Linear Math 2

Meant to be read and evaluated while watching the LinearMath2 video.
Please pause along the way and explore by changing parameters.

Of course Mathematica can rotate in three dimensions:

$Post=MatrixForm
rvec = {x,y,z}
RotationMatrix[θ,{0,0,1}] (* the familiar Chapter 8 version *)
%.rvec
%%.{0,0,1} (* just noting that zhat is an eigenvector with eigenvalue=1 *)
RotationMatrix[θ,{1,1,1}] 
%.{1,1,1}//Simplify (* and now {1,1,1} is an eigenvector, i.e on the spin axis *)

Here we create a square frame of point masses and spin it with constant ω-vector.
Note that unless we choose to spin in a principal direction, L will be changing and something is providing a torque.

(sA={0,0,0};
sB={1,0,0};
sC={1,1,0};
sD={0,1,0};) (* the round parenthesis here suppresses some output *)
sqframe={{Black,Thick,Line[{sA,sB,sC,sD,sA}]},{Red,(Sphere[#,.03]&/@ {sA,sB,sC,sD})}}
Graphics3D[sqframe,PlotRange->Sqrt[2]]

ωvec={0,0,1} (*Change this one and try different directions *)
Animate[Graphics3D[{{Green,Thick,Arrow[{sA,sA+ωvec}]},Rotate[sqframe,θ,ωvec]},PlotRange->Sqrt[2]],{θ,0,2π}]

Or instead let’s draw a triangular frame, animate its rotation, and compute Lvec.
There are two m’s at A and B, and one 2m at C

rA={-2,-1,0}
rB = {2,-1,0}
rC = {0,1,0}
rA+rB+2rC  (* check that the CM is at 0 *)
triframe = {{Thick,Black,Line[{{rA,rB,rC,rA}}]},{Red,Sphere[#,.1]&/@{rA,rB,rC}}}
Graphics3D[triframe]

Now choose an omega and rotate.

ωvec={0,1,0} (* again explore by changing ω *)
Animate[Graphics3D[{{Thick,Green,Arrow[{{0,0,0},ωvec}]},{Rotate[triframe,θ,ωvec]}},PlotRange->3],{θ,0,2π}]

Next let’s add in angular momentum.  Note that unless ω-vec is exactly in a special principal direction,
L-vec will precess, meaning something will have to be providing a torque for this motion to be possible.

ωvec={2,2,0} (* again explore by changing ω *)
LA = Cross[a rA,m Cross[ω ωvec, a rA]] (* here we restore m ω and a just for a moment *)
LB = Cross[a rB,m Cross[ω ωvec, a rB]] (* to make the point that there really are units *)
LC = Cross[a rC,2m Cross[ω ωvec, a rC]]
Lvec = LA+LB+LC /.{m->1,a->1,ω->1/6};  (* but for drawing we choose numbers *)
Animate[
Graphics3D[{{Thick,Green,Arrow[{{0,0,0},ωvec}]},{Rotate[{triframe,{Thick,Blue,Arrow[{{0,0,0},Lvec}]}},θ,ωvec]}},PlotRange->3],{θ,0,2π}]

And now let’s find the full-blown moment of inertia matrix
by adding up the matrices for each point mass.

rvec={{x},{y},{z}} ; (* this is an explicitly 3x1 column vector *)
Transpose[rvec] (* which transposes to a row vector *)
inner=Transpose[rvec].rvec (* this one gives a 1x1 matrix containing the dot product *)
outer=rvec.Transpose[rvec] (* this one gives a 3x3 matrix of all possible products *)
Tr[inner] IdentityMatrix[3] - outer (* and this one is the moment of inertia, with m=1 *)

Let’s make a function which takes in a vector and returns the associated Imatrix,
and use it for each of the point masses:

Imat[{x_,y_,z_}]=Tr[inner] IdentityMatrix[3] - outer
ImatA= Imat[rA]
ImatB= Imat[rB]  (* mB = mA = 1 *)
ImatC= 2Imat[rC] (* mC = 2 *)
Itot = ImatA+ImatB+ImatC  (* and note the principal directions are exactly as expected *)

We didn’t really need the matrix for the symmetric frame,when our heuristics tell us to consider spinning about
(a) symmetry axes, (b) perpendicular to the plane, and (c) perpendicular to the others.
But lastly let’s create a truly random triangular frame.

(qqA=4*{Random[],Random[],Random[]}; (*Random[] give a number between 0 and 1 *)
qqB=4*{Random[],Random[],Random[]} ;
qqC=4*{Random[],Random[],Random[]} ;
qqCM =(qqA+qqB+qqC)/6;  (*find the CM, and define qA, qB and qC relative to it *)
qA=qqA-qqCM;
qB = qqB - qqCM;
qC = qqC - qqCM;)
qframe = {{Thick,Black,Line[{qA,qB,qC,qA}]},{Red,Sphere[#,.1]&/@{qA,qB,qC}}}
Graphics3D[qframe]

Here we form the Imatrix and ask for its Eigenvectors and Eigenvalues.

qmat = Imat[qA]+Imat[qB]+Imat[qC]  (*all m=1*)
eigvecs=Eigenvectors[qmat]
Eigenvalues[qmat]
Eigensystem[qmat] (*this one gives eigenvectors together with eigenvalues *)

And lastly let’s animate as we spin about some fixed ω

ωvec={0,0,1};
zero={0,0,0};
Lvec = qmat.ωvec/2 (*scale by a factor here just for the sake of the picture *)
Animate[Graphics3D[{
{Rotate[{qframe,{Blue,Thick,Arrow[{zero,Lvec}]},
{Black,Arrow[{zero,eigvecs[[1]]}],Arrow[{zero,eigvecs[[2]]}],Arrow[{zero,eigvecs[[3]]}]
}},θ,ωvec]},
{Green,Arrow[{zero,ωvec}]}},PlotRange->3],{θ,0,2π}]