Physics 230 -- Lab 10 (Linear Algebra)

Fall 2015 - Transtrum

Labs developed by John S. Colton, BYU Physics & Astronomy, based in part on past labs by Branton Cambell and Steve Turley

In most fields of science and engineering, the easiest multi-variable problems are the "linear" ones, where the unknown quantities appear in linear combinations with constant coefficients.  Linear problems are usually solved by straightforward methods and have exact solutions.  Non-linear problems, on the other hand, can become very difficult very quickly.  In this lab, we will explore matrix and vector properties and operations, practice formulating problems in terms of matrix/vector equations, and learn how to use Mathematica's linear-algebra capabilities to solve simple and moderately-challenging math and physics problems.

Clear["`*" ]

Review

Defining matrices

As a reminder, when you define matrices, do NOT define them in matrix form. Instead, use MatrixForm after the definition to view them.

incorrect={{1,2,3},{4,5,7},{8,9,10}}   //MatrixForm  (* this will produce bad results *)
incorrect //Det  (* The determinant function fails to work *)

correct={{1,2,3},{4,5,7},{8,9,10}}  
% //MatrixForm  
correct //Det (* this works! *)

Matrix multiplication

Mathematica's matrix multiplication operation is called Dot (shortcut "."), and is the same operator used for vector dot products. 

In order for two matrices to be multipliable, the number of columns in the first matrix must equal the number of rows in the second matrix.  If, for example, A is an m×n matrix and B is an n×p matrix, then A.B will be an m×p matrix.

A= {{1,2},{3,4},{5,6}};  (* a 3x2 matrix *)
% //MatrixForm
B = {{8,7,6,5},{4,3,2,1}};  (* a 2x4 matrix *)
%//MatrixForm
A.B //MatrixForm  (* a 3x4 matrix *)  (* also, it's OK to do MaxtrixForm like this here since I'm not defining a new matrix at the same time. If I wanted C=A.B, then I couldn't use the //MatrixForm bit on the same line *)
B.A // MatrixForm  (* doesn't work; wrong dimensions *)

Map

In a previous lab (the “functions and lists” lab, to be specific), we introduced you to the Map function, shortcut /@. As a reminder, it can be used to apply a function to a list of objects. Here it’s used to apply the Sin function to each number in a list of numbers.

list={.1, .2, .3, .4, .5};
Sin[#]&  /@list  (* "for each element x of list, calculate Sin[x]" *)

The same thing could be accomplished in a seemingly more simple form via Sin[list], but I myself would nearly always use the Map function for this because the translation, “for each element x of list, calculate Sin[x]” is exactly how I myself would think of the situation. That is, I would think: “I have a collection of objects, and I want to apply a function to each one... how to do that? Oh, yes, I should use Map!”

If you can internalize that translation and that way of thinking, you will have a powerful tool at your disposal. It’s closely related to something called a “for each loop” which exists in nearly all programming languages, see http://en.wikipedia.org/wiki/Foreach_loop.

Similarly, you can use Map to apply a function to each matrix in a list of matrices.

A= RandomReal[{0,10},{3,3}] ;(* a random 3x3 matrix with elements from 0-10 *)
B= RandomReal[{0,10},{3,3}];
mlist = {A,B,A.B};  (* a list of 3 matrices *)

MatrixForm[#]&  /@mlist  (* "for each matrix m in mlist, display MatrixForm[m]" *)

Det[#]&  /@ mlist  (* "for each matrix m in mlist, calculate Det[m]" *)
Det /@ mlist (* same thing *)

In this particular case, you can’t use MatrixForm[mlist] or Det[mlist] to get the desired results for a list of matrices like you could with Sin[list] in the previous example, given a list of numbers. (Try it!) So if you want to do a function on a list of matrices and end up with a list of results, Map is really the clear way to go. Again, try to start thinking of things in terms of the “for each” translations that I provided.

Column/row/regular vectors

Evaluate the examples below, which contain 3x1 and 1x3 matrices as A and B, respectively. Those are often called “column vectors” and “row vectors”, but Mathematica treats them a little differently than just “regular vectors” such as {1,2,3}. Specifically, these column and row vectors are 2D matrices whereas a regular vector is a 1D list of three elements. Usually Mathematica is smart enough to figure out what you mean if you use a “regular vector” where a 3x1 or 1x3 vector technically should be used. But be aware that I’ve found in some cases where you technically should be using 3x1 or 1x3 vectors, Mathematica may actually require regular vectors instead!

A = {{1},{2},{3}};  (* column vector *)
B = {{3,2,1}};  (* row vector *)
MatrixForm/@{A,B,A.B, B.A}

vec={11,12,13}; (* regular vector *)
MatrixForm[vec]

Solving Systems Of Equations

Application: Linear systems of equations

A vitally important application of matrices is to solve systems of linear equations. For example, consider these simulataneous equations: 
2x + 2y + 3z = 4
3x + 4y - 5z = 6
2x - 2y + z = 1

They can be written in this equivalent matrix form:

(({{2, 2, 3}, {3, 4, -5}, {2, -2, 1}}))(({{x}, {y}, {z}}))=(({{4}, {6}, {1}}))

If you’re not familiar with this, take a moment to do some matrix multiplication by hand to verify that the matrix equation really does contain the three simulateous equations. Each row of the multiplied-out matrix is equivalent to one of those three equations.

Written in matrix form, and with a nice program such as Mathematica around to do the matrix inversion, the solution is easily obtained by multiplying both sides by the inverse of the matrix: 

(({{x}, {y}, {z}})) = Inverse(({{2, 2, 3}, {3, 4, -5}, {2, -2, 1}})) .  (({{4}, {6}, {1}}))

M = {{2,2,3},{3,4,-5},{2,-2,1}};
M //MatrixForm
v = {4,6,1};  (* I'm using a regular vector instead of a 3x1 column vector *)
v //MatrixForm

answer=Inverse[M]. v ;  (* IMPORTANT! In Mathematica the inverse of a matrix M is *not* M^-1.  It's Inverse[M] *)
answer //MatrixForm //N

So the solution to the three simultaneous equations above is x = 1.175, y = 0.7125, and z = 0.075.  In this particular case, the Solve command also could have been used to solve the system of equations, but the technique of inverting matrices to solve systems of equations is very generally applicable and one you should know.

Assignment 1. A question of ages.

Alice is twice as old as Bob. Alice's age added to Bob's age is the same as Christa's age. Bob’s age added to Christa’s age is the same as David’s age. The sum of their four ages is 60. How old is each person?  
(a) Write down four equations based on these statements, where Alice's, Bob's, Christa's, and David’s ages are called a, b,  c, and d, respectively. 
(b) Put the equations in matrix form and solve by multiplying by the inverse of a matrix. You should obtain all four ages, which fit all of the given conditions!

Optional: You can also do this type of thing with the LinearSolve command. I believe that command solves linear systems like this via a faster technique than matrix inversion. Try it out if you want. I myself nearly always use matrix inversion, though, even though LinearSolve may be faster, because matrix inversion is fairly intuitive to me and I rarely deal with huge matrices where there could be a noticeable speed difference.

Matrix Properties and Functions

Associative? Commutative?

Matrix multiplication is always associative. That is, for any two matrices a and b,  a.(b.c) = (a.b).c. However, it’s important to remember that matrix multiplication is not in general commutative. That is, only rarely does a.b = b.a. 

You can see that with these three 4x4 matrices containing random numbers:

a = RandomReal[{-10,10},{4,4}];
%//MatrixForm
b = RandomReal[{-10,10},{4,4}];
%//MatrixForm
c = RandomReal[{-10,10},{4,4}];
%//MatrixForm

a.(b.c) == (a.b).c  (* the double equals sign here queries to see if the two sides are equal *)
a.b == b.a  (* if this turns out to be True for anyone using the random matrices, I'll give you a dollar! :-) *)

Sometimes, however, matrix multiplication does turn out to be commutative. That’s always the case, for example, when multiplying two diagonal matrices together (ones with zeroes everywhere except down the diagonal).

a = DiagonalMatrix[RandomReal[{-10,10},4]];
%//MatrixForm
b= DiagonalMatrix[RandomReal[{-10,10},4]];
%//MatrixForm
a.b== b.a   (* Now I'll give a dollar to anyone for whom this doesn't turn out to be true. :-)  *)

Trace and determinant

First, evaluate the cell below to define a 3x3 matrix M for use in this section.

M = {{1,2,3},{3,1,2},{2,3,1}};
%//MatrixForm

The two most important functions which take a matrix as the argument and return a single number are Trace and Determinant. A matrix’s trace is defined as the sum of its diagonal elements. It can be important in eigenvalue problems (see below), since the trace always equals the sum of the eigenvalues. The determinant of a matrix is more complicated to calculate and (when done by hand) is usually calculated by writing the determinant of a large matrix in terms of determinants of its sub-matrices. If you are unfamiliar with that technique, here’s a link to Wikipedia’s page on determinants.  The single-most important property of determinants is that a matrix whose determinant = 0 cannot be inverted. It’s just like you can’t calculate 1/x if x = 0. To rephrase; if det(A) = 0,  there is no matrix B such that B.A = identity matrix. That has very important ramifications in solving simultaneous equations: if the matrix you worked out for “A question of ages” assignment above had had determinant = 0, there would have been no solution.

In other “fun facts”, the determinant of a 3x3 matrix equals the volume of the parallelepiped whose edges are formed by the vectors making up the three columns (or also the three rows) of the matrix.

Assignment 2. Tr and Det calculations.

(a) Use Diagonal to extract the diagonal elements from the matrix M, and then Apply the Plus operation to obtain the sum of the diagonal elements. (The Total function would be easier, but I wanted you to see Apply. Its short cut is @@.)  Verify that this calculation agrees with Tr[M].
(b) Use Det to compute the determinant of M.  Then compute the determinant of (M^(-1)) (meaning the inverse of M, which is not the Mathematica command M^-1).  How are the two determinants related?
(c) Create a 3x3 matrix whose determinant is 0. An easy way to do this is to force two of the column (or row) vectors that make up the matrix to be multiples of each other. Compute its determinant, then try to take the inverse.

Eigenvalues and eigenvectors

The functions Eigenvalues and Eigenvectors give you the eigenvalues and eigenvectors of a matrix. That is, they give you the only possible numbers and vectors which are solutions of the “eigenvalue equation”  λ x = M.x, where λ is a number, x is a generalized “vector” (i.e. not necessarily representing 3D space), and M is the matrix. The matrix must be square. There are the same number of eigenvalues and eigenvectors as there are rows (or columns) of the matrix, although some eigenvalues could be repeated. The solutions are sometimes written as (Subscript[λ, 1]), (Subscript[λ, 2]), (Subscript[λ, 3]), ..., and (Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3]), ..., meaning that (Subscript[λ, 1])(Subscript[x, 1])= M.(Subscript[x, 1]) and so forth. 

The command Eigenvalues tells you the allowed values of the constant λ. 

The command Eigenvectors tells you the allowed values of the vector x (in matrix form, with each allowed vector given as one row of the matrix). Eigenvectors are only defined up to a multiplicative constant. That is, if you have a 2x2 matrix where (1,3) is an eigenvector, then Mathematica might report (1,3) for the eigenvector... or (2,6), or even (1/sqrt(10), 3/sqrt(10)). They are all equivalent.

Each eigenvalue has a particular eigenvector associated with it. Unfortunately, Mathematica doesn’t always sort the eigenvectors returned by the Eigenvectors command in the same order that the eigenvalues are given in results of the Eigenvalues command. The additional function called Eigensystem allows you to see which ones correspond to which.

Eigenvalues and Eigenvectors: What Do They Mean?

The concept of eigenvalues and eigenvectors is a hard one for many students to grasp. I myself had great difficulty understanding what they represented and what they were good for until I learned about coupled harmonic oscillators in a physics class several years after I took linear algebra. I’ve simplified that application into this example; I hope it similarly helps you to understand how eigenvalues and eigenvectors can arise in physics contexts and why they are extremely important to physicists.

Coupled harmonic oscillators: normal modes

The following system of two masses and three springs is an example of coupled harmonic oscillators. The masses happen to be equal and two of the springs happen to have the same spring constant, but that wouldn’t necessarily need to be the case.



Before we solve this system, let’s review how you likely solved the situation of a single mass connected to a single spring in Physics 121. That’s called a simple harmonic oscillator, or SHO. Picture the above diagram with only the left spring and the left mass. When you analyzed the SHO in Physics 121, you undoubtedly did so with Newton’s Second Law: (Subscript[F, net])=m a . The only force acting on the mass is from the spring, a negative kx (from Hooke’s law), so the equation becomes (reversing the usual left and right hand sides):

m(((d^2)x)/(dt^2))=-k x

Then, you (or your teacher) probably guessed a solution of the form x=A cos(ωt), or perhaps x=(Subscript[A,  ])sin(ωt), depending on whether it starts at rest or in motion. Plugging that trial solution into the equation, it’s easy to show that the equation is true as long as ω=(Sqrt[k/m]). Problem solved!

Back to coupled masses & springs. If the middle spring were not present the two masses would each oscillate at ω=(Sqrt[(Subscript[k,  b])/m]). However, the middle spring couples the oscillations, and (as you will see) there are now two different fundamental modes of oscillation at two different frequencies.

To analyze the coupled system, we’ll do it like the SHO and start with Newton’s Second Law--but applied to each mass individually. Let’s call “u” the displacement of the left mass and “v” the displacement of the right mass. The net displacement of the center spring is then v-u. Each spring supplies a force according to Hooke’s Law, (Subscript[F, spring])=kΔx. Here are the two equations (ma written on left, (Subscript[F, net]) written on right), with signs chosen very carefully so that the springs push/pull in the right directions:

m(((d^2)u)/(dt^2))=(Subscript[k, s])(v-u)-(Subscript[k, b])u
m(((d^2)v)/(dt^2))=-(Subscript[k, b])v-(Subscript[k, s])(v-u)

Grouping terms a little differently, we have:

m(((d^2)u)/(dt^2))=u(-(Subscript[k, b])-(Subscript[k, s])) + v((Subscript[k, s]))
m(((d^2)v)/(dt^2))=u((Subscript[k, s]))+v(-(Subscript[k, b])-(Subscript[k, s]))

Since a single mass on a spring oscillates sinusoidally, let’s make a guess and assume that each mass here will just oscillate back and forth in a sinusoidal fashion, both of them at the same ω and each with a constant amplitude (not necessarily the same). This assumption turns out to NOT be true in general, as you will see in the next section. But it is true if the starting positions of the masses are chosen carefully. The oscillation modes for which this assumption is true are called the “normal modes” of the system. 

With this assumption, the displacements u and v will then have solutions of the form u=(Subscript[u, 0])cos(ωt) and v=(Subscript[v, 0])cos(ωt), where (Subscript[u, 0]) and (Subscript[v, 0]) are the amplitudes of oscillation. As with solving the SHO, the next thing to do is to plug those guesses into the equations and see under what conditions the equations can be true.

The time derivative on the left hand side of the top equation becomes: ((d^2)/(dt^2))((Subscript[u, 0])cosωt) =-(ω^2)(Subscript[u, 0])cos(ωt), and the time derivative for the bottom equation is -(ω^2)(Subscript[v, 0])cos(ωt). Multiplying by -1 and dividing by cos(ωt), the simultaneous equations now simplify to:

m (ω^2)(Subscript[u, 0])=(Subscript[u, 0])((Subscript[k, b])+(Subscript[k, s])) + (Subscript[v, 0])(-(Subscript[k, s]))
m (ω^2)(Subscript[v, 0])=(Subscript[u, 0])(-(Subscript[k, s]))+(Subscript[v, 0])((Subscript[k, b])+(Subscript[k, s]))

Or, in matrix form, with λ=m (ω^2):

λ(({{(Subscript[u, 0])}, {(Subscript[v, 0])}})) = (({{(Subscript[k, b])+(Subscript[k, s]), -(Subscript[k, s])}, {-(Subscript[k, s]), (Subscript[k, b])+(Subscript[k, s])}}))(({{(Subscript[u, 0])}, {(Subscript[v, 0])}}))

This is now an eigenvalue equation:  λ x = M.x, where λ is a constant, x is a two element “vector”, and M is a matrix. Not just any λ and x can work! In fact, in this case only two λ’s will yield solutions to the equation. (Remember there are the same number of eigenvalues as there are simultaneous equations.) Find the eigenvalues will tell you what the allowed oscillation frequencies are: one allowed ω for each λ (because ω=(Sqrt[λ/m])). 

Each eigenvalue & frequency will also have a particular x (eigenvector) associated with it. In this example x = ((Subscript[u, 0]),(Subscript[v, 0]))... that means the eigenvectors give you information about the amplitudes of oscillation that correspond to the allowed frequencies. As mentioned above, eigenvectors can always be multiplied by a constant, so if a particular (u,v) combination “works”, then any constant times that (u,v) pair will also work and will produce the same λ--so the eigenvectors really only tell you about the ratios of the amplitudes, not about their absolute values. As an example, this means that if you get (1,3) as an eigenvector, then for that mode of oscillation the second mass would always have three times the amplitude of the first mass. It wouldn’t matter if the amplitudes of oscillation were 1 cm vs. 3 cm, 2 cm vs. 6 cm, etc., all such amplitudes will give rise to the same frequency. Or, if you get (-1,3) as an eigenvector, then that means the oscillations are out of phase: the ratio of amplitudes is still 3 to 1, but one is moving to the left while the other is moving to the right.

Enough discussion, let’s solve the problem already!

Assignment 3. Coupled harmonic oscillators, part 1

For simplicity, let's set m = 1 and (Subscript[k, s])=(Subscript[k, b])=1 for some choice of units. Use Mathematica to solve the eigenvalue equation and determine what two frequencies are allowed, along with what the relationship is between (Subscript[u, 0]) and (Subscript[v, 0]) for each allowed frequency. Compare the two frequencies of the two modes to the frequency the masses would have if they were not coupled together. Qualitatively explain why the allowed frequencies are higher, lower, or the same as the uncoupled case.

By the way, notice how easily this approach could straight-forwardly be applied to a MUCH larger system of coupled masses and springs. You just have to be able to specify the appropriate “spring constant matrix”. In my Physics 581 class, I had students solve the problem for 100 masses and 101 springs, where there are 100 allowed frequencies. The answer to that problem turns out to be very closely connected to the allowed vibrational modes in a solid (you can think of the masses as being vibrating atoms and the springs as being the bonds between the atoms).

Coupled harmonic oscillators: time dependence

When we were doing the analysis of the coupled harmonic oscillators just above, we made a critical assumption in saying that the two masses would be oscillating sinusoidally. Such a situation occurs if the masses are released from rest in just the right fashion to force them into one of the two “normal modes” (one of the oscillation modes described by the eigenvectors you just found). In general, however, this might not occur... and the motion then won’t be purely sinusoidal.

Solving the general situation is actually very similar to solving the normal mode situation, and we can use most of the work we did above. Returning to the equations just before the critical assumption, we have:

m(((d^2)u)/(dt^2))=u(-(Subscript[k, b])-(Subscript[k, s])) + v((Subscript[k, s]))
m(((d^2)v)/(dt^2))=u((Subscript[k, s]))+v(-(Subscript[k, b])-(Subscript[k, s]))

In matrix form the equations are:

m ((d^2)/(d (t^2)))x + M·x = 0
 
where m is the mass, x = (({{u}, {v}})) and M = (({{(Subscript[k, b])+(Subscript[k, s]), -(Subscript[k, s])}, {-(Subscript[k, s]), (Subscript[k, b])+(Subscript[k, s])}})), the same matrix we used above.
 
Here’s some advanced code to solve that matrix equation. It uses DSolve, which you don’t need to know how to use for this course--so focus on the results, not the code. The output is horrendous because we’ve left it very general--no values for the spring constants and no values for the initial conditions. Bear with me for now because pretty graphs are coming up next.

Clear[m,kb,ks,u,v,M] (* in case these variables were set to values above *)
M = {{kb+ks,-ks},{-ks,kb+ks}};  (* the matrix of spring constants *)
unknowns ={u[t],v[t]};    (* This is the vector x, containing u(t) and v(t) *)
equations = (#==0&)/@ (m D[unknowns,{t,2}]+M.unknowns);
initialconditions = {u[0]==u0,u'[0]==uvel0,v[0]==v0,v'[0]==vvel0};  (* u0 and uvel0 are initial position and velocity for mass 1; v0 and vvel0 are same thing for mass 2 *)
positions = unknowns/.DSolve[{equations,initialconditions},unknowns,t][[1]]

Here is some code that will use the ugly output from DSolve to plot the displacements of the two masses, u and v. This specific code is for the situation where (Subscript[k, s]) = very small compared to (Subscript[k, b]) (”weak coupling”) and where the left mass has some initial displacement but the right mass does not. Neither has any initial velocity. Notice how there’s a beating effect, and the energy from the left mass gets transferred to the right mass and vice versa. This is an indication that it’s not oscillating in one of the normal modes--in normal mode oscillation both masses oscillate purely sinusoidally.

timedep = positions /.{m->1,kb->1,ks->0.1,u0->1, v0->0,uvel0->0,vvel0->0}//Re//ComplexExpand//Chop;
Plot[timedep,{t,0,100},PlotStyle->{{Thick,Red},{Thick,Blue}}]

Assignment 4. Coupled harmonic oscillators, part 2

Modify the code for the previous graph so that the springs match the situation for which you found the normal modes above (“strong coupling”). 
(a) Play around with the initial conditions and make plots of the displacement of the two masses for a few different initial conditions. Among other things, try the case where only one mass is given an initial displacement and neither has any initial velocity, as I just did for the weak coupling case above.
(b) Now give both masses initial displacements as in normal mode 1. Verify that the period you observe matches your predictions of the normal mode frequency from the last assignment. In case you’ve forgotten, period = 2π/ω. Change the time scale as appropriate. 
(c) Repeat, giving both masses initial displacements as in normal mode 2. 

No other set of initial conditions other than the ones you just tried will give rise to pure sinusoidal behavior (well, none with no initial velocity, anyway). Try it! And no purely sinusoidal behavior will exist at any other frequencies than the ones predicted by the normal modes. And the eigenvalues and eigenvectors are directly connected to the normal mode frequencies and patterns of oscillation.

Rotation Matrices

2D Rotation matrices

When a 2x2 matrix operates on a 2 component vector (via standard matrix multiplication) or when a 3x3 matrix operates on a 3 component vector, the result is another 2 or 3 component vector. If the resulting vector has the same length as the original vector, the matrix multiplication can be visualized as having rotated the vector to a new angle. Matrices that do this are called “rotation matrices”.

In 2 dimensions, the standard form for a rotation matrix is this, with phi being the counter-clockwise angle of rotation in the x-y plane.

R2d = {{Cos[phi] , -Sin[phi]},{Sin[phi],Cos[phi]}};
% //MatrixForm

Remember that Sin and Cos in Mathematica require angles to be in radians, not degrees, so use Degree or /Degree as appropriate.

Incidentally, rotation matrices are part of a special class of matrices called “unitary” matrices because all of their eigenvalues (which can be complex numbers) have a magnitude of 1. Unitary matrices are particular important in some quantum mechanics applications, so you’ll see eigenvalues in that context in Phys 451 and/or 452 if you take those classes in the future.

Eigenvalues[R2d/.phi->RandomReal[360] Degree ]   (* eigenvalues of a rotation matrix containing a random angle between 0 and 360 *)
% //Abs  (* Abs gives the magnitude of complex numbers. Both eigenvalues have magnitudes of 1. *)

Assignment 5. 2D rotation by an angle

(a) Create a 2x1 column vector and call it v. It can be any vector you want. Calculate its length (via Norm). Then operate on your vector with the R2d matrix, using some angle between 30 and 90 deg. That is: vnew =  R2d . v.  Calculate the length of the resulting vector vnew... it should have the same length as the vector you started with!
(b) Use the Graphics and Arrow commands from Lab 3 to graphically view the vector rotation. Use the options Axes->True, AspectRatio->1, and use PlotRange to have the same scale on the x- and y-axes. You will likely need to use Flatten to get v and vnew into the right format for the Arrow command (look at the help page on Flatten to see what it does). Use different colors for your vectors.
(c) Verify that the angle (from the x-axis) for vnew compared to v has indeeed increased by the angle you chose. Use one of these two methods:  ArcTan with a pair of numbers, or Arg with a complex argument.
(d) Verify these properties of rotation matrices for a general two-dimensional rotation matrix (no angle specified): (1) det R = 1.  (2) inverse(R) = transpose(R). The Simplify command might be useful.

3D rotation matrices

In 3 dimensions, the form of the rotation matrix for arbitrary angles is much more complicated to define, and there are varying conventions about how the angles should be specified. You can look at the Wikipedia article on rotation matrices if you are interested.

Although we won’t tackle 3D rotation matrices in general, here are two specific 3D rotation matrices for you to consider. (Subscript[R, 1]) performs a 120° rotation around an axis in the (1,1,1) direction, which permutes the vector’s x, y, and z values. (Subscript[R, 2]) performs a 180° rotation around the z axis, which inverts the x and y values. Make sure you can visualize why those rotations perform those specific effects.

R1 = {{0,0,1},{1,0,0},{0,1,0}};(* 120° rotation around {1,1,1} direction  *)
MatrixForm[%]
R2 = {{-1,0,0},{0,-1,0},{0,0,1}};(* 180° rotation around {0,0,1} direction *)
MatrixForm[%]

You can see the effects mentioned by doing matrix multiplication with a general (x,y,z) vector:

R1. {x,y,z}

R2 . {x,y,z}

These types of matrices can also be created directly by Mathematica using the RotationMatrix command: (notice RotationMatrix can also be used to produce 2D rotation matrices)

R1 = RotationMatrix[2 Pi/3, {1,1,1}];
MatrixForm[%]
R2 = RotationMatrix[Pi, {0,0,1}];
MatrixForm[%]

“Similar” Matrices

Similarity transforms

Rotation matrices can be important in what are called “similarity transformations”. Two matrices (Subscript[M, 1]) and (Subscript[M, 2]) are called “similar matrices” if (Subscript[M, 2])=(S^(-1))(Subscript[M, 1])S, for some matrix S. The process of multiplying by (S^(-1)) on the left and by S on the right is called a transformation because it transforms (Subscript[M, 1]) into a different matrix. Although the matrix S doesn’t necessarily need to have any special properties beyond being invertible, if it is a rotation matrix then the similarity transformation has a special meaning as outlined in the next paragraph.

Consider the product of a matrix M and a vector v, namely w = M.v. A way of interpreting this equation is as follows: the matrix product M.v produces a new vector, w, which is a “linear transformation” of v. That means that the three components of w are linear combinations of the three components of v. We sometimes call M a linear “operator” because it operates on v to produce w. Sometimes a linear operator M is known in one coordinate system (“1”) but not in another (“2”). To find the linear operator in the new coordinate system, you just have to do this similarity transform: (Subscript[M, 2])=(S^(-1))(Subscript[M, 1])S, where S is the rotation matrix that rotates the new x-y-z axes into the old x-y-z axes. Conceptually, that means that the effect of (Subscript[M, 2]) on a vector is this (going from right to left):  (a) S rotates the vector (in the new coordinate system) back into the old coordinate system, (b) the linear operator (Subscript[M, 1]) is applied to transform the vector (in the old coordinate system), then (c) then (S^(-1)) rotates the transformed vector back into the new coordinate system.

Matrix invariants

As just discussed, we can often think of a pair of similar matrices as being the same matrix operation displayed in two different coordinate systems. This means that many important matrix properties are invariant under similarity transformations (i.e. they are the same for both (Subscript[M, 1]) and (Subscript[M, 2])). Here, you will explore a few of these invariant properties.

Assignment 6. Some matrix invariants

I have constructed a matrix (Subscript[M, 1]) in the cell below for you to use in this problem.

M1 = {{7,2,3},{2,6,2},{3,2,1}}

(a) Construct a 3×3 matrix of random real numbers named S and use it to calculate the matrix (Subscript[M, 2])=(S^(-1))(Subscript[M, 1])S, which will be similar to (Subscript[M, 1]) by definition. Hint: you can construct a matrix of random real numbers with a single RandomReal command, as I did above in the “Associative? Commutative?” section.
(b) Compute the traces of both (Subscript[M, 1]) and (Subscript[M, 2]). Are they equal?
(c) Compute the determinants of both (Subscript[M, 1]) and (Subscript[M, 2]). Are they equal?
(d) Compute the eigenvalues of both (Subscript[M, 1]) and (Subscript[M, 2]). Are they equal?
(e) Compute the eigenvectors of both (Subscript[M, 1]) and (Subscript[M, 2]). Are they equal? (The eigenvectors are the rows of the matrix returned by the Eigenvectors command.)

Diagonalization

There is often a similarity transformation that can transform a matrix (Subscript[M, 1]) into a diagonal matrix (one that is zero everywhere off the diagonal).  As you saw in the previous section, the eigenvalues are unchanged in a similarity transformation, so the eigenvalues of  (Subscript[M, 1]) are the same as the eigenvalues of this diagonal matrix.  Also, as you should know, the eigenvalues of a diagonal matrix are just the numbers along the diagonal. Therefore you can produce a diagonal matrix similar to (Subscript[M, 1]) by finding the eigenvalues of (Subscript[M, 1]) and placing them along the diagonal, even if you have no idea what the matrix S actually is that transforms (Subscript[M, 1]) into the diagonal matrix via  (S^(-1))(Subscript[M, 1])S.  Producing a diagonal matrix that contains the eigenvalues of (Subscript[M, 1]) along the diagonal is often called diagonalizing the matrix.

It turns out that such a similarity transformation exists if and only if the eigenvectors are linearly independent, meaning that none of them can be written as linear combinations of the others. And in fact, the matrix S turns out to be a matrix whose columns are constructed from the eigenvectors--which is exactly the Transpose of what the Eigenvectors command returns.

Assignment 7. Two methods for diagonalization

(a) Find the eigenvalues of the matrix M below and use them to produce a diagonal matrix which is similar to M.  No copying & pasting allowed. The command DiagonalMatrix will be helpful. (Note that the eigenvalues can be complex numbers.)
(b) Diagonalize M by defining the correct matrix S from the result of the Eigenvectors command (re-read the paragraph just above this assignment), and calculating (S^(-1)).M.S. The result should match part (a). I found the Chop command to be helpful to remove the very small numbers that appeared in the off-diagonal elements due to internal rounding errors.
(c) Oddly enough, Mathematica does not have a built-in command to diagonalize a matrix. Let’s make one! Use Module to create a function called diagonalize, which takes a square numerical matrix as input and returns a diagonalized version of the matrix as output. Use either of the two methods from (a) or (b) as the guts for your Module. Test your function on the matrix M.

M = RandomReal[{-10,10},{3,3}];
% //MatrixForm

Additional Assignments

Rotated plots

In a previous lab we talked briefly about using ContourPlot to plot conic sections. For example, the parabola y=(x^2)-2 can be plotted like this:

ContourPlot[y==x^2 - 2, {x,-3,3},{y,-3,3}]  (* note the double equals sign *)

The 2D rotation matrices described above can be used to rotate such plots. To rotate that parabola a certain angle about (0,0), I would (a) find the appropriate rotation matrix, (b) apply the matrix to {x,y} to find the new y and x values in terms of the old x and y values, then (c) substitute the new x and y values into the old equation for the parabola to get a new equation for a parabola. Here it is all together, rotating the parabola 30 degrees CCW:

m = RotationMatrix[-30 Degree];  (* see below as to why this is -30 deg *)
{newx, newy} = m . {x,y};
newlhs = y /. y-> newy
newrhs = x^2 - 2  /. x-> newx
ContourPlot[newlhs==newrhs, {x,-3,3},{y,-3,3}]

Things to notice: 
(1) The reason -30° needed to be used in the rotation matrix, is because the parabola that has been rotated 30° CCW, is the one that has the property that when rotated BACK by 30° (i.e. 30° CW), it satisfies the equation, y=(x^2).

(2) This is no longer a “function” any more, in that some x-values give you more than one y value. That’s why ContourPlot needed to be used rather than Plot.

(3) ContourPlot is apparently very picky about needing to have the entire function (including double equals sign) defined inside the function call itself. That’s why I had to define the left hand side and right hand side separately. For example, this command doesn’t work, even though it seems to me like it should:

ContourPlot[y == x^2  /. {x->newx, y-> newy}, {x,-3,3},{y,-3,3} ]   (* doesn't work! *)

Assignment 8. A rotated ellipse

(a) ContourPlot the ellipse, (x^2)+((y^2)/8)=1.
(b) Rotate the ellipse 10° counter-clockwise, and plot the rotated ellipse.

Assignment 9, Extra credit. Rotated eigenvectors

This is a continuation of Assignment 6. Although the eigenvectors of (Subscript[M, 1]) were not equal to the eigenvectors of (Subscript[M, 2]), they are closely related. Show that the three eigenvectors of (Subscript[M, 2]) are equal to the three eigenvectors of (Subscript[M, 1]), just rotated by S.  (Your random matrix S will not technically be a rotation matrix, but treat it as if it were.)

Wrap-up

The following functions were discussed in today's lab, or are very closely related:
Getting rid of small numbers that you know should be zero: Chop
Applying the same function repeatedly to matrices grouped in a list: Map
Matrix functions: Inverse, Tr, Det, Dot, Transpose,  Diagonals
Adding lists: Apply with Plus, or Total
Producing specific matrices: RotationMatrix, DiagonalMatrix 
Solving simultaneous equations: multiplying by a matrix Inverse, or LinearSolve
Eigen stuff: Eigenvalues, Eigenvectors, Eigensystem

Potentially useful summary/tutorial pages:
(I don’t think there were any mentioned today.)