The Fourth Season:  Spring in Four Acts

Everything you didn't realize you wanted to know about springs but never got around to asking

Marcus McGuff, A.Mus.D.
Austin Community College

Introduction

One of the simplest second order differential equations can be used to model the oscillation of a spring with a mass attached.  I will be the first to admit that this doesn't really sound all that interesting.  I have good news:  It is, and more importantly, there are several important characteristics of this system that show up in much more complicated systems.  We are going to take a look at some of those characteristics and I will point out some of their wider applications.  (Like many other things in math, we can start to understand things with a simpler system, where things don't get all that complicated.  From there, you can work your way up to reality...)

So, to start off, we need to derive the equation of motion of a mass suspended from a spring; we will ignore the mass of the spring itself for simplicity.  As usual, when trying to find the equation of motion of an object, we need to find the total forces acting upon it, so we can use Newton's second law of motion (the ever popular F=m a).  In this case, there are two forces:  the force exerted by the spring and the force exerted by gravity due to the weight of the mass.  The force exerted by the spring can be modeled by Hooke's Law F=-k y, where y is the distance the spring has been stretched beyond equilibrium (I have arbitrarily decided to consider +y to be upwards in the figure below), where k is called the spring constant, which measures the "stiffness" of the spring; note that it is negative so that it acts in opposition to the displacement.

					
					
Now, it turns out that we can ignore the force due to gravity here:  Since the force due to gravity remains constant, that means that it will cause the spring to stretch an extra amount even when at rest; that extra stretch remains the same, no matter what the mass is doing (just solve (Subscript[F, gravity])=-m g=-k y=(Subscript[F, spring]) for y).  So, the basic effect of this is to simply make our "equilibrium length" a bit longer.  (Okay, this is only strictly true for a perfectly linear spring.  If you really feel the urge, you can enter in a gravity term to your "spring force" in the widget...)  So, our equilibrium point will be measured from this gravitationally stretched spring.

Now, if we put this all together, we get m y''=F=-k y or:   m y''+k y=0.

In the lab that follows, many of the questions can actually be answered without knowing how to use Mathematica beyond using the spring simulation widget I included with the lab.  Some of the things I ask you to compute using Mathematica can be done fairly easily by hand (certainly all of Acts I and II); if you are not required to know how to know Mathematica in your course, you should be able to answer most of the questions in the lab without using anything beyond the spring widget.  If you are required to know how to use Mathematica some, then you should use it to compute the things I ask for.  In either case, when you describe the impact of some parameter on the system, you need to give me examples; you can do that by copying the entire widget or just clicking on the particular graph you want to copy and then copy/paste that.  (Don't get carried away, but I do want some examples showing what behavior you are seeing.)

Act I:  Simple Harmonic Oscillation

When the forces mentioned above are the only ones considered (i.e., no damping or driving force), the motion of the system is an example of what is called "simple harmonic oscillation".  This equation is quite simple to solve by hand, but I would like you to start actually doing some work within Mathematica, so:

Use DSolve to solve the equation for this simple harmonic oscillator, with initial conditions y(0)=y0 and y'(0)=v0, in terms of m and k.

Now, if you recall your trig identities (or even if you don't), we can convert a sum of sin and cos into a single sin (or cos) with a new amplitude and a phase shift, as long as the angles are the same.  Again, you can do this by hand without too much trouble, but I'll save you the trouble by giving you a Mathematica function that will convert this for you (you can thank me later - this is actually kind of a cool trick):

combineSinCos[stuff_]:=Apart[stuff]/.(a_ Cos[θ_]+b_ Sin[θ_])->((Sqrt[(a^2)+(b^2)])Sin[θ+ArcTan[(a/b)]])

So, you can use this to display the formula:

combineSinCos[A Sin[x]+B Cos[x]]

(Notice that the order of the sin and cos doesn't matter here...  I did cheat a little here, since I assumed both a and b are positive; that works fine for our purposes, but this needs a little more work for generic use.)

Use combineSinCos to convert the solution you found above into a single sin function.  What is the frequency (the angular frequency) of this solution?  This is sometimes called the "natural frequency" of the system.

So, note that all solutions to this system are given by a sin curve, with an amplitude, period, and phase angle determined by the various constants involved.

Now, take a look at the "spring simulation" widget at the end of this lab.  Start it running and play around with it for a little bit.  Note the graphs of y(t) and y'(t) on the left, but also pay attention to the phase graph on the lower right (y vs. y'); it can be especially helpful in helping you see the relationship between the position and velocity of the spring (and you can also set your initial conditions by clicking on it).

Use the spring simulation to answer the following questions (you may want to run your time control out to the right so you can see how the graphs change; pay attention to the scaling of the graphs and start most of these with both initial conditions non-zero):

What is the impact on the motion/velocity of the system if the mass is increased/decreased?

What is the impact on the motion/velocity of the system if the spring constant is increased/decreased (increasing this number is equivalent to a "stiffer" spring)?

What is the impact on the motion/velocity of the system if the initial conditions are changed?

Now, go back and look at the solution you found above and explain mathematically why each of these 3 things behave as they do.  Also, try to give a physical explanation as well.  (So, if the spring gets stiffer, it makes sense that the spring will...)

What is the relationship between the phase of y(t) and y'(t)?

Act II:  Reality lends a hand - Damping

Of course, in reality, there will be some energy loss over time, so the spring can't remain a perfect simple harmonic oscillator forever.  There are various ways this loss can be modeled, but we are going to consider the case where the spring is "damped" with a force proportional to its velocity.  That means the forces acting on the spring become F=-k y-b y' (the damping acts in opposition to the velocity), so the equation of motion becomes:  m y''+b y'+k y=0.

Solve the new equation with damping included.  (Hint:  Don't bother setting the initial conditions; just get the general solution, it's much easier to understand).

Why is (b^2)-4k m important?  Carefully examine the behavior of the spring for the different values of (b^2)-4k m, explain how they are different, and give examples of each case (there are 3 cases).  The system is said to be "critically damped" if (b^2)-4k m=0.

Sometimes, it can be useful to rewrite the differential equation by renaming a few of the parameters.  For example, the natural frequency of the undamped system is often called (Subscript[ω, 0])=(Sqrt[(k/m)]).  The damping ratio of the system is defined to be the ratio of the actual damping coefficient b to the damping coefficient that would yield critical damping (Subscript[b, crit]).  So, for our system (Subscript[b, crit])=2(Sqrt[k m]) (just solve (b^2)-4k m=0 for b).  So the damping ratio is given by ζ=(b/(Subscript[b, crit]))=(b/(2(Sqrt[k m]))) (note that this means b=2ζ(Sqrt[k m]); but using our formula for (Subscript[ω, 0]), we also know that (Sqrt[k m])=m (Subscript[ω, 0]) so we get b=2ζ m (Subscript[ω, 0])).  So with a little substitution/algebra, we can take the exponents from our previous solution and rewrite them in terms of these new parameters.  First, we substitute in for b:

step1={(((-b-(Sqrt[(b^2)-4 k m])))/(2 m)),(((-b+(Sqrt[(b^2)-4 k m])))/(2 m))}/.b->2ζ m ω0

Remember that /.b-> 7a means replace all the b's with 7a.  (Very big, very important warning: Don't use (Subscript[ω, 0]) within a Mathematica calculation.  It gets confused between (Subscript[ω, 0]) and ω; it is much safer to use ω0 instead.)  Notice, we still have k under the radical; we can solve (Sqrt[k m])=m (Subscript[ω, 0]) for k:

step2=step1/.k ->m ((ω0)^2)

Now, clearly there is some factoring and canceling that should happen right about now.  This is a good time to try the Simplify command; however, to really get it to factor things out from under the radicals, you have to tell Mathematica that we know that all the parameters are positive (this is a useful trick...):

ωNew=Simplify[step2,Assumptions-> m>0&&ω0>0&&ζ>0]

So, what's the point?  If it seems like we have just traded one complicated set of formulas for another, that's true.  However, by using the specific substitutions given here, we have discovered something important about the solution to our system:  These exponents are just multiples of the natural frequency of the system, with the coefficient determined by the damping ratio.  We have three cases:

If ζ=1, the system is said to be "critically damped" (by definition).  In this case, the system returns towards the equilibrium position with no more than one "overshoot", i.e., there is no oscillation.  This marks the dividing line between over-damped and under-damped.  So, in a sense, critical damping is the minimum amount of damping/energy loss needed to keep the system from oscillating.  Critical damping is often desirable in practical engineering applications, such as the design of shock absorbers and automatic door closing mechanisms.

Notice that the solution Mathematica gave you to the original differential equation isn't actually complete in this case.  (Why is this?)  Since ζ=1 means b=2(Sqrt[k m]), solve the original equation again to find the complete solution for this case.

If ζ>1, the system is said to be "over-damped".  The system behaves very similarly to critical damping, i.e., there is no oscillation and the system heads back towards equilibrium.

If ζ<1, the system is said to be "under-damped".  In this case, the exponential powers become complex, so there is an oscillatory component to the motion.  This is mathematically the most interesting case, but is generally considered bad in many physical situations (think shock absorbers...).

Consider that the solution is of the form (e^((-ζ±(Sqrt[-1+(ζ^2)])) (Subscript[ω, 0]) t)), find a formula for the exponential and sinusoidal parts of the solution.  (You do not have to combine the sin and cos parts of the answer.)  Look at several examples of this (different initial conditions) in the widget; notice that the graph of y(t) is basically a sin (or cos) curve with an amplitude modulated by an exponential.  Take a look at the formula you just found and tell me what determines the rate at which the oscillations die out.  How does the frequency of this damped oscillation compare to that of the "natural frequency" of the spring?  (Again, check out that formula you just found...)

Act III:  Driven to succeed  - The forced spring

Another important case to consider is when the spring is "driven" by some external force.  This has quite a few direct applications, but some of the behaviors we look at here show up in many other non-spring-related fields (acoustics and electronics, for example).  Let's start by going back to a simple harmonic oscillator (no damping), but with some external driving force applied to it; for our model, we are injecting that force by moving a piston on which the spring rests:
                                                                      

Okay, I will admit that this would have been simpler if I had kept the spring stationary and driven it by a string attached to the mass or something.  This way just seemed like the obvious way to drive the spring, however.  Sorry.  (Well, not really...)  Why is this more complicated?  Well, y(t) is measured from the equilibrium length of the spring; that means as the piston moves, so does the origin.  That's why I added the "Eq" arrow; notice that it moves with the piston - y(t) is measured relative to that point.  Notice that I do give you an extra "Position: spring+piston" graph at the bottom; this is simply a plot of y(t)+P(t), where P(t) is the position of the piston at time t.  Note that P(t)≠F(t); F(t) is a force, so it is proportional to P''(t).  You might want to manually move the time slider back and forth for a bit and look at how all the graphs relate to the spring until you get a feel for it.  (Notice that I graph the forcing function F(t) on the same graph as y(t).  This is a bit bogus, since they are totally different units, but I think it makes it much easier to see how y(t) is responding to the forcing function this way.)

We will consider a sinusoidal forcing function F(t)=A sin(ω t) (there is a second one with a phase shift you can experiment around with as well), since that gives some really interesting behavior.  The equation for this system is pretty simple:  m y''+k y=A sin(ω t).

Solve the new equation with this forcing function.  (Hint:  Don't bother setting the initial conditions; just get the general solution, it's much easier to understand.  Also, you might want to run it through FullSimplify, or just use some basic trig identities.)  Notice that there are three terms here; two of them have the same frequency (the terms with the C[1] and C[2]) and are the same as what we got for the case of simple harmonic oscillation.  These could be combined into a single sin function with frequency (Subscript[ω, 0]) (the amplitude and phase angle won't be quite the same, however).  The third term is the forcing term with its amplitude modified by k, m, and ω.  (Notice that at least part of it comes from integrating F(t) to get P(t).)

Use the widget to investigate the behavior of the system in the following situations:

If the driving frequency ω is much higher than the natural frequency of the system (Subscript[ω, 0]), how will the system behave?  Start with a small amplitude A and then gradually increase it.  What impact does it have on the overall frequency and amplitude of the system?  Is the overall period of the solution closer to ω or (Subscript[ω, 0])?  (If you "uncheck" Driving force, you can flip back and forth between the "natural" system and the driven one with the same parameters.)  Try to explain physically why this happens.

If the driving frequency ω is much smaller than the natural frequency of the system (Subscript[ω, 0]), how will the system behave?  Start with a small amplitude A and then gradually increase it.  What impact does it have on the overall frequency and amplitude of the system?  Is the overall period of the solution closer to ω or (Subscript[ω, 0])?  (You will probably need to look at a fairly long time interval to answer this one.)  Try to explain physically why this happens.

If the driving frequency ω is exactly the same as the natural frequency of the system (Subscript[ω, 0]), how will the system behave?  Start with a small amplitude A and then gradually increase it.  What impact does it have on the overall frequency and amplitude of the system?  (You can set an exact value for a parameter if you click on the little + symbol to the right of the "slider".)   This is an example of what is called "resonance"; it can be a very bad think bridge oscillations) or a very good thing (think signal amplifier).  Try to explain physically why this happens.  How can this occur mathematically?  (Hint:  Mathematica is not always your friend...)  You do not have to find the equation for the solution in this situation; just explain what's wrong with the other solution here.  (The reason the solution behaves as it does is that, since the driving function has the same frequency as the natural frequency, it is itself a solution of the equation; this means you need a t sin(ω t) term in your solution.)

Now, change the driving frequency by just a small amount (say 0.1 or 0.2) from being exactly equal to the natural frequency; what happens?  Start with a small amplitude A and then gradually increase it.  (You may need to look at a longer time interval to really see what is happening.)  This is an example of a very odd behavior called "beats" that shows up when two sinusoids of very close frequency are combined together.

What impact does making the mass larger/smaller have on the system?  What impact does making the spring constant larger/smaller have?  (Hint:  Don't spend forever playing with the widget on this one; just think about the math.

Act IV:  Happily ever after - Damping and forcing get together

Of course, there will usually be some damping component active in the case of a driven spring as well.  So we need to combine Acts II and III together to see how the damping and driving forces interact.  The equation is just a combination of the two:  m y''+b y'+k y==A sin(ω t) (again, we'll stick to a sinusoidal driving force).

Solve the new equation.  (Hint:  Don't bother setting the initial conditions; just get the general solution, it's much easier to understand.  Also, you might want to run it through Simplify.)  Notice that there are three terms here; two of them are the same as in the undriven damped case (the C[1] and C[2] aren't necessarily the same).  These have the same solution types (critical/under/over damping) as before. The third term (actually two terms) has both a sin and a cos function with the same period as your driving function; you can use combineSinCos to rewrite that term as a single sin term with the same period as the driving force, but a phase shift.

If you use the substitutions we used earlier (for (Subscript[ω, 0]) and ζ), you can find the phase angle to be ArcTan[((2 ζ ω ω0)/((ω^2)-((ω0)^2)))].  What does this tell you about the impact of ω's relation to (Subscript[ω, 0])'s will be?  (In other words, what does it mean if ω>(Subscript[ω, 0]), etc.?)  Check your mathematical conclusion against the simulation.

Also, notice that this means the system will always have an oscillatory component.  If the system is critically damped or over-damped, what will the long-term behavior of the system be?  Notice, it will take a little while to settle down to this behavior; the early behavior is sometimes called the "transient" behavior of the system; this can be quite complicated, but will eventually damp out, leaving a "steady-state" solution.

If the system is under-damped, does it also have a similar steady-state solution?  Explain/give examples.

What is the effect of damping on resonance?  First, make a guess based on some sort of physical/mathematical reasoning and explain why.  Now, use the spring widget to see what happens.  (You can set up the resonance with damping turned off, then switch damping on and vary it to see what happens; be sure to examine both side of critical damping.  You'll want to look at some longer time periods on this one.)  If this didn't match your guess, how might you modify your reasoning to account for the observed virtual facts?

This is not necessarily a question about forced/damped oscillation, but rather a more general question:  Explain the usefulness of the phase plane (y vs y') in figuring out what is going on with the system?  Please stay away from vague pretend answers ("It shows me the relationship between y and y'."  This is of course indisputably true, but almost totally useless...).  Pick a few examples and explain specifically what it is about the behavior that the phase plane is particularly helpful in showing.

Encore:  The spring reexamined

So far, we have assumed that our spring is perfectly linear.  However, not all springs are linear and even those that are sometimes change over time (they can age and lose elasticity or become more brittle).  One way model this "aging" spring is to assume that the spring constant gradually changes from one constant (Subscript[k, 0])to another one k.  (If (Subscript[k, 0])>k, the spring is becoming less elastic, while if (Subscript[k, 0])<k it is becoming "stiffer".)  We could try a linear transition (Subscript[k, 0])(1-t)+k t, which would transition between t=0 and t=1.  We could scale this so it happens over a specific time period, but we would need to be careful not to model beyond the end time period.  If we replace t with (e^(-t)), then we don't have to worry about our time domain.  If we throw in a decay rate, we get something like (Subscript[k, 0])(e^(-ϵ t))+(1-(e^(-ϵ t)))k (included as a choice for the spring in the widget).  For a simple harmonic oscillator, this gives us:  m y''+((Subscript[k, 0])(e^(-ϵ t))+(1-(e^(-ϵ t)))k)y=0.

This equation can be solved, but you will not enjoy it...  Take a look at this with the spring widget.  To make it easier to see, set ϵ very small so the transitions happens much more gradually and be sure k and (Subscript[k, 0]) are far apart.  What happens if (Subscript[k, 0])>k?  What happens if (Subscript[k, 0])<k?  Without knowing anything about the form of the solution, explain why the things you see here happen.

How would you expect this to effect the damped, driven, and driven/damped cases?

Another change that can be made to the spring is to assume that if you stretch/compress it too far, it will lose elasticity.  So, this isn't a matter of the spring aging, but tries to reflect what happens when you pull it too far, its return force won't be as strong as it should be.  Again, this could be modeled in many ways, but one method might be to provide a nonlinear term to reduce the spring return force, like F=-(k y-ϵ (y^3)).  For this to work k y-ϵ (y^3)≥0 for positive y and k y-ϵ (y^3)≤0 for negative y.  This has intercepts at y=0, y=±(Sqrt[(k/ϵ)])and (think about your algebra) it is positive 0≤y≤(Sqrt[(k/ϵ)]) and negative -(Sqrt[(k/ϵ)])≤y≤0.  So, this only makes sense if  -(Sqrt[(k/ϵ)])≤y≤ (Sqrt[(k/ϵ)]).   (This is not a nice thing...)

Take a look at this with the spring widget.  To make it easier to see, start with ϵ very small and k very large. This will be very tricky, since as you change either of these, the amplitude will also be changing, and once the amplitude gets bigger than (Sqrt[(k/ϵ)]), everything will go up in smoke.  Don't freak out when the widget goes crazy here; just move the k or ϵ back a bit until the graph reappears.  In order to really see how this is different, you will need to push ϵ and k as close to each other as you can (i.e., maximize the non-linearity).  How does the non-linearity seem to effect the behavior of the system?  Look at all three graphs carefully to really understand this.

The spring simulation widget

If you want to enter your own driving, damping, or spring force function and you are running a full version of Mathematica, you can click on the little triangle next to each of those to open up a box you can type directly into (either using standard mathematical form or Mathematica notation).  Be warned that if you choose something too complex, you may run into some errors (so feel free to experiment, but don't complain to me if something bizarre happens).  You can set initial conditions using either the y(0)/v(0) sliders or by clicking your mouse on the y/y' phase plane.  If you want to "zoom in" to look more carefully at the transient response, click on "shorter time"; if you want to "zoom out" to examine more of the long-term behavior, click on "longer time".  The driving force must be only a function of time, t; any other "variables" that are in that formula will be treated as parameters.