Watch it fly by as the pendulum swings
First things first: Yes, that’s a Linkin Park lyric. I’m not proud of it but there aren’t enough pop-culture pendulum references.
I think most people understand what a pendulum is. It’s when you have a heavy thing that swings on a cable or rod, like an old grandfather clock. In this post we’ll assume it’s an “idealized” pendulum. That means the rod is massless, the heavy thing is a point mass, and there’s no friction or air resistance. It almost seems like cheating but it makes the math a lot easier.
Let’s skip to the relevant equations.
1. Theory.
The angle theta (θ) that a pendulum makes with the vertical axis, as a function of time, can be expressed like this:
Where…
- θ0 is the initial angle (radians)
- g is the “little g” gravitational constant (9.81 m/s2)
- L is the rod length (meters)
- t is time (seconds)
The catch is that this formulation uses the small angle approximation, which is indeed an approximation. And it only holds true up to 0.5 radians or so. If we pulled the pendulum back further and released it, our model would begin to diverge from reality.
Here’s a GIF to illustrate exactly what theta is measuring:
Theta is the angle formed by the pendulum rod and the vertical axis. Usually we measure angles with respect to a horizontal axis, so keep that in mind for later.
2. Code.
The plan is to draw two animated plots side-by-side. On the left will be a pendulum swinging back and forth like the GIF above. On the right will be theta as a function of time. As we saw from the formula, it’s a cosine function, so we’ll be plotting a waveform.
Most of the code will be nested inside a large for loop. We’ll create one frame of the animation during each iteration of the loop. Once all the images are generated, we can open them in image editing software—I’ll use GIMP—and save them as an animated GIF.
Start by declaring the constants. We’ll use a rod length of 2 meters to keep things interesting. Make the initial angle (THETA_INITIAL) 0.5 radians to satisfy the small angle approximation. The gravitational constant is always 9.81 m/s2. T_STEP is how much time (t) will elapse between frames. The smaller the number, the smoother the animation, but file size increases as well. A value of 0.02 implies 50 frames per second, which looks good to me.
284 frames might seem like a random number but it’s not. A pendulum’s period, i.e. one full cycle from left to right and back again, is defined like this:
If you do the math, one cycle takes about 2.84 seconds. I think it will be good to animate two complete cycles, which at 0.02 seconds per frame works out to 284 frames.
t is our time variable. Naturally it will begin at 0.0 seconds.
from math import sqrt, sin, cos
from numpy import arange
import matplotlib.pyplot as plt
LENGTH = 2
THETA_INITIAL = 0.5
G_CONSTANT = 9.81
T_STEP = 0.02
NUMBER_OF_FRAMES = 284
plt.style.use("wollen_dark.mplstyle")
t = 0.0
t_list = []
theta_list = []
for frame_number in range(NUMBER_OF_FRAMES):
print(f"{frame_number + 1}/{NUMBER_OF_FRAMES}")
[...]
Now we can begin plotting inside our loop. Create a figure with two Axes in a horizontal layout—1 row, 2 columns.
Implement the equation above to calculate theta as a function of time. For the plot on the right, ax1, I’d like to trace out theta’s curve piece by piece. To do that, append the current values of theta and t to their corresponding lists. This will save all of theta’s history up to the current point in time.
for frame_number in range(NUMBER_OF_FRAMES):
[...]
fig, (ax0, ax1) = plt.subplots(1, 2)
theta = THETA_INITIAL * cos(sqrt(G_CONSTANT / LENGTH) * t)
t_list.append(t)
theta_list.append(theta)
[...]
With theta’s value in hand, let’s jump back to ax0 and plot the pendulum. It will hang from the origin (0, 0). We’ll do this in two parts: first the rod, then the mass.
Since we know the rod’s angle and length, we can use basic trigonometry to calculate its position. Notice that the x-component of the rod uses sin and the y-component uses cos. It’s the reverse of how this usually works. That’s because theta forms an angle with the vertical axis rather than the horizontal one. It’s a little unusual but it’s perfectly fine.
I’m hard-coding axis arguments for brevity’s sake. We could parameterize them based on LENGTH and THETA_INITIAL but the code would quickly become difficult to read.
for frame_number in range(NUMBER_OF_FRAMES):
[...]
ax0.plot([0, 0 - LENGTH * sin(theta)],
[0, 0 - LENGTH * cos(theta)],
color="#988ED5")
ax0.scatter([0 - LENGTH * sin(theta)],
[0 - LENGTH * cos(theta)],
color="#988ED5",
s=400)
ax0.set(xticks=arange(-1.5, 2.0, 0.5),
xlim=(-1.6, 1.6),
yticks=arange(-2.5, 1.0, 0.5),
ylim=(-2.6, 0.6))
[...]
Now jump to the right-hand side and plot theta. Remember that we’re continuously updating a ledger of values so we just need to pass those lists to plot.
I think it looks nice to call scatter and place a marker at the leading edge of the curve. Simply pass the most recent values of t and theta, repackaged as lists.
The y-axis is in units of radians so it would be nice to set ticks in multiples of π. But we’re dealing with very small angles so I think it’s better to stick with basic floats.
for frame_number in range(NUMBER_OF_FRAMES):
[...]
ax1.plot(t_list,
theta_list,
color="#B0E441")
ax1.scatter([t_list[-1]],
[theta_list[-1]],
color="#B0E441",
s=50)
ax1.set(xticks=arange(0.0, 7.0, 1.0),
xlim=(-0.2, 6.2),
yticks=arange(-0.5, 0.75, 0.25),
ylim=(-0.53, 0.53))
[...]
Set axis labels. We’re plotting a lot of different units so it’s good to be very explicit.
Then save the figure. I’m dumping image files into a “frames” folder.
Before exiting the loop, remember to increment the time variable, t.
for frame_number in range(NUMBER_OF_FRAMES):
[...]
ax0.set(xlabel="Distance (m)",
ylabel="Distance (m)",
title="Pendulum Position")
ax1.set(xlabel="Time (s)",
ylabel="Radians",
title="Angle θ")
plt.savefig(f"frames/{frame_number}.png")
t += T_STEP
After running the script, we’ll have 284 PNG images in a “frames” folder. I’ll open them in GIMP and export an animated GIF.
Another option is to use FFMPEG and automate the GIF creation process. It’s a command line tool but there are Python packages to help you call it directly from a script.
3. The output.
We initially released the pendulum from an angle of 0.5 radians. There’s no friction in our idealized system so it swings exactly that far in the other direction. Repeat ad infinitum.
Download the Matplotlib style.
Full code:
from math import sqrt, sin, cos
from numpy import arange
import matplotlib.pyplot as plt
LENGTH = 2
THETA_INITIAL = 0.5
G_CONSTANT = 9.81
T_STEP = 0.02
NUMBER_OF_FRAMES = 284
plt.style.use("wollen_dark.mplstyle")
t = 0.0
t_list = []
theta_list = []
for frame_number in range(NUMBER_OF_FRAMES):
print(f"{frame_number + 1}/{NUMBER_OF_FRAMES}")
fig, (ax0, ax1) = plt.subplots(1, 2)
theta = THETA_INITIAL * cos(sqrt(G_CONSTANT / LENGTH) * t)
t_list.append(t)
theta_list.append(theta)
ax0.plot([0, 0 - LENGTH * sin(theta)],
[0, 0 - LENGTH * cos(theta)],
color="#988ED5")
ax0.scatter([0 - LENGTH * sin(theta)],
[0 - LENGTH * cos(theta)],
color="#988ED5",
s=400)
ax0.set(xticks=arange(-1.5, 2.0, 0.5),
xlim=(-1.6, 1.6),
yticks=arange(-2.5, 1.0, 0.5),
ylim=(-2.6, 0.6))
ax1.plot(t_list,
theta_list,
color="#B0E441")
ax1.scatter([t_list[-1]],
[theta_list[-1]],
color="#B0E441",
s=50)
ax1.set(xticks=arange(0.0, 7.0, 1.0),
xlim=(-0.2, 6.2),
yticks=arange(-0.5, 0.75, 0.25),
ylim=(-0.53, 0.53))
ax0.set(xlabel="Distance (m)",
ylabel="Distance (m)",
title="Pendulum Position")
ax1.set(xlabel="Time (s)",
ylabel="Radians",
title="Angle θ")
plt.savefig(f"frames/{frame_number}.png")
t += T_STEP