Simple Harmonic Motion

(copyright Kieran Mullen 2007)

In order to become a bit more familiar with simple harmonic motion (SHM) 
	it helps to plot a variety of different functions, and investigate the nature of the  		olutions to the dquation of motion.  First, let's fire  up the kernel by evaluating the
	simple arithmetic statement below.

2+2

You should make your screen as wide as possible, by clicking the mouse button on
	the lower left hand corner of this window, and stretching the window to fit your 
	screen.

It pays to first become a bit familar with the trigonometric functions.  Recall our
	basic starting points are the sine and cosine functions.  These are plotted in the next
	commands.

Plot[Sin[t],{t,0,4 π}]

Plot[Cos[t],{t,0,4 π}]

Note that the fundamental period of each function is 2Pi.  We have plotted two 
 	periods of each.  We can put both curves on the same graph, by placing the 
 	functions to be plotted in a list, {fxn1[t], fxn2[t]......}.  Below is an example of
 	this.

Plot[{Cos[t],Sin[t]},{t,0,4 π}]

When we consider oscillatory motion, we often refer to the amplitude, the
	phase, and the frequency of the oscillation.  For example, if we consider the
	function:
	
	        a   Sin[ w t + phi]
	       
	 The amplitude is "a", the phase is "phi" and the frequency is "w".

(1) Plot the following three graphs over the range 0 to 4Pi.  (each graph will have
	two curves on it). 
  	    (a) Sin[t] and a sine function with twice the amplitude.
   	    (b) Sin[t] and a sine function with twice the frequency.
  	    (c)  Sin[t] and a sine function with a  phase shift of Pi/6

As we have discussed in class, the sine and cosine are handy functions because
	their derivatives are sines and cosines.  You may verify your understanding of
	Mathematica by calculating the following derivatives:

(∂_t)Sin[w t]

(∂_t)Cos[w t]

From these two results we can see that the second derivatives of sine and
	cosine return the negative of the original function:

(∂_t)((∂_t)Sin[w t])

(∂_t)((∂_t)Cos[w t])

This means that the solution of the differential equation:

                  h''[t] + w^2 h[t] = 0
                  
	is that h[t] can be either the sine function, or a cosine function, or some linear
	combination of the two.   We can see this by asking Mathematica to solve the
	differential equation of a simple harmonic oscillator:

DSolve[h''[t]+(w^2) h[t]==0,h[t],t]

You will notice the constants C[1] and C[2].  They appear because the above 
	differential equation has two solutions, namely the sine and the cosine.  Any 
	combnation of the form:

   		   C[1]  Sin[w t]  + C[2] Cos[w t]
      
	for any constants C[1] and C[2] will be a solution.  We can only determine the
	constants by knowing more information, called the boundary conditions.  We 
	usually specify propertes of the desired solution at the start, t=0.  These are 
	called initial conditions. Below we look for the solution that has a value of 0
	at t=0, and next we look for the solution with a derivative of zero at t=0.

Simplify[DSolve[{h''[t]+(w^2) h[t]==0,h[0]==0},h[t],t]]

Simplify[DSolve[{h''[t]+(w^2) h[t]==0,h'[0]==0},h[t],t]]

Note that we put the initial conidition and the differential equations together in
	a list "{...}".  The DSolve function tries to solve this as a system of equations.

Ok.  Here is where the  fun begins.  We want to solve not just
	plain harmonic oscillators (springs) but we want to add damping.  Note we've
	added a damping term (g h'[t]) to our differential equation.  The constant "g"
	is the strength of the damping, larger g means more friction.  Remember that "w"
	is the natural frequency of our oscillator, the angular frequency it would have if
	there were no damping.   Shift-Retrun on the expression below to get the 
	general solution.

DSolve[h''[t]+g h'[t]+(w^2) h[t]==0,h[t],t]

Again we see our friends C[1], and C[2].  We still have two independent solutions
	to our differential equation.  We can see that there are exponentials buried in there,
	too.  If the terms inside the Sqrt are negative, we'll get a factor of "I", and this will
	give us oscillations. 

 	    We can pick out two independent solutions by again looking for a solution
 	that starts out with zero value, but non-zero derivative, and a second that starts
 	out with a non-zero value, but zero derivative.  Look at the expressions below,
	and make sure you understand what they are doing.  Then shift-return on them.

DSolve[{h''[t]+g h'[t]+(w^2) h[t]==0,h[0]==0,h'[0]==1},h[t],t]

DSolve[{h''[t]+g h'[t]+(w^2) h[t]==0,h[0]==1,h'[0]==0},h[t],t]

Ugh!  Let's define this a function, and plot it, so we can see what is going on.  We	choose the solution that starts at h[0]=1.  We will call this solution H[t, g, w].   Don't  worry about the complexity of the next line.  We are persuading Mathematica to take the output of one  command (DSolve) and turn it into a new function.  This takes a bit 	of finesse.  (The command translates as "Define the function H[t,g,w] to be the function h that is the first solution of the differential equation X with boundary conditions Y."  It happens to be the only solution, not just the first, but Mathematica  is a bit stupid sometimes.)

H[t_,g_,w_]=h[t]/.DSolve[{h''[t]+g h'[t]+(w^2) h[t]==0,h[0]==1,h'[0]==0},h[t],t][[1]]

Below is a command to plot the function H over the range t=0 to t= 6Pi, for a 	damping (g) of  zero and a natural angular frequency (w) of 1. We've also told 	Mathematica to limit the  vertical axis to be between -1 and 1.

Plot[H[t,0,1],{t,0,6 π},PlotRange->{-1,1}]

(2)  Copy the above commad and replot the function over the same interval and keeping
	the angular frequency w=1.  Choose a different values of g such that the quantity 
	
	                g^2 - 4 w^2
	                
	 (a)  is less than zero.
	 (b) is very close to, but not equal to zero.
	 (c) is greater than zero.
	 
	 Which case is underdamped, which case is overdamped, and which case is
	 critically damped?   Assume that H is the height of an object tied to a spring
	submerged in water.   Give a physical interpretaion of the graphs you plotted.  
	
	(3) Use Plot3D to plot a range of solutions for g between 0 and 3, and t over
	the range  0 to 4 Pi .  Explain how to interpret this surface plot.

Ok.  Now the REAL fun begins.  The damped harmonic oscillator will eventually
 come to rest,  The friction term consumes energy, until the object eventually comes to
 rest.  In order to avoid this sedentary fate, we can drive the oscillator with a periodic
 driving force.   This driving force continues to do work, putting energy back into the 
 system.   Note however that we can drive the system at a frequency that is different
 from its natural, undamped frequency.  Below is the general equation for such a system.
 Shift-return on it, and wait for a bit.  This will take a minute or two. 
 
      While you wait, here's a bit of trivia.  This differential equation is described
  as "linear" because it depends only on a function and its derivatives, but never the
  square of the function, or the square of a derivative.  It is also described as 
  inhomogeneous", because the right hand side is a function that doesn't depend on
  y[t] or any of its derivatives.  Inhomogeneous equations are usually much more of
  a pain to solve.

DSolve[y''[t]+g y'[t]+(w^2) y[t]==Sin[w0 t],y[t],t]

Ugh!  Mathematica doesn't have enough common sense!  We can combine a lot of terms
	here.  Below we force it to do so.  Execute the next statement, and again wait for a bit.

Simplify[PowerExpand[DSolve[y''[t]+g y'[t]+(w^2) y[t]==Sin[w0 t],y[t],t]]]

If you look closely, you will see our two friends, C[1] and C[2].  They again multiply exponential functions.   However, you will see a third term, too!  This is called the
inhomogenous part.  It is present because we have a periodic driving force.

 As before, we'd like to define this as a function, and then plot it for a variety of	different cases.  We'll call the resulting function s[t, g, w].  This is done in the next 	Mathematica line below.  Execute it now, and read the next paragraph while you wait            for it to do the calculation.   The semicolon at the end  of the line will stop it from printing 	out the solution to the screen, which would wasteat lof of room. You'll have to wait 	 a bit.

      While you wait, you should take a look at the expression closely.  See anything      different?  We have set w = 1.  We could have set it equal to anything, but one  is a nice number.  This is the natural frequency of the oscillator.   We have left the frequency of the driving force (w0) as a variable, so that we can drive the system faster or slower than it would naturally oscillate.   You should also note that we are looking for the solution where the oscillator is initially at rest, at the origin. The driving force starts it  moving.

s[t_,g_,w0_]=y[t]/.DSolve[{y''[t]+g y'[t]+y[t]==Sin[w0 t],y[0]==0,y'[0]==0},y,t][[1]];

This next command will tell it to simply this result as best it can.  Again, execute it and wait for a bit.

s[t_,g_,w0_]=Simplify[PowerExpand[s[t,g,w0]]]

Ok!  Now we are ready to start.  Let's  choose g=10 (very high damping) and w0=1. We'll plot it for a range of t from 0 to 20pi (that's ten oscillations).  We actually don't start exactly at t=0, because Mathematica gets a little confused there.  We start at t=0.1    You'll also note this funny Evaluate[expression].  It actually isn't required; it just makes things go faster.  It tells Mathematica to stick in the values for g and w0 in the above expression first, and generate a new function.  It then uses that function for each value of t when it plots the curve.  If we didn't do this, it would use the big ugly expression above for every point it plots, recalculating it at every single darn   point.  That's a waste of time!

Plot[Evaluate[s[t,10,1]],{t,.1,20 π}]

Note that it takes a while to settle down to a steady oscillation.  The initial bit is called the transient response.  The steady oscillations are called the steady state behavior.

(4)  Copy the above plotting command ard replot the function for a variety of
	  different values of w0.   Does the steady state behavior oscillate at w, the
	  natural frequency of the oscillator (which we have set equal to one), does it
	  oscillate at w0, the frequency you select, or at some other frequency?   Be 
	  sure you plot enough of the function to see its steady state behavior.

First we'd like to study the case when w0=w.  This is called resonance. The oscillator and the driving force are at the same frequency.  In our case this means choosing w0=1.  Below is a sample plot.

Plot[Evaluate[s[t,.5,1]],{t,0.1,20 π}]

(5)  Plot this function for g=0.20, g=0.10, and g=.05.  Describe the effect in
	words.  You should describe how the amplitude of the oscillations varies
	with g, and how the transient response varies with g.

     	     Next, fix g=0.20, and vary w0.  Try values of w0 slightly greater and
	slightly less than w (which we have set equal to  one).   How does the 
	amplitude of the steady state response depend upon w0?

Ok, let's put this altogether in one 3D plot.  Below is a command that looks at the solution  for t between 20π and 24π.  This should be far enough in time that the transient behavior has settled down, and we are only looking at the steady state behavior.  In this graph we are fixing the driving frequency so that we are on resonance, and changing the dissipation  constant g.  The "PlotPoints" option specifies how many points are plotted in the g and t intervals in the range of the plot.  The "BoxRatios" option  allows us to stretch the box so that it is longer than it is wide or tall.   Finally, we can specify the location of our viewpoint.  (This is done in  an abstract coordinate system, not measured in the units of g and t).  You
may wish to play with the options.

Plot3D[Evaluate[s[t,g,1]],{t,20 π,24 π},{g,0.1`,1},PlotPoints->{30,15},BoxRatios->{2,1,1},ViewPoint->{1,2,2}]

(6)  Explain in words why the above graph has the shape it does.  Why is it
	flat where it is flat, where is it more oscillatory and why?

This next graph fixes the dissipation at a small value, g=0.005, and then plots t he steady state behavior for different values of the driving freqency.

Plot3D[Evaluate[s[t,0.005`,w0]],{t,20 π,30 π},{w0,0.7`,1.3`},PlotPoints->{40,100},BoxRatios->{4,2,1},ViewPoint->{-2.692`,-0.795`,1.889`},PlotRange->All]

(7)  Explain the shape of this plot.  In particular you must explain why the
	ridge down the center is so high, and it's physical meaning.  How could you
	make it rise even sharper?  How can you flatten it?