Phys 263 session for 10 May 2012

BTM stuff

Let's warm up with some Shankar homework.
BTM 9.5.2  asks about the rotation matrix

RotationMatrix[θ]//MatrixForm

More precisely, it asks about the 3D version of this.
Your first task: read about RotationMatrix, and ask for the one which rotates about the z-axis by θ.
Then ask for the Eigensystem.

And BTM 9.5.3 asks you to find the eigenvalues and eigenvectors for some 3x3 matrices. e.g. this one:

(M = {{1,3,1},{0,2,0},{0,1,4}})//MatrixForm

To do this "by hand", you first form the characteristic polynomial:

p[λ_] = Det[M - λ IdentityMatrix[3]]

Then you need find the roots.  This may seem hard "by hand", if you get one root, then you can divide
that one out and get a quadratic.   Here you might notice that λ=1 is a root (i.e. 8-14+7-1==0).  Then...

Simplify[p[λ]/(λ-1)]

Or just ask Mathematica for it:

Solve[p[λ]==0]

To find the eigenvectors we pick one of the eigenvalues, e.g. 2, and write down the equations the components must satisfy:

vec = {x,y,z}

eqns = M.vec == 2 vec

And solve:

s=Solve[eqns]

Note the warning: the system is underdetermined.  But the result is still useful.
Choosing, e.g. z=1, we get:

vec /. s[[1]] /. z->1

Or again, we just ask for

Eigensystem[M]

2 masses / 3 springs

And the utility of eigenvectors comes to the fore when we study systems of springs.
Here is a function to draw it:

state[x_,y_,L_:12]:= (Show[{ParametricPlot[{0+(t*x),Sin[2 Pi 6 t]},{t,0,1},Axes->True,PlotStyle->Red,PlotRange->{{0,L},{-1,1}}],Graphics[Rectangle[{x-1/4,-1/4},{x+1/4,1/4}]],
ParametricPlot[{x+t*(y-x),Sin[2 Pi 6 t]},{t,0,1},PlotStyle->Blue],
Graphics[Rectangle[{y-1/4,-1/4},{y+1/4,1/4}]],
ParametricPlot[{y+t*(L-y),Sin[2 Pi 6 t]},{t,0,1},PlotStyle->Orange]}])

Here is what we'll take to be the equilibrium:

state[4,8]

Your task: draw force diagrams and complete this system of equations:

eq1 = m1 x1''[t]  == -(k1+k2)  x1[t] + k2 x2[t]
eq2 = m2 x2''[t]  ==  +k2 x1[t] - (k2+k3) x2[t]

Start with the simplest case, with identical masses and springs:

p= {k1->1,k2->1,k3->1,m1->1,m2->1}

{eq1,eq2}/.p //MatrixForm

Let's solve with certain boundary conditions and solve for the motion:

s=DSolve[{eq1,eq2,x1[0]==1,x2[0]==0,x1'[0]==0,x2'[0]==0}/.p,{x1[t],x2[t]},t]

And then Animate:

y1[t_]=x1[t]/.s[[1]]
y2[t_] = x2[t]/.s[[1]]

Animate[state[4+y1[t],8+y2[t]],{t,0,40}]

Your first task: play with the boundary conditions until you get the "nice" modes with only one frequency.

Next task: find these "nice" modes via matrix technology.
For the case at hand the matrix of couplings is:

p= {k1->1,k2->1,k3->1,m1->1,m2->1}

(m = {{-(k1+k2)/m1,k2/m1},{k2/m2,-(k2+k3)/m2}}/.p) // MatrixForm

Ask for the Eigenvectors etc.  Observe that the "nice" modes are exactly the ones you found above (we assume), i.e.
one mode where the displacements are equal, and one where they are opposite.

Tweak #1--change the middle spring

BTM 9.6.2 asks what happens when the middle spring is changed to k2 = 2 k.  What does happen?
First go up and change the parameter in the rule above and just see what happens.
Then appeal to your intuition, and note that if we launch with equal displacements, the middle spring never gets
stretched, so the precise value of its spring constant is irrelevant.  And if we launch with opposite displacements,
we'll oscillate with the higher frequency appropriate to a mass coupled to two springs, one of  strength k and
the other of 4k (i.e a 2k spring cut in half).
And of course find this result by forming the matrix and asking for Eigenvalues etc.:

p1= {k1->1,k2->2,k3->1,m1->1,m2->1}

(m = {{-(k1+k2)/m1,k2/m1},{k2/m2,-(k2+k3)/m2}}/.p1) // MatrixForm

Eigensystem[m]

s=DSolve[{eq1,eq2,x1[0]==1,x2[0]==-1,x1'[0]==0,x2'[0]==0}/.p1,{x1[t],x2[t]},t]

s=DSolve[{eq1,eq2,x1[0]==1,x2[0]==+1,x1'[0]==0,x2'[0]==0}/.p1,{x1[t],x2[t]},t]

Tweak #2--change one of the masses

Morin 4.9 (solved in the book) asks about the 2 mass case where m1 = m and m2 = 2 m.  Modify the rule p accordingly,
and ask for the movie and the Eigensystem.

p2= {k1->1,k2->1,k3->1,m1->1,m2->2}
(m = {{-(k1+k2)/m1,k2/m1},{k2/m2,-(k2+k3)/m2}}/.p2) // MatrixForm

Eigensystem[m]

s=DSolve[{eq1,eq2,x1[0]==-1-Sqrt[3],x2[0]==1,x1'[0]==0,x2'[0]==0}/.p2,{x1[t],x2[t]},t]

FullSimplify[s,Assumptions->t>0]

Tweak #3--remove one spring

Morin 4.31 asks about the 2 mass case where there is no third spring.  This is accomplished by changing k3 to 0.
Your task: modify the equations and get the movie.  Here is the matrix again:

p3= {k1->1,k2->1,k3->0,m1->1,m2->1}
(m = {{-(k1+k2)/m1,k2/m1},{k2/m2,-(k2+k3)/m2}}/.p3) // MatrixForm

Eigensystem[m]//Simplify

s=DSolve[{eq1,eq2,x1[0]==(-1-Sqrt[5])/2,x2[0]==1,x1'[0]==0,x2'[0]==0}/.p3,{x1[t],x2[t]},t]

s3=FullSimplify[s,Assumptions->t>0]

And here is a drawing function sans spring 3:

state2[x_,y_]:= (L=12;Show[{ParametricPlot[{0+(t*x),Sin[2 Pi 6 t]},{t,0,1},Axes->True,PlotStyle->Red,PlotRange->{{0,L},{-1,1}}],Graphics[Rectangle[{x-1/4,-1/4},{x+1/4,1/4}]],
ParametricPlot[{x+t*(y-x),Sin[2 Pi 6 t]},{t,0,1},PlotStyle->Blue],
Graphics[Rectangle[{y-1/4,-1/4},{y+1/4,1/4}]]}])

q1[t_] = (x1[t]/.s3[[1]])
q2[t_] = (x2[t]/.s3[[1]])
Animate[state2[4+q1[t],8+q2[t]],{t,0,20}]

Tweak #4--damping

Morin 4.34 asks you to add damping, i.e. a friction force - b x'[t] = -2 γ m x'[t].
Choosing some reasonably small γ (i.e. underdamping), add the -2γ term to the systems above and watch them slow down.
Observe that since the damping is the same for both masses the system has the same eigenmodes---check it out by
playing with the boundary conditions.
To spell it out, if we take the usual ansatz x(t) = x0 Exp[I β t], the matrix equation becomes:

-β^2 x0 + (2 γ) I β x0 + M.x0 ==0

i.e. if x0 is an eigenvector of M with eigenvalue ω0^2, we need to choose β to satisfy the usual damped oscillator quadratic:

-β^2 + 2 I β γ + ω0^2 ==0

More generally, if the damping were different for the two masses, or if it depends on the relative coordinates (as would happen
e.g. if two masses are coupled by a spring and a dashpot) then we get a more complicated matrix equation:

-β^2 x0 + 2  I β Γ.x0 + M.x0 ==0

where Γ is the matrix of damping terms.  To solve we'd need to find the eigenvalues of the matrix (M + 2IβΓ), leaving β as an unknown symbol, and then solve the (potentially complicated) algebraic equation for the frequencies β.

Tweak #5--more masses

Your final task: add a third mass.  That is, write the equations for three coordinates, DSolve, and animate with:

state[x_,y_,z_]:= (L=12;Show[{
ParametricPlot[{0+(t*x),Sin[2 Pi 6 t]},{t,0,1},Axes->True,PlotStyle->Red,PlotRange->{{0,L},{-1,1}}],Graphics[Rectangle[{x-1/4,-1/4},{x+1/4,1/4}]],
ParametricPlot[{x+t*(y-x),Sin[2 Pi 6 t]},{t,0,1},PlotStyle->Blue],
Graphics[Rectangle[{y-1/4,-1/4},{y+1/4,1/4}]],
ParametricPlot[{y+t*(z-y),Sin[2 Pi 6 t]},{t,0,1},PlotStyle->Orange],
Graphics[Rectangle[{z-1/4,-1/4},{z+1/4,1/4}]],
ParametricPlot[{z+t*(L-z),Sin[2 Pi 6 t]},{t,0,1},PlotStyle->Green]}])

eq1 = m1 x1''[t]  == -(k1+k2)  x1[t] + k2 x2[t]
eq2 = m2 x2''[t]  ==  +k2 x1[t] - (k2+k3) x2[t] + k3 x3[t]
eq3 = m3 x3''[t] == k3 x2[t] - (k3+k4) x3[t]

Start with the simplest case, with identical masses and springs:

p4= {k1->1,k2->1,k3->1,k4->1,m1->1,m2->1,m3->1}

system=({eq1,eq2,eq3,x1[0]==1,x2[0]==0,x3[0]==0,x1'[0]==0,x2'[0]==0,x3'[0]==0}/.p4)

s=DSolve[system,{x1[t],x2[t],x3[t]},t]

s4 = Simplify[s[[1]]]

s=DSolve[{eq1,eq2,eq3,x1[0]==1,x2[0]==-Sqrt[2],x3[0]==1,x1'[0]==0,x2'[0]==0,x3'[0]==0}/.p,{x1[t],x2[t],x3[t]},t]

y1[t_]=x1[t]/. s4
y2[t_] = x2[t]/.s4
y3[t_] = x3[t]/.s4

Animate[state[3+y1[t],6+y2[t],9+y3[t]],{t,0,50}]

Note that to find the "nice" starting conditions which excite only one frequency,
you'll want to find the eigenvectors of the appropriate 3x3 matrix.

Generalizing to more masses

For more masses we need even bigger matrices, and we don't want to bother writing them out explicitly. 
Instead ... here is a useful Mathematica function:

?If

Here is how one might use it, e.g. to make an absolute value function:

f[x_] := If[x>0,x,-x]

Plot[f[x],{x,-2,2}]

And here we use it to construct the matrix elements of our spring system for a larger number of masses.
First pick a number of masses:

num = 12

Next decide if you want fixed or periodic boundary conditions, i.e. whether masses 1 and N
have springs which attached to a fixed wall as usual, or if there is no wall, and the masses
and springs form a circle. Try one and then the other.

(*fixed wall*) m[i_,j_] = If[i==j,2.0,If[(i==j+1)||(i==j-1),-1.0,0]]

(* periodic *)m[i_,j_] = If[i==j,2.0,If[(i==j+1)||(i==j-1)||(i==num&&j==1)||(i==1&&j==num),-1.0,0]]

(mm = Table[ m[i,j],{i,1,num},{j,1,num}])//MatrixForm

For full generality, let's leave a hook where we can put in damping later:

γ=0.0

And define a list of variables we want to solve for:

(vars = Table[x[i,t],{i,1,num}])//MatrixForm

A list of equations of motion for them:

(eq =Table[D[x[i,t],{t,2}]== -Sum[m[i,j] x[j,t],{j,1,num}]-γ D[x[i,t],t],{i,1,num}])//MatrixForm

A list of boundary conditions saying that everyone starts at rest:

bc1=Table[((D[x[i,t],{t,1}])/.t->0)== 0,{i,num}]

And that almost everyone starts at equilibrium:

bc2 = Table[ x[i,0]== 0,{i,2,num}]

Lastly we pick one special mass (the first one) which is displaced:

bc3 = {x[1,0]== 1}

Then we can solve---might as well do it numerically for efficiency:

s1=NDSolve[Join[eq,bc1,bc2,bc3],vars,{t,0,30}][[1]];
Do[xx[i,t_]= Re[x[i,t]/.s1],{i,num}]

Here is a low-tech way of plotting---just make a ListPlot:

state[t_] := ListPlot[Table[xx[i,t],{i,1,num}],PlotStyle->{PointSize[.03],Hue[.6]},PlotRange->{{0,num+1},{-1,1}}]

state[0]

Animate[state[t],{t,0,10,.1}]

The "nice" modes are of course the eigenvectors of our matrix:

(eig=Eigenvectors[mm])//Transpose//MatrixForm

Let's change the starting conditions, picking mode number 10 in the list
Obviously you might change the 10 to something else and explore.

bc4 = Table[ x[i,0]== eig[[10,i]],{i,1,num}];
s2=NDSolve[Join[eq,bc1,bc4],vars,{t,0,30}][[1]];
Do[xx[i,t_]= Re[x[i,t]/.s2],{i,num}] ;
Animate[state[t],{t,0,30,.2}]