Using VPython simulations to model chaotic motion and investigate what defines a chaotic system.
Any physics student has examined the lowly pendulum, probably more than they wanted to. Along with the mass on a spring, pendulums are the quintessential example of simple harmonic motion. The keyword there is simple. They, as you may know, are incredibly easy to model mathematically. The period of the pendulum can be represented by the formula below, where g is the force of gravity and L is the length of the pendulum.
The angular displacement is somewhat more interesting but still not too deep. The equation itself is more complex, but in reality it is just a cosine wave. The equation below represents the angular displacement of a single pendulum, where θ is the angular displacement, A is the maximum angular displacement (in most cases this is just the release angle), ω is the angular frequency, and t is time.
With that out of the way, let’s get to the more interesting, albeit effectively impossible (you’ll see), mathematics and physics.
Odds are, if you’ve worked with a regular pendulum, you’re at least somewhat familiar with the double pendulum, which is reminiscent of the infamous three-body problem. The three-body problem represents the idea that an analytical solution to the motion of three bodies in a system where all bodies’ forces affect each other is impossible to solve, since the forces of each body interact with the others simultaneously. In short, the double pendulum is impossible to solve analytically. In terms of the pendulum equations stated above, the one that you really can’t solve is the equation for θ, or the angular displacement of both pendulums. You also cannot solve for the period, but that’s mostly because a double pendulum doesn’t exactly have a regular period. Below is a visual representation of the double pendulum system we’ve discussed, where the pivot point is anchored to the wall, the angles are marked with θ1 and θ2, and the hanging masses are marked m1 and m2.
To reiterate– in our current physics, we cannot predict the motion of a double pendulum analytically. But why? Well, the motion is not simply random, but rather chaotic. How do we prove that chaotic and random aren’t the same? That is the question we will seek to answer later on.
You’ve likely seen some YouTube title about chaos or Chaos Theory, and I’m willing to bet the thumbnail image looked something like this:
That image, along with many other butterfly wing-esque drawings, are most commonly associated with chaos theory. Something not commonly thought of as part of chaos theory is the double, triple, and any other pendulum with more than one arm. Let’s explore some of the properties of chaotic motion and try to learn more about it as best we can, with help from a computer program that we will dive into later.
One property we can investigate is the effect that making a small change in starting conditions has on the end result. Will the change be negligible in the end, or will it yield an answer entirely outcome from the original?Below is a computer simulation of a double pendulum at t=50, where the starting angles were 80° and 0° respectively for θ1 and θ2. (Angles are the same as shown in the diagram above)
Now, let’s run the exact same simulation, the image is of t=50, all the masses and arm lengths are completely identical. However, the angles are 79° and 0° for θ1 and θ2. There’s a 1° change in the starting value of θ1, let’s see how it affects the outcome.
Unsurprisingly (or surprisingly, I don’t know what you expected), with a change of just 1° the system has dramatically changed to the point where the final position of the system and the line traced by it don’t resemble the first image at all. That’s independent of any outside variables, as in our simulation there is no friction, air resistance, or any other hindrances. This example shows the idea that when a system exhibits chaotic motion, as in a double pendulum, even the most minute change in conditions can cause an enormous change in the result.
So, we’ve nailed down one interesting aspect of chaotic motion– a tiny change in starting condition can cause an enormous change in result. Let’s try to identify a few more. What would happen if we were to run the same double pendulum experiment again? Experimentally, this is all but impossible; it would involve having an infinite degree of precision when measuring things like gravitational field, mass, air resistance, and every single other variable in the experiment, and then using your magical infinite precision data to set up the experiment again, identically. We, however, have a nice simulation to use which makes that challenge disappear.
Let’s run the double pendulum simulation similar to before, but with some variables changed so the image isn’t the same: t=50, θ1=80°, θ2=15°, and the physical properties of the pendulum are the same (length, mass). Here is the result:
Spoiler alert, it was the exact same the second time… And the third, fourth fifth, and sixth time. Just to prove a point, let’s also try this exercise with a triple pendulum. We’ll use the same starting conditions, although now we have a third angle, θ3, so we’ll use θ1=80°, θ2=15°, θ3=0°. A triple pendulum is much more computationally intensive, though, so the image will be taken at t=10.
Although it is even more complex than the double pendulum simulation, the fact that the triple pendulum exhibits the same tendency to repeat its pattern under identical starting conditions is important. That idea of an identical outcome under identical conditions proves that chaotic motion and chaotic systems are deterministic, thus not random. Again, since this is so hard to prove experimentally or use in practice it doesn’t affect our understanding of the real world very much. We do, however, now know that chaotic systems are not random and can be predicted with the right knowledge. Additionally, since we proved that a small change can have a large effect, we know that an incredibly high degree of precision is needed to effectively predict the outcome of a chaotic system.
We’ve now identified two significant properties of chaotic motion: the first is that a tiny change can have huge implications in the outcome of the system, and the second is that chaotic motion is not random but rather it is deterministic. Let’s try to find one more– as we discussed at the very beginning of the article, the variable you really can’t solve for in a double pendulum is the angular displacement. However, we also said that solving for period was futile because a double pendulum does not repeat– let’s test that theory.
First, let’s run the double pendulum experiment with the following conditions: θ1=30°, θ2=30°, with the same masses and lengths as before. This time, let’s use VPython’s graphing capabilities to record a chart of the height of the bottom mass as a function of time:
If you ask me, that looks pretty periodic. So, does this mean that the system described above is not a chaotic system? And if the motion is not chaotic, does that mean we could reach an analytical solution to the motion of the double pendulum when it has less than some threshold energy needed to become chaotic? I’ll leave that as an exercise to the reader. To do that though, you might need the help of the program I wrote and used in this article, so let’s walk through it. This program is written in VPython using GlowScript, so be sure that you’re running it on the WebVPython system, and not just any python editor.
First, let’s define the physical constants we’ll use in the program. These include the masses, lengths of arms, and force of gravity. Additionally, we’ll define the starting conditions of the system, where theta1 and theta2 are the starting angles, and theta1dot and theta2dot are the starting angular velocities of the arms. Finally we’ll also set the starting time to 0, and define dt, which is the amount of time we’ll step forwards in each iteration of our loop later.
#lengths of the strings
L1 = 1.5 #top string
L2 = 1 #middle string#masses
M1 = 2 #top moving ball
M2 = 1 #middle moving ball
#g
g = 9.8
#starting angles and velocities
theta1=30*pi/180 #release angle top ball
theta2= 30*pi/180 #release angle middle ball
theta1dot = 0
theta2dot = 0
t = 0
dt = 0.001
Next we define all the parts of our system, the pivot point, mass 1 (m1), arm 1 (stick1), and the same for the second set of mass and stick (m2, stick2).
pivot = sphere(pos=vector(0,L1,0), radius=0.05)m1 = sphere(pos=pivot.pos-vector(0,L1,0),radius=0.05,color=color.white)
stick1 = cylinder(pos=pivot.pos,axis=m1.pos-pivot.pos,radius=0.015,color=color.yellow)
m2 = sphere(pos=m1.pos-vector(0,L2,0),radius=0.05,color=color.white)
stick2 = cylinder(pos=m1.pos, axis=m2.pos-m1.pos, radius=0.015,color=color.yellow)
Now we position the masses and sticks so that they reflect the starting angles, as if you ran the code before it would just be two vertical rods. To do this we just use some trigonometry to put the mass the end of where the stick will go, then set the stick to go the distance between each mass. Here we also use VPython’s ‘attach_trail’ method to make the bottom mass generate a trail.
m1.pos = pivot.pos+vector(L1*sin(theta1),-L1*cos(theta1),0)
m2.pos = m1.pos+vector(L2*sin(theta2),-L2*cos(theta2),0)stick1.axis = m1.pos - pivot.pos
stick2.pos = m1.pos
stick2.axis = m2.pos - m1.pos
attach_trail(m2, retain=200, color=color.red)
Next we begin the main loop of our program; in this simulation, as in many physics simulations we’ll use a while loop with an exit condition at t=50. We’ll also use the the rate() method to set the speed that the program will run at.
while t <= 50:
rate(1000)
Now we have to dive into the real math. We’ll use a bit of calculus (not really, but we do deal with derivatives), trigonometry, and algebra to first calculate the angular acceleration, which we will then use to increment the angular velocity, and then apply that angular velocity to the position of the masses and rods to move them according to their calculated accelerations.
To calculate the top mass’s acceleration, let’s use this formula to represent the angular acceleration, or second derivative of θ1:
You can also represent the same equation using the following code:
theta1ddot = (-g*(2*M1 + M2)*sin(theta1) -M2*g*sin(theta1 -2*theta2)-2*sin(theta1 - theta2)*M2*(L2*theta2dot**2 + L1*cos(theta1 - theta2)*theta1dot**2))/(L1*(2*M1 + M2 - M2*cos(2*theta1 -2*theta2)))
For the second mass’s acceleration, we must use a different formula since it represents a different part of the system.
That equation is written in the program as:
theta2ddot = (2*sin(theta1 - theta2)*((M1 + M2)*L1*theta1dot**2 + g*(M1 + M2)*cos(theta1) + L2*M2*cos(theta1 - theta2)*theta2dot**2))/(L2*(2*M1 + M2 - M2*cos(2*theta1 -2*theta2)))
Next let’s increment the variables for angular velocity, ‘theta1dot’ and ‘theta2dot.’ To do this we simply use the equation below, which increments the angular velocity by the angular acceleration times dt. Although the image only states the equation for theta1dot, it is the same for both θ1 and θ2.
Let’s represent this with code for our purposes:
theta1dot = theta1dot + theta1ddot*dt
theta2dot = theta2dot + theta2ddot*dt
For the final part of our calculations for the motion of the double pendulum must increment the angles θ1 and θ2 by the angular velocity. This is represented by the equation below, where we set the angle θ equal to the previous θ plus the angular velocity times dt. Again, the calculation is the exact same for θ1 and θ2.
In code, we write the equation as:
theta1 = theta1 + theta1dot*dt
theta2 = theta2 + theta2dot*dt
Now we must update the position of the masses and the rods in the program, so that we can see the outcome of the program. First we will update the position of the masses. We’ll take an intuitive approach to updating these positions; since one thing that will never change in the system is the length of the rods, the middle mass (m1) will be a distance equal to L1 away from the pivot point. Logically we can then construct the following equation for the position of mass m1:
The same formula can be used to reflect the formula for the position of m2, only with θ2 and L2, and it is based off of the position of m1, not the pivot. See this equation below.
These two formulas are written in our program very simply, since one feature of VPython is that it makes vector calculations very simple. All objects already have their positions stored as vectors, so all you need to do is use the vector() method to define a new vector, and use it in your computation.
m1.pos = pivot.pos + vector(L1*sin(theta1),-L1*cos(theta1),0)
m2.pos = m1.pos + vector(L2*sin(theta2),-L2*cos(theta2),0)
Now we must take one final step to realize the outcome of our calculations by updating the position of the rods. Technically this isn’t totally necessary, as the masses will still show and follow the pattern they should. For visual effect though, let’s write this code. The math isn’t too interesting, so the code is as follows. First we set the axis of the first stick so that it fills the distance between the pivot point and mass m1. Next, we set the second stick’s position equal to mass m1. Finally, we set the axis of the second stick equal to the distance between the first and second mass. The final line is incrementing t by dt to make the program go forwards.
stick1.axis = m1.pos - pivot.pos
stick2.pos = m1.pos
stick2.axis = m2.pos - m1.post += dt
Here is the final code for ease of copying and pasting:
Web VPython 3.2#lengths of the strings
L1 = 1.5 #top string
L2 = 1 #middle string
#masses
M1 = 2 #top moving ball
M2 = 1 #middle moving ball
#g
g = 9.8
t = 0
dt = 0.001
#starting angles and velocities
theta1=80.5*pi/180 #release angle top ball
theta2= 0*pi/180 #release angle middle ball
theta1dot = 0
theta2dot = 0
pivot = sphere(pos=vector(0,L1,0), radius=0.05)
m1 = sphere(pos=pivot.pos-vector(0,L1,0),radius=0.05,color=color.white)
stick1 = cylinder(pos=pivot.pos,axis=m1.pos-pivot.pos,radius=0.015,color=color.yellow)
m2 = sphere(pos=m1.pos-vector(0,L2,0),radius=0.05,color=color.white)
stick2 = cylinder(pos=m1.pos, axis=m2.pos-m1.pos, radius=0.015,color=color.yellow)
m1.pos = pivot.pos+vector(L1*sin(theta1),-L1*cos(theta1),0)
m2.pos = m1.pos+vector(L2*sin(theta2),-L2*cos(theta2),0)
stick1.axis = m1.pos - pivot.pos
stick2.pos = m1.pos
stick2.axis = m2.pos - m1.pos
attach_trail(m2, retain=200, color=color.red)
eChart = gdisplay(x=500, y=0, width=600, height=400, title="", xtitle="", ytitle="", foreground=color.black, background=color.white)
ePlot = gdots(color=color.red)
while t < 50:
rate(1000)
#angular acceleration
theta1ddot = (-g*(2*M1 + M2)*sin(theta1) -M2*g*sin(theta1 -2*theta2)-2*sin(theta1 - theta2)*M2*(L2*theta2dot**2 + L1*cos(theta1 - theta2)*theta1dot**2))/(L1*(2*M1 + M2 - M2*cos(2*theta1 -2*theta2)))
theta2ddot = (2*sin(theta1 - theta2)*((M1 + M2)*L1*theta1dot**2 + g*(M1 + M2)*cos(theta1) + L2*M2*cos(theta1 - theta2)*theta2dot**2))/(L2*(2*M1 + M2 - M2*cos(2*theta1 -2*theta2)))
#angular velocity
theta1dot = theta1dot + theta1ddot*dt
theta2dot = theta2dot + theta2ddot*dt
#angular position
theta1 = theta1 + theta1dot*dt
theta2 = theta2 + theta2dot*dt
m1.pos = pivot.pos + vector(L1*sin(theta1),-L1*cos(theta1),0)
m2.pos = m1.pos + vector(L2*sin(theta2),-L2*cos(theta2),0)
stick1.axis = m1.pos - pivot.pos
stick2.pos = m1.pos
stick2.axis = m2.pos - m1.pos
t = t+dt
Thank you for making it this far, I hope you enjoyed my analysis of chaotic motion through the double and triple pendulums. In the end, I’m neither a professional physicist nor programmer, so if you have suggestions for how I can improve my work please don’t hesitate to let me know.
Additionally, all imaged used in the article were created by me, using Python and Microsoft Word, both of which have incredible capabilities.
Finally, I’d like to end on a chewy thought, which will play into my next project. If you took a multi-armed pendulum which was already experiencing chaotic motion, how would its motion change if you were to seamlessly add another arm of very high mass (maybe 50x the largest already in a system), at a starting θ of 0°?
