Resonance with a Damped-Driven Mass-Spring System in Python

Image for post
Image for post
Photo: Rhett Allain

Physicists like to say that pretty much everything can be modeled as just a mass connected to a spring. I mean, it’s not completely untrue. The mass on a spring appears in lots of different problems. But what about an oscillating mass WITH some type of frictional force? What about a mass-spring with a driving force? This is what we are going to build — in python. It’s going to be awesome.

Mass with a Spring

If you want to build a complicated model, the best option is to start with something simple and make sure that works correctly and then add the complicated parts. That’s what I’m going to do here. So, let’s consider a mass connected to a spring without friction and without a driving force. Just a plain simple harmonic oscillator.

Of course, this all starts with the force due to a spring. Here is a mass connected to a spring. I’ve connected the other end of the spring to some “anchor” point (this will be important later).

Here are the key details for the force due to a spring.

  • The force is in the same direction as the spring. It can either push or pull — which ever way is back towards it’s equilibrium length (the unstretched length).
  • The magnitude of this spring force is proportional to the change in length of the spring.
  • This proportionality constant for the force is called the spring constant (k) and has units of Newtons per meter.

With that, we have the following expression for the force from a spring.

In this expression, L is a vector from the beginning to the end of the spring and L0 is the unstretched length of the spring. Oh, the L with the hat over it (L-hat) is a unit vector in the direction of the length of the spring. This unit vector is needed so that the whole thing is still a vector. If you have the spring force as a vector (and not just a scalar), you can do some seriously powerful things.

If you know the force on an object, you can determine its motion. In this case, we can use the momentum principle.

Where the momentum is defined as:

Of course there is a problem here — the force from the spring is not constant so that it’s not so straight forward to find the new momentum. One way to fix this is to break the motion into many very small time steps. During each of these short time steps, we can assume the force is constant instead of changing. With this assumed-constant force, we can then use the momentum principle and find the momentum (and position) of the mass at the end of the time interval. Then you just do this same thing again for the next time interval.

This process is called a numerical calculation. Here is the recipe. First, use the position of the mass (represented by the vector r) to calculate the spring force. Note: I also have to calculate the length-vector for the spring.

Here I am using ra as the vector location of the “anchor” point so that the length of the spring is from the anchor to the mass (position r). This step might not make much sense now — but just you wait and you will see why it’s useful.

With this force, we can use the momentum principle to find the momentum at the end of the time interval.

Just to be clear — with the momentum at the start of the interval (p1), you can find the momentum at the end of the interval (p2). But next, the position of the mass needs to be updated. We can find that using the definition of average velocity.

OK, there’s a trick here. I started with the definition of the average velocity, but in the second line I use p2 (the momentum at the end of the time interval) instead of p-average. Yes, this is wrong. However, if the time interval is small, this is only “sort of wrong”. Trust me, this will work out fine.

Of course, no one wants to do all of these calculations by had — so I’m going to make a computer do it. Instead of setting up every part of the code, I’m just going to show it to you and then explain the important parts. It’s a little long, so here is part 1.

Yup, there’s lot’s of cool stuff here — it might look bad, but it’s not that crazy. Here are my comments (using the line numbers above). Oh, here’s the full code.

  • 1: This line is sort of like importing a module in python. This is Glowscript Vpython, an online version of python that has this built in module. Anyway, this line just loads all that stuff. Don’t delete it.
  • 2–5: These lines just set up the graph that I’m going to make. Here’s a different post that has everything you need to know about graphs.
  • 9–11: The most awesome thing about Glowscript are the 3D objects. sphere() is one of these pre-built objects. For the parameters in the “sphere” object, you need to give it a vector location of the center — this is the pos. You also should give it a radius. It’s that simple. Notice that the anchor is at an x-value of -L0/2 and the mass is at x=L0/2+0.02. Since the unstretch length of the spring is L0, this means the spring is going to start off stretched 0.02 meters.
  • 12–17: Mostly just setting up constants and initial conditions. Notice that for the mass of the “mass”, I use mass.m = 0.02. This assigns the property “m” to the object mass. It’s just a nice way to keep up with all your variables and stuff. Notice line 17: ra0=anchor.pos. This creates the variable ra0 and sets its value to the vector position of the anchor. Why? I’m going to need this later.
  • 19: This is another object, the helix. It’s used to make the 3D representation of the spring. Don’t worry about this too much.
  • 21–22: Since I’m breaking stuff into small time steps, I need to know the size of the steps. That’s the dt. Also, time starts at t = 0.
  • 23: Since I’m going to use it inside the loop (to come below), I need an initial value for the vector length of the spring.

Now for part 2.

  • 26: This is the start of the loop. Everything indented after this while statement is part of the loop. This tells the computer to do all that stuff as long as the time variable (t) is less that 5. Oh, I guess I should add that the program is a numerical calculation. It uses numbers. It doesn’t use units. You (the programmer) have to keep track of the units.
  • 27: The rate(1000) tells the program to not do MORE than 1000 loops per second. Why does this matter? Well, it’s important for the 3D visualization. Since the time step is 0.001 seconds, a rate of 1000 would be in “real time”.
  • 28–30: Calculate the spring force and set it equal to the total force (I will add other forces later). Note that mag(L) returns the magnitude of the vector L and norm(L) is a unit vector in the direction of L.
  • 31: Update the momentum (from the momentum principle above). Note that the = is not an algebraic equal sign. Instead it’s a “make equal to” sign. It says to take the old momentum, add F*dt and then make that the new momentum. That means we don’t have to keep track of p1, p2, p3…we just need one p variable that keeps changing.
  • 32: Update the position. Note that this updates the vector value of the location of the mass — but it also actually moves the 3D representation of the mass.
  • 37–40: Calculate the energy. One thing I can check is to see if the total energy is constant. I’ll show that below.
  • 41–43: These plot data points for the kinetic energy, the spring potential energy, and the total energy.

OK, here’s what this looks like when you run the program.

That’s pretty, isn’t it? But is it legit? Here is the plot of energy for this system.

The green curve (total energy) changes slightly — but it looks good enough for now. I originally had a time step of 0.01 seconds, but the energy wasn’t constant enough.

Mass with a Spring AND Viscous Drag

Now suppose the mass is in something like oil. In this case, there will be an extra force pushing on the mass that is in the opposite direction as the velocity of the mass (and the opposite direction as the momentum vector).

In the case of something like oil, we would call this viscous drag. The magnitude of this drag force (labeled Fv above) would be proportional the the velocity of the mass. Since I’m using the momentum for the other calculations, I can write this in terms of momentum.

In this expression, c is some type of coefficient that depends on the size and shape of the object and the thickness of the oil. Of course I’m going to need a value for this drag coefficient. But it’s already in the code. I’m going to start with c = 0.01 Ns/m.

Now I just need to add a line into the code. Here is the modified version.

The only new parts are line 29 and 30. Line 29 is a calculation of the drag force and then I add this into the Fnet calculation. The rest of the code is identical. When you run it, this is what it looks like.

The graphs look better. Here is a plot of the energy as a function of time.

Noice. Of course the total energy isn’t constant since there is work done on the system by the viscous force. If you want, you can change the drag coefficient to something much higher, you get something like this.

But you should go ahead and put your own values for the drag coefficient. It’s fun. Here is a link to the actual code. Notice that we did a lot of work for the first part without viscous drag such that it was fairly trivial to add this force and get the code to work.

Mass with Viscous Force and Driving Force

See how this works? We build a model and get parts of it to work before adding another part. Now we just need to add this driving force. Here’s how it’s going to work. I’m going to move the “anchor” point back and forth. This means there won’t be a separate driving force on the mass but instead will be caused by the compression and stretching of the spring due to the movement of the anchor.

Remember how I calculated the vector length (L) of the spring? Yup, that was based on the position of the anchor. So, all I need to do is to move the anchor back and forth.

Here is a function to calculate the location of the anchor in each time interval.

This says that the x-position of the anchor will oscillate back and forth using this sine function. The amplitude of oscillation is D and the angular frequency is ωd. I’m using the d-subscript to indicate this is the angular frequency of the “driving” force. Oh, I have to then add this to the center position of the anchor.

For the actual value of the driving frequency, I want it to be somewhere around the natural frequency of the oscillating mass. For a mass on a spring by itself, it should have the following angular frequency.

We can easily change this driving frequency, but I’ll start with an initial value that is 0.9 times this natural frequency.

OK, let’s fix the code.

I only need to add line 26. This is exactly the same equation as above for the motion of the anchor. Everything else should work out. Let’s do it. I’m not going to include an animation (because it’s not very interesting) — but here is the code if you want to play with it.

Check out the energy plot.

There’s so much more to do. I can’t do it all. That means it’s homework for you.

  • What happens when you change the driving frequency?
  • What is critical dampening? Can you get it to do that?
  • Can you make a plot of amplitude as a function of driving frequency?

Physics faculty, science blogger of all things geek. Technical Consultant for CBS MacGyver and MythBusters. Former WIRED blogger.

Get the Medium app

A button that says 'Download on the App Store', and if clicked it will lead you to the iOS App store
A button that says 'Get it on, Google Play', and if clicked it will lead you to the Google Play store