Mathematical Drawings with Pendulums

Look at this harmonograph; every time I do, it makes me laugh

A dynamical system is a mathematical model that shows how a set of variables evolves over time from their initial conditions according to an internal ruleset. The variables’ values at a single snapshot in time uniquely define the system’s state; the internal ruleset, written as a system of differential equations, describe how the variables will change in relation to their current values.

Dynamical systems play roles in lots of applied math fields like biology, epidemiology, physics, and engineering. Pendulums are great examples of dynamical systems;¹ let’s use some to understand how they work, and draw some pretty pictures in the process.

using Dates
using OrdinaryDiffEq
using Plots
using Printf

Torque Derivation for a Single Simple Pendulum

Let’s imagine using a string with negligible weight and length \(\ell\) to hang a small, dense plumb bob with mass \(m\) from the ceiling. If we pull the plumb bob back and let it go, we’ve formed a simple pendulum like the figure below. What information do we know about it that we can use to model its movement as a system of differential equations?

Let’s think of the pivot point where the pendulum swings from as the origin of a coordinate system. With this, we know that the plumb bob will experience gravitational force \(F_{g} = -mg\). Let’s also call the angle the string makes with the y-axis \(\theta\); keep in mind that this angle \(\theta\) changes over time \(t\).²

The simple pendulum: The point mass \(m\) connected to a massless rod with length \(\ell\) swings back and forth. The mass \(m\) is pulled down by the gravitational force \(F = mg\). The angle the string makes with the vertical axis (moving counter-clockwise) is \(\theta\).

There are several ways to derive the formulas for a simple pendulum, but the one that makes the most sense to me uses two definitions of torque to do it. Torque is the twisting effectiveness of a force around a center point. The torque pseudo-vector \(\vec{\tau}\) at a point relative to an origin is equal to the cross product of the force vector \(\vec{F}\) at that point and the displacement vector \(\vec{r}\) from the origin. Using the definition of the cross product, we can also define the magnitude (or size) of the torque pseudo-vector in terms of the distance from the origin \(r\), force magnitude \(F\), and the angle \(\theta\) between them.³

\[ \begin{aligned} \vec{\tau} &= \vec{r} \times \vec{F} \\ \left\|\vec{\tau}\right\| &= r F \sin(\theta) \end{aligned} \]

NoteSide Note

Torque is a pseudo-vector because it behaves differently than normal vectors. A torque pseudo-vector must be parallel to its axis of rotation, and the direction in which it points depends on the frame of reference (i.e., the way in which you view it). When you view the force vector pushing in a counter-clockwise direction, the torque points towards you; when the force vector pushes in a clockwise direction, the torque points away from you. The right-hand rule is an easy way to remember which direction the torque vector points.

A good example of torque is using a hex key to screw in a bolt. Turning the hex key to screw the bolt in with the long side out is easier than turning it with the short part out. Increasing the distance between the applied force and the origin increases the torque.

Since torque is the rotational equivalent of force, it can also be derived with the moment of inertia \(I\) and angular acceleration \(\alpha\) using Newton’s Second Law for Rotation. The plumb bob’s angular acceleration \(\alpha\) is the derivative (instantaneous rate of change) of its angular velocity \(\omega\), which in turn is the rate of change of its angular displacement \(\theta\).⁴

\[ \left\|\vec{\tau}\right\| = I \alpha \]

An object’s moment of inertia is its resistance to angular acceleration, like how mass is to normal acceleration. The formula for an object’s moment of inertia is dependent on its shape. If we assume that our pendulum’s plumb bob is a point mass, that will make our moment of inertia equal to \(mr^2\).⁵

\[ \text{For point mass:} \quad I = m r^2 \]

To recap:

Symbol	Definition	Units
\(\ell\)	Length of the string	meters, \(\text{m}\)
\(m\)	Mass	kilograms, \(\text{kg}\)
\(\theta\)	Angular displacement	radians, \(\text{rad}\)
\(F\)	Force magnitude	newtons, \(\text{N}\)
\(g\)	Acceleration due to gravity	meters per second squared, \(\text{m}/\text{s}^2\)
\(r\)	Displacement vector from rotation axis to applied force location	meters, \(\text{m}\)
\(\vec{\tau}\)	Torque pseudo-vector	netwton meters, \(\text{N}\cdot\text{m}\)
\(\alpha\)	Angular acceleration magnitude	radians per second squared, \(\text{rad}/\text{s}^2\)
\(\omega\)	Angular velocity magnitude	radians per second, \(\text{rad}/\text{s}\)
\(I\)	Moment of inertia	kilogram meters squared, \(\text{kg}\cdot\text{m}^2\)

We now have all the parts needed to derive the motion of the simple pendulum.

Statement	Reason
\(\vec{\tau} = \vec{r} \times \vec{F}\)	Definition of Torque Vector
\(\left\|\vec{\tau}\right\| = rF \sin(\theta)\)	Torque Vector Magnitude
\(\left\|\vec{\tau}\right\| = \ell (-mg) \sin(\theta)\)	Substitute \(F = -mg\)
\(I \alpha = -\ell mg \sin(\theta)\)	Newton’s Second Law of Rotation
\(I \frac{d^2 \theta}{dt^2} = -\ell mg \sin(\theta)\)	Definition of Angular Acceleration
\(m\ell^2 \frac{d^2 \theta}{dt^2} = -\ell mg \sin(\theta)\)	Moment of Inertia for a point mass
\(\frac{d^2 \theta}{dt^2} = -\frac{g}{\ell}\sin(\theta)\)	Simplify
\(\frac{d^2 \theta}{dt^2} + \frac{g}{\ell}\sin(\theta) = 0\)	Rearrange

Our final equation shows us that the pendulum’s mass has no effect on the path that it takes over time; the mass \(m\) dropped out of the equation entirely!

Visualizing a Damped Pendulum in Julia

Our assumptions in the previous section didn’t include how friction (like air resistance or energy loss from the pivot point) affects our pendulum. Let’s assume that our pendulum is damped by some drag amount \(k\) that is proportional to the angular velocity \(\omega = \frac{d\theta}{dt}\).⁶

\[ \underbrace{\frac{d^2 \theta}{dt^2}}_{\frac{d\omega}{dt}} + k\underbrace{\frac{d\theta}{dt}}_\omega + \frac{g}{\ell}\sin(\theta) = 0 \]

This equation is a second-order ordinary differential equation (ODE) because its highest derivative is the second derivative and the equation is only dependent on the angle \(\theta\). Let’s define our constants and initial conditions here. We’ll run a pendulum simulation for six seconds.

# Constants
g = 9.81
ℓ = 1.0
k = 0.5
plumbbob_rad = 0.125

# Initial Conditions
θ₀ = π / 2
ω₀ = 0
u₀ = [θ₀, ω₀]
tspan = (0.0, 6.0)
dt = 1 / 20

We need to rearrange our second-order ODE equation it into a system of first-order ODEs to input into our ODE solving program. Let’s substitute \(\omega\) back into our equation and rewrite using \(\omega\), \(\theta\), their derivatives, and constants.

\[ \begin{aligned} \frac{d\omega}{dt} + k\omega + \frac{g}{\ell}\sin(\theta) &= 0 \\ \frac{d\omega}{dt} &= -\frac{g}{\ell}\sin(\theta) - k\omega \end{aligned} \]

Now, with both derivatives defined, we can use linear algebra to package them into a system of first-order ODEs. The vector \(\vec{u}\) contains our dependent variables. We can see the structure of the ODE system by taking the derivative of \(\vec{u}\).

\[ \begin{aligned} \vec{u} &= \begin{bmatrix} \theta \\ \omega \end{bmatrix} \\ \frac{d\vec{u}}{dt} &= \begin{bmatrix} \frac{d \theta}{dt} \\ \frac{d \omega}{dt} \end{bmatrix} &= \begin{bmatrix} \omega \\ -\frac{g}{\ell} \sin(\theta) - k\omega \end{bmatrix} \end{aligned} \]

We can define the relationship between \(\vec{u}\) and \(\frac{d\vec{u}}{dt}\) in the function below.⁷

"Describe the movement of a simple pendulum as an ODE system."
function simple_pendulum!(du, u, _, _)
    θ = u[1]
    ω = u[2]
    du[1] = ω
    du[2] = -(g / ℓ) * sin(θ) - k * ω
    return nothing
end

simple_pendulum!

There are many numerical methods to solve ODE systems. The easiest is Euler’s method: it estimates the equation’s solutions by computing the dependent variables’ changes in incremental discrete-time steps.⁸ I used it in Genuary 2024 prompt 8 to visualize the Lorenz system. The Runge-Kutta (RK) family of methods is more popular; it improves on how Euler’s method approximates these discrete-time step changes. I’ll use the Tsit5 solver to solve this equation, which further improves upon the RK4 method.

To visualize the pendulum, we’ll need its \((x, y)\) coordinates. We can calculate these using \(\theta\) and \(\ell\).

sp_prob = ODEProblem(simple_pendulum!, u₀, tspan)
sp_sol = solve(sp_prob, Tsit5(); saveat=dt)

u = reduce(hcat, sp_sol.u)'
t = sp_sol.t
θ = u[:, 1]
ω = u[:, 2]

x = ℓ * cos.(θ .- π / 2)
y = ℓ * sin.(θ .- π / 2)

Now, let’s see the pendulum in action!

# Draw the pendulum pivot
pivot_width = 0.25
pivot_height = 0.125
pivot = Shape([
    (-pivot_width / 2, 0),
    (-pivot_width / 2, pivot_height),
    (pivot_width / 2, pivot_height),
    (pivot_width / 2, 0),
])

# Initialize angular displacement plot
angle = plot(
    1;
    legend=false,
    xlims=tspan,
    ylabel="θ(t) [rad]",
    ylims=(minimum(θ), maximum(θ)),
    guidefont=(10, :black),
    framestyle=:origin,
)

# Initialize angular velocity plot
ang_vel = plot(
    1;
    legend=false,
    xlabel="Time, t [s]",
    xlims=tspan,
    ylabel="ω(t) [rad/s]",
    ylims=(minimum(ω), maximum(ω)),
    guidefont=(10, :black),
    framestyle=:origin,
)

anim = Animation()
for i in eachindex(t)
    # Side note: moving this layout macro outside the for-loop will break the animation
    pendulum_layout = @layout [
        a [grid(2, 1)]
    ]
    # Update angular displacement & velocity plots
    push!(angle, t[i], θ[i])
    push!(ang_vel, t[i], ω[i])
    # Pendulum Motion Plot
    pendulum = plot(
        [0; x[i]],
        [0; y[i]];
        linewidth=:7,
        aspect_ratio=:equal,
        xlabel="Horiz. Displacement, x [m]",
        xlims=(-ℓ - plumbbob_rad, ℓ + plumbbob_rad),
        ylabel="Vert. Displacement, y [m]",
        ylims=(-ℓ - plumbbob_rad, ℓ + plumbbob_rad),
        guidefont=(10, :black),
        framestyle=:box,
        title="Pendulum Simulation",
    )
    scatter!([x[i]], [y[i]]; markersize=plumbbob_rad * 96, legend=false)
    plot!(pivot; fillcolor="gray")
    # Bring it all together!
    plot(
        pendulum,
        angle,
        ang_vel;
        layout=pendulum_layout,
        size=(800, 400),
        guidefont=(10, :black),
        left_margin=4Plots.mm,
        bottom_margin=4Plots.mm,
    )
    frame(anim)
end
gif(anim, "pendulum_plot.gif"; fps=Int(1 / dt))

Completed pendulum graph. Notice what the pendulum is doing when the angular displacement \(\theta(t)\) and angular velocity \(\omega(t)\) are equal to zero.

2+ Pendulums Make A Harmonograph

If we plot the motions from two independently operating pendulums on different axes over time, we’ll get complex harmonic motion. These plots are harmonograms, and the devices that make them are harmonographs. If the ratio between the two pendulums’ oscillation frequencies are (close enough to) neat whole-number ratios and the drag constants \(k_a, k_b\) are low enough, the resulting harmonograms resemble Lissajous curves.⁹ With higher drag constants, the pendulums’ swings shrink over time and settle at the center point.

Harmonograph designs usually have three pendulums. Let’s instead simulate a simpler two-pendulum harmonograph by expanding upon our previously-derived pendulum formula. Now, instead of one pendulum to keep track of, we have two; I’ve called them “Pendulum A” and “Pendulum B”. We can modify our definition of \(\vec{u}\) to run both pendulum systems in parallel.

\[ \begin{aligned} \vec{u} &= \begin{bmatrix} \theta_a \\ \omega_a \\ \theta_b \\ \omega_b \end{bmatrix} \\ \frac{d\vec{u}}{dt} &= \begin{bmatrix} \frac{d\theta_a}{dt} \\ \frac{d\omega_a}{dt} \\ \frac{d\theta_b}{dt} \\ \frac{d\omega_b}{dt} \end{bmatrix} &= \begin{bmatrix} \omega_a \\ -\frac{g}{\ell_a} \sin(\theta_a) - k\omega_a \\ \omega_b \\ -\frac{g}{\ell_b} \sin(\theta_b) - k\omega_b \end{bmatrix} \end{aligned} \]

# Constants
g = 9.81
ℓa, ℓb = 1.0, 0.8
ka, kb = 0.09, 0.1
tspan = (0.0, 25.0)

# Initial Conditions
u₀ = [π / 3, 0, -π / 3, 0]
dt = 1 / 20

"Describe the movement of a two-pendulum harmonograph as an ODE system."
function harmonograph!(du, u, _, _)
    θa = u[1]
    ωa = u[2]
    θb = u[3]
    ωb = u[4]
    du[1] = ωa
    du[2] = -(g / ℓa) * sin(θa) - ka * ωa
    du[3] = ωb
    du[4] = -(g / ℓb) * sin(θb) - kb * ωb
    return nothing
end

# Solve the ODE system
h_prob = ODEProblem(harmonograph!, u₀, tspan)
h_sol = solve(h_prob, Tsit5(); saveat=dt)

# Get the ODE systems' time series
u = reduce(hcat, h_sol.u)'
t = h_sol.t
θa = u[:, 1]
ωa = u[:, 2]
θb = u[:, 3]
ωb = u[:, 4]

xa = ℓa * cos.(θa .- π / 2)
xb = ℓb * cos.(θb .- π / 2)

# Get range of x-values for plotting bounds
x_range = maximum(abs.(vcat(xa, xb)))

# Initialize harmonograph plot and animate it
h_plot = plot(
    1;
    linewidth=2,
    xlim=(-x_range, x_range),
    ylim=(-x_range, x_range),
    legend=false,
    aspect_ratio=:equal,
    size=(600, 600),
    framestyle=:box,
    xlabel="Pendulum A Horiz. Displacement, xa [m]",
    ylabel="Pendulum B Horiz. Displacement, xb [m]",
    title="Harmonograph",
)
anim = @animate for i in eachindex(t)
    push!(h_plot, 1, (xa[i], xb[i]))
end
gif(anim, "harmonograph_plot.gif"; fps=Int(1 / dt))

This parametric graph of a harmonogram is made by matching up two pendulums’ \(x\)-coordinates over time.

The Spherical Pendulum

Let’s once again imagine using a string with negligible weight and length \(\ell\) to hang a small, dense plumb bob with mass \(m\) from the ceiling. Instead of letting it go, let’s give it a small push in some direction different from the one it would swing in. Now, instead of a simple pendulum, we have a spherical pendulum.¹⁰

The spherical pendulum: The point mass \(m\) connected to a massless rod with length \(\ell\) swings in multiple directions. The azimuthal angle \(\phi\) denotes its rotation around the vertical axis. The polar angle the string makes moving outward from the vertical axis is \(\theta\). Image Source: Wikimedia Commons; original image by Maschen.

Deriving the equations for the spherical pendulum is a lot harder than the two-dimensional simple pendulum case; in fact, I couldn’t do it myself. The way that made the most sense to me requires Hamiltonian mechanics, which uses how the system’s total energy is constant to keep track of the system’s generalized momentum components. If you want to learn more about this, I highly recommend checking out Prof. Martin Houde’s lecture notes in his classical mechanics course. Out of all the sources that I could find on the spherical pendulum, his was the most understandable.¹¹

The generalized momentum components for the azimuth angle \(\phi\) and polar angle \(\theta\) are \(p_\phi\) and \(p_\theta\), respectively.¹² These momenta are required for the ODE system to function, and are defined below.

\[ \begin{aligned} p_\theta &= m \ell^2 \theta^{\prime} \\ p_\phi &= m \ell^2 \sin^2(\theta) \phi^{\prime} \end{aligned} \]

These equations are in terms of the pendulum’s mass \(m\), the string length \(\ell\), the spherical coordinate angles \(\theta\) and \(\phi\), and their derivatives \(\theta^{\prime}\) and \(\phi^{\prime}\).¹³ I’ll use them to determine part of the system’s initial conditions.

We can define the system in its entirety below.

\[ \begin{aligned} \vec{u} &= \begin{bmatrix} \theta \\ \phi \\ p_\theta \\ p_\phi \end{bmatrix} \\ \frac{d\vec{u}}{dt} &= \begin{bmatrix} \theta^{\prime} \\ \phi^{\prime} \\ p^{\prime}_\theta \\ p^{\prime}_\phi \end{bmatrix} &= \begin{bmatrix} \frac{p_\theta}{m \ell^2} \\ \frac{p_\phi}{m \ell^2 \sin^2(\theta)} \\ \frac{p_\theta^2 \cos(\theta)}{m \ell^2 \sin^3(\theta)} - m g \ell \sin(\theta) \end{bmatrix} \end{aligned} \]

In addition to solving the system, I’ve also added the shape that the \(x\)- and \(y\)-coordinates trace over time.

# Constants
g = 9.81
ℓ = 1.0
m = 1.0

# Initial Conditions
θ₀ = π / 4
dθ₀ = 1
ϕ₀ = π / 4
dϕ₀ = 1
pθ₀ = m * ℓ^2 * dθ₀
pϕ₀ = m * ℓ^2 * sin(θ₀)^2 * dϕ₀
u₀ = [θ₀, ϕ₀, pθ₀, pϕ₀]
tspan = (0.0, 10.0)
dt = 1 / 20

"Describe the movement of a spherical pendulum as an ODE system."
function spherical_pendulum!(du, u, _, _)
    θ = u[1]
    ϕ = u[2]
    pθ = u[3]
    pϕ = u[4]
    du[1] = pθ / (m * ℓ^2)
    du[2] = pϕ / (m * ℓ^2 * sin(θ)^2)
    du[3] = (pϕ^2 * cos(θ)) / (m * ℓ^2 * sin(θ)^3) - m * g * ℓ * sin(θ)
    du[4] = 0
    return nothing
end

# Solve the ODE system
h_prob = ODEProblem(spherical_pendulum!, u₀, tspan)
h_sol = solve(h_prob, Tsit5(); saveat=dt)

# Obtain the ODE systems' time series
u = reduce(hcat, h_sol.u)'
t = h_sol.t
θ = u[:, 1]
ϕ = u[:, 2]
x = ℓ .* sin.(θ) .* cos.(ϕ)
y = ℓ .* sin.(θ) .* sin.(ϕ)
z = -ℓ .* cos.(θ)

# Get the current color palette for plotting
palette = theme_palette(:auto)

anim = @animate for i in eachindex(t)
    plot(
        x[1:i],
        y[1:i],
        fill(-1, i);
        linecolor=palette[3],
        linestyle=:dot,
        xlims=(-1, 1),
        ylims=(-1, 1),
        zlims=(-1, 1),
        xlabel="x [m]",
        ylabel="y [m]",
        zlabel="z [m]",
        legend=false,
    )
    plot!([0; x[i]], [0; y[i]], [0; z[i]]; linecolor=palette[1])
    scatter!([x[i]], [y[i]], [z[i]]; color=palette[2])
end
gif(anim, "spherical_pendulum_plot.gif"; fps=Int(1 / dt))

Animated spherical Pendulum Plot with additional xy-plane trace.

Conclusion

In this article, we’ve demonstrated how both simple and spherical pendulums are dynamical systems. They can be represented as systems of differential equations that control how properties change over time. While the simple pendulum’s differential equations can be derived with high-school level physics, the spherical pendulum’s motions are much harder to model.

Package Versions & Last Run Date

Julia 1.11.7
OrdinaryDiffEq 6.102.1
Plots 1.41.1
Last Run 2025-10-04T12:51:03.435

References

Houde, Martin. “Hamiltonian Dynamics.” In Physics 350/Applied Math 353: Classical Mechanics I, 91–118, 2023. https://physics.uwo.ca/~mhoude2/courses/physics_350.html.

Lebl, Jiří. Notes on Diffy Qs: Differential Equations for Engineers. 6.6 ed. Lulu, 2023. https://www.jirka.org/diffyqs/.

Moebs, William, Samuel J. Ling, and Jeff Sanny. University Physics. Vol. 1. University Physics. Houston, Texas: OpenStax College, Rice University, 2016. https://openstax.org/details/books/university-physics-volume-1.

Sayama, Hiroki. Introduction to the Modeling and Analysis of Complex Systems. Deluxe color ed. Geneseo, NY: Open SUNY Textbooks, Milne Library, 2015. https://milneopentextbooks.org/introduction-to-the-modeling-and-analysis-of-complex-systems/.

Sims, Karl. “How to Make a Three-Pendulum Rotary Harmonograph,” 2009. http://www.karlsims.com/harmonograph/index.html.

Footnotes

1. Sayama, Introduction to the Modeling and Analysis of Complex Systems.

2. Moebs, Ling, and Sanny, University Physics, vol. 1, chap. 15.4.

3. Moebs, Ling, and Sanny, vol. 1, chap. 15.4.

4. Moebs, Ling, and Sanny, vol. 1, chap. 10.6.

5. Moebs, Ling, and Sanny, University Physics.

6. Lebl, Notes on Diffy Qs.

7. This function has four inputs to match the DifferentialEquations.jl ODE problem definition. The two other inputs that I’m not using are for parameters p (which I’m defining as constants) and time t (which isn’t needed because this ODE isn’t time-dependent).

8. Lebl, Notes on Diffy Qs, chap. 1.7.

9. Sims, “How to Make a Three-Pendulum Rotary Harmonograph”.

10. Houde, “Hamiltonian Dynamics”.

11. Houde.

12. Houde.

13. Houde.