Using VPython simulations to model chaotic motion and investigate what defines a chaotic system.

Any physics student has examined the lowly pendulum, probably greater than they desired to. Together with the mass on a spring, pendulums are the quintessential example of easy harmonic motion. The keyword there may be easy. They, as you could know, are incredibly easy to model mathematically. The period of the pendulum might 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 point of fact it’s only a cosine wave. The equation below represents the angular displacement of a single pendulum, where θ is the angular displacement, A is the utmost angular displacement (usually that is just the discharge angle), ω is the angular frequency, and t is time.

With that out of the way in which, let’s get to the more interesting, albeit effectively inconceivable (you’ll see), mathematics and physics.

Odds are, if you happen to’ve worked with a daily pendulum, you’re not less than somewhat conversant in the double pendulum, which is paying homage to the infamous three-body problem. The three-body problem represents the concept that an analytical solution to the motion of three bodies in a system where all bodies’ forces affect one another is inconceivable to resolve, because the forces of every body interact with the others concurrently. In brief, the double pendulum is inconceivable to resolve analytically. By way of the pendulum equations stated above, the one that you simply really can’t solve is the equation for θ, or the angular displacement of each pendulums. You furthermore may cannot solve for the period, but that’s mostly because a double pendulum doesn’t exactly have a daily period. Below is a visible 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 will not be simply random, but slightly chaotic. How can we prove that chaotic and random aren’t the identical? That’s the query we’ll seek to reply afterward.

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, together with many other butterfly wing-esque drawings, are mostly related to chaos theory. Something not commonly regarded as a part of chaos theory is the double, triple, and every other pendulum with a couple of arm. Let’s explore a number of the properties of chaotic motion and take a look at to learn more about it as best we are able to, with help from a pc program that we are going to dive into later.

One property we are able to investigate is the effect that making a small change in starting conditions has on the top result. Will the change be negligible ultimately, or will it yield a solution entirely consequence from the unique?Below is a pc simulation of a double pendulum at t=50, where the starting angles were 80° and 0° respectively for θ1 and θ2. (Angles are the identical as shown within the diagram above)

Now, let’s run the exact same simulation, the image is of t=50, all of the masses and arm lengths are completely equivalent. Nonetheless, the angles are 79° and 0° for θ1 and θ2. There’s a 1° change within the starting value of θ1, let’s see the way it affects the consequence.

Unsurprisingly (or surprisingly, I don’t know what you expected), with a change of just 1° the system has dramatically modified to the purpose where the ultimate position of the system and the road traced by it don’t resemble the primary image in any respect. That’s independent of any outside variables, as in our simulation there isn’t any friction, air resistance, or every other hindrances. This instance shows the concept that when a system exhibits chaotic motion, as in a double pendulum, even essentially the most minute change in conditions could cause an unlimited change within the result.

So, we’ve nailed down one interesting aspect of chaotic motion– a tiny change in starting condition could cause an unlimited change in result. Let’s attempt to discover a number of more. What would occur if we were to run the identical double pendulum experiment again? Experimentally, that is all but inconceivable; it could involve having an infinite degree of precision when measuring things like gravitational field, mass, air resistance, and each other variable within the experiment, after which using your magical infinite precision data to establish the experiment again, identically. We, nevertheless, have a pleasant simulation to make use of which makes that challenge disappear.

Let’s run the double pendulum simulation just like before, but with some variables modified so the image isn’t the identical: t=50, θ1=80°, θ2=15°, and the physical properties of the pendulum are the identical (length, mass). Here is the result:

Spoiler alert, it was the very same the second time… And the third, fourth fifth, and sixth time. Simply to prove a degree, let’s also do this exercise with a triple pendulum. We’ll use the identical starting conditions, although now we now have a 3rd angle, θ3, so we’ll use θ1=80°, θ2=15°, θ3=0°. A triple pendulum is rather more computationally intensive, though, so the image might be taken at t=10.

Even though it is much more complex than the double pendulum simulation, the incontrovertible fact that the triple pendulum exhibits the identical tendency to repeat its pattern under equivalent starting conditions is significant. That concept of the same consequence under equivalent conditions proves that chaotic motion and chaotic systems are deterministic, thus not random. Again, since that is so hard to prove experimentally or use in practice it doesn’t affect our understanding of the actual world very much. We do, nevertheless, now know that chaotic systems are usually not random and might be predicted with the appropriate knowledge. Moreover, since we proved that a small change can have a big effect, we all know that an incredibly high degree of precision is required to effectively predict the consequence of a chaotic system.

We’ve now identified two significant properties of chaotic motion: the primary is that a tiny change can have huge implications within the consequence of the system, and the second is that chaotic motion will not be random but slightly it’s deterministic. Let’s try to seek out another– as we discussed on the very starting of the article, the variable you actually can’t solve for in a double pendulum is the angular displacement. Nonetheless, we also said that solving for period was futile because a double pendulum doesn’t repeat– let’s test that theory.

First, let’s run the double pendulum experiment with the next conditions: θ1=30°, θ2=30°, with the identical masses and lengths as before. This time, let’s use VPython’s graphing capabilities to record a chart of the peak of the underside mass as a function of time:

When you ask me, that appears pretty periodic. So, does this mean that the system described above will not be a chaotic system? And if the motion will not be chaotic, does that mean we could reach an analytical solution to the motion of the double pendulum when it has lower than some threshold energy needed to change into chaotic? I’ll leave that as an exercise to the reader. To try this though, you may need the assistance of this system I wrote and utilized in this text, so let’s walk through it. This program is written in VPython using GlowScript, so make sure that you’re running it on the WebVPython system, and not only any python editor.

First, let’s define the physical constants we’ll use in this system. These include the masses, lengths of arms, and force of gravity. Moreover, 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 period 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 of the parts of our system, the pivot point, mass 1 (m1), arm 1 (stick1), and the identical 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 in order that they reflect the starting angles, as if you happen to ran the code before it could just be two vertical rods. To do that we just use some trigonometry to place the mass the top of where the stick will go, then set the follow go the gap between each mass. Here we also use VPython’s ‘attach_trail’ method to make the underside 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 start the essential loop of our program; on this simulation, as in lots of physics simulations we’ll use some time loop with an exit condition at t=50. We’ll also use the the speed() method to set the speed that this system will run at.

while t <= 50:  
rate(1000)

Now we now have to dive into the actual math. We’ll use a little bit of calculus (probably not, but we do take care of derivatives), trigonometry, and algebra to first calculate the angular acceleration, which we’ll then use to increment the angular velocity, after which apply that angular velocity to the position of the masses and rods to maneuver them in accordance with their calculated accelerations.

To calculate the highest mass’s acceleration, let’s use this formula to represent the angular acceleration, or second derivative of θ1:

It’s also possible to represent the identical equation using the next 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 distinct formula because it represents a distinct a part of the system.

That equation is written in this system 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 that 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 similar for each θ1 and θ2.

Let’s represent this with code for our purposes:

theta1dot = theta1dot + theta1ddot*dt
theta2dot = theta2dot + theta2ddot*dt

For the ultimate a part of our calculations for the motion of the double pendulum must increment the angles θ1 and θ2 by the angular velocity. That 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 very 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 this system, in order that we are able to see the consequence of this system. First we’ll update the position of the masses. We’ll take an intuitive approach to updating these positions; since one thing that may never change within the system is the length of the rods, the center mass (m1) might be a distance equal to L1 away from the pivot point. Logically we are able to then construct the next equation for the position of mass m1:

The identical formula might be used to reflect the formula for the position of m2, only with θ2 and L2, and it is predicated 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 quite simple. All objects have already got their positions stored as vectors, so all you could do is use the vector() method to define a brand 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 comprehend the consequence of our calculations by updating the position of the rods. Technically this isn’t totally vital, because the masses will still show and follow the pattern they need to. For visual effect though, let’s write this code. The maths isn’t too interesting, so the code is as follows. First we set the axis of the primary stick in order that it fills the gap 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 gap between the primary and second mass. The ultimate line is incrementing t by dt to make this system go forwards.

stick1.axis = m1.pos - pivot.pos
stick2.pos = m1.pos
stick2.axis = m2.pos - m1.pos
t += dt

Here is the ultimate 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

Thanks for making it this far, I hope you enjoyed my evaluation of chaotic motion through the double and triple pendulums. In the long run, I’m neither an expert physicist nor programmer, so if you’ve gotten suggestions for the way I can improve my work please don’t hesitate to let me know.

Moreover, all imaged utilized in the article were created by me, using Python and Microsoft Word, each of which have incredible capabilities.

Finally, I’d wish to end on a chewy thought, which is able to play into my next project. When you took a multi-armed pendulum which was already experiencing chaotic motion, how would its motion change if you happen to were to seamlessly add one other arm of very high mass (perhaps 50x the biggest already in a system), at a starting θ of 0°?