Numerically modelling the simple pendulum I

Activity 1 | Week 3 | Lecture

1 What to expect

1.1 Story so far

Over the past two weeks, we have spent a great deal of time measuring and modelling the period of a pendulum. Our measurements provided evidence that the periods at \(10^\circ\) and \(20^\circ\) are different.

To strengthen this experimental result, we developed the mathematical model

\[ \frac{d^2\theta}{dt^2} = -\frac{g}{L}\sin\theta. \tag{1}\]

We then explored the small-angle approximation, which predicts

\[ T=2\pi\sqrt{\frac{L}{g}}. \tag{2}\]

This is a wonderfully simple result, but the amplitude does not appear in it. The approximation therefore predicts that \(T_{10^\circ}\) and \(T_{20^\circ}\) are identical. Yet our measurements provide some evidence to the contrary.

So, we must return to Equation 1. Solving it analytically requires mathematics that is beyond what we need here. Instead, I will show you how to construct a very accessible numerical solution that can support and complement our measurements.

1.2 Euler to the rescue

Euler’s method, which you first encountered in 73 and will use again for Daisyworld in 75, begins with a simple idea. Over a small interval \(\Delta x\), we can write

\[ \frac{dy}{dx} \approx \frac{\Delta y}{\Delta x}. \]

Rearranging,

\[ \Delta y \approx \frac{dy}{dx}\Delta x. \tag{3}\]

In words: if we know the present value of \(y\) and its rate of change, we can estimate what \(y\) will be after a small step.

For the pendulum, let

\[ \omega=\frac{d\theta}{dt} \]

be the instantaneous angular velocity. We can then rewrite Equation 1 as

\[ \begin{aligned} \frac{d\omega}{dt} &= -\frac{g}{L}\sin\theta,\\ \frac{d\theta}{dt} &= \omega. \end{aligned} \]

Over a small time step \(\Delta t\),

\[ \begin{aligned} \Delta\omega &\approx -\frac{g}{L}\sin\theta\,\Delta t,\\ \Delta\theta &\approx \omega\,\Delta t. \end{aligned} \tag{4}\]

These two expressions give us a numerical strategy: update the angular velocity, update the angle, and repeat.

2 Warming up

Before tackling the pendulum, let’s warm up by first using Euler’s method on a simpler problem.

Imagine a box of mass \(m\) pushed by a constant force \(F\). The box has a constant acceleration, so its speed continually increases. As it speeds up, it also travels further during each time step.

Newton’s second law gives

\[ F = ma = m\frac{dv}{dt} = m\frac{d^2x}{dt^2}. \]

2.1 Solving analytically

I chose this toy problem because we can easily solve it analytically.

Since

\[ a=\frac{F}{m}, \tag{5}\]

and the box begins at position \(x_0\) with velocity \(u\), the solution is

\[ \begin{aligned} v(t)&=u+at,\\[1em] x(t)&=x_0+ut+\frac{1}{2}at^2. \end{aligned} \tag{6}\]

Velocity

Beginning with the acceleration,

\[ \frac{dv}{dt}=a. \]

Integrating from the initial velocity \(u\) at \(t=0\) to the velocity \(v(t)\) at time \(t\),

\[ \begin{aligned} \int_u^{v(t)}dv' &= \int_0^t a\,d{t},\\ v(t)-u &= at. \end{aligned} \]

Therefore,

\[ v(t)=u+at. \]

Displacement

Since \(v=dx/dt\),

\[ \frac{dx}{dt}=u+at. \]

Integrating again,

\[ \begin{aligned} \int_{x_0}^{x(t)}dx &= \int_0^t (u+at)\,dt,\\ x(t)-x_0 &= ut+\frac{1}{2}at^2. \end{aligned} \]

Therefore,

\[ x(t)=x_0+ut+\frac{1}{2}at^2. \]

2.2 Solving using Euler

To use Euler’s method:

Step 1: Re-write the motion as two first-order equations:

\[ \begin{align} a&=\frac{dv}{dt}=\frac{F}{m}\\[2em] v&=\frac{dx}{dt}. \end{align} \tag{7}\]

Step 2: Over a small time step \(\Delta t\) approximate,

\[ \begin{align} \Delta v & \approx\frac{F}{m}\Delta t,\\[2em] \Delta x & \approx v\Delta t. \end{align} \tag{8}\]

Implementing in Python

Let’s implement Euler in Python for our toy problem. Please note how we use the state now (speed and acceleration) to estimate the next state.

Here is the typical code.

force = 10.0       # N
mass = 2.0         # kg
dt = 0.01          # s

acceleration = force / mass

time = 0.0
initial_position = 0.0
initial_velocity = 0.0

position = initial_position
velocity = initial_velocity

# =========== Simulate using  Euler =========== #
history = [[time, position, velocity]]

while position < 1000:  # Stop after travelling 1 km
    position += velocity * dt
    velocity += acceleration * dt
    time += dt

    history.append([time, position, velocity])
# ============================================ #

# Extract the simulation results
t, x, v = np.array(history).T

# Calculate the analytical result at the same times
x_analytical = (
    initial_position
    + initial_velocity * t
    + 0.5 * acceleration * t**2
)

# Compare the results
plt.plot(t, x, label="Displacement (Euler)")
plt.plot(t, x_analytical, linestyle="--", label="Displacement (analytical)")

plt.xlabel("Time (s)")
plt.ylabel("Displacement (m)")
plt.legend()
plt.show()

If you run this with different, values for time-steps (dt) you will get the following results shown in Figure 1. Notice how the accuracy of the result depends on the value of dt.

Figure 1: Euler with different values of dt

2.3 Beyond Euler

I am introducing Euler’s method first because it allows us ‘get a feel’ for the system we are studying by building a simulation one step at a time. Once the model makes sense and we want greater accuracy and efficiency, we need alternatives.

Standard Python packages such as SciPy provide highly optimised, higher-order numerical solvers that you can unleash on your problems. You encountered one of these, odeint(), in 73. Let me now introduce you to its more modern cousin: solve_ivp().

Let me demonstrate it using our toy box problem.

solve_ivp() expects us to provide:

  1. the differential equations
  2. the initial conditions
  3. the time interval over which we want a solution

We express these equations as a Python function.

First, lets load the solver.

from scipy.integrate import solve_ivp

The differential equations

force = 10.0  # N
mass = 2.0    # kg

initial_position = 0.0  # m
initial_velocity = 0.0  # m/s

def box_model_derivatives(time, state):
    position, velocity = state

    acceleration = force / mass

    return [velocity, acceleration]

The first returned value is \(dx/dt\), while the second is \(dv/dt\).

Time interval of the simulation

I would like the simulation to stop when the box reaches \(1\,\mathrm{km}\). We can define this as an event reach_one_kilometre() which must return a zero at the event.

def reach_one_kilometre(time, state):
    position, velocity = state
    return position - 1000                

reach_one_kilometre.terminal = True       # stop when the event triggers
reach_one_kilometre.direction = 1         # trigger when zero is crossed -'ve to +'ve

Initial conditions

Now we can solve the equations:

solution = solve_ivp(
    fun=box_model_derivatives,
    t_span=(0, 100),
    y0=[initial_position, initial_velocity],
    events=reach_one_kilometre,
    max_step=0.1
)

t = solution.t           # time values
x, v = solution.y        # x, v values

Unlike our Euler code, solve_ivp() chooses its time steps adaptively. I.e. it takes smaller steps when the solution needs more care and larger steps when the solution is changing smoothly. Here, max_step limits the largest step it may take; it does not force every step to have that size.

Since we know the analytical solution, we can check the results. This is shown in Figure 2

Figure 2: Comparing solve_ivp() and the analytical solution.
WarningNot a panacea

Numerical solvers such as solve_ivp() are sophisticated, accurate, and efficient, but they are not infallible.

For instance, most general-purpose solvers are not designed to conserve energy even when the numerical error at each step is small. There is inevitable some drift.

So, we must always check whether a numerical solution behaves physically.

When we want energy conservation, there is a family of methods called symplectic integrators. These methods do not conserve energy exactly, but they usually keep the energy error bounded (i.e. varying about a mean) rather than allowing it to drift continuously.

Let me show you one of the simplest: symplectic Euler, also called the Euler–Cromer method.

2.4 A symplectic Euler model

Euler–Cromer makes one tiny but crucial change to forward Euler:

  1. calculate the new velocity using the old position
  2. calculate the new position using the new velocity

That’s it!

NoteSymplectic Euler (Euler-Cromer) 
acceleration = acceleration_at(position)

velocity += acceleration * dt
position += velocity * dt
NoteForward Euler 
acceleration = acceleration_at(position)

position += velocity * dt
velocity += acceleration * dt

Notice the order: I update the velocity first and then use that new velocity to update the position. This small change generally produces much more sensible long-term behaviour for an oscillating system.

3 Finding the period

Unlike an analytical function, a numerical solution contains only a finite set of points. It is unlikely that a peak or zero crossing will always fall exactly one of the points in this finite set of numbers. So, before we dig into the pendulum simulation, I want to show you two strategies for finding the period of a numerical solution.

  1. writing our own zero-crossing function
  2. using find_peaks() from SciPy

To see how they work, let us use a simple sine curve, whose period we know is \(T=2\pi = 360^\circ\).

3.1 Using zero crossings

One way to determine the period is to find successive events where the signal crosses zero in the same direction.

The following function finds positive-to-negative crossings by interpolation between the two points as seen in Figure 3. \[ \f{y_1}{x_0-x_1} = \f{-y_2}{x_2-x_0} \Rightarrow x_0 = \f{y_1x_2-y_2x_1}{y_1-y_2} \tag{9}\]

Figure 3: We approximate the zero-crossing point by interpolating along a straight line joining the adjoining points. 
def downward_zero_crossings(x, y):
    crossing_x = []

    for i in range(len(y) - 1):
        if y[i] > 0 and y[i + 1] <= 0:
            # Interpolate between the two surrounding points
            crossing = y[i]*x[i+1]-y[i+1]*x[i]
            crossing /= y[i ] - y[i+1]
            crossing_x.append(crossing)
    return np.array(crossing_x)

We can now estimate the period:

crossings = downward_zero_crossings(x, y)
periods = np.diff(crossings)

period_mean = periods.mean()
period_std = periods.std(ddof=1)

print(f"Expected period: {np.rad2deg(2 * np.pi):.6f}")
print(f"Mean period: {np.rad2deg(period_mean): .6f}")
print(f"Standard deviation: {np.rad2deg(period_std):.2e}")
Expected period: 360.000000
Mean period:  360.000000
Standard deviation: 3.17e-06

3.2 Using find_peaks()

SciPy’s find_peaks() function is more convenient approach because it identifies the local maxima of the signal:

from scipy.signal import find_peaks

peaks, properties = find_peaks(y)

peak_positions = x[peaks]
periods = np.diff(peak_positions)

period_mean = periods.mean()
period_std = periods.std(ddof=1)

print(f"Expected period: {np.rad2deg(2 * np.pi):.6f}")
print(f"Mean period: {np.rad2deg(period_mean): .6f}")
print(f"Standard deviation: {np.rad2deg(period_std):.2e}")
Expected period: 360.000000
Mean period:  359.965883
Standard deviation: 2.96e-01

We can check what find_peaks() detected by adding the peaks to the plot:

plt.plot(x, y)
plt.scatter(peak_positions, y[peaks], color="red", label="Detected peaks")
plt.xlabel("$x$")
plt.ylabel(r"$\sin x$")
plt.legend(loc='lower right')
plt.show()
Figure 4: The \(\sin \theta\) curve with the peaks located using find_peaks()

Note:

  • For the pendulum, replace x with the simulation time and y with the simulated angle.
  • Whichever method you use, estimate the period across several cycles by running the simulation for longer.
  • Although find_peaks() performs less well for this particular dataset, it offers useful options for dealing with noisy data. So, do not rush to discard it.

4 Simulating the pendulum; your turn!

It is now your turn to simulate the pendulum. We will proceed in stages so that you can see what works, what fails, and how to improve it.

4.1 Build the simulation

  1. Implement Equation 4 using forward Euler.
    Use:

    • \(L=1\) m
    • \(g\) from scipy.constants (from scipy.constants import g)
    • \(\theta(0)=10^\circ\)
    • \(\omega(0)=0\)
    • a simulation time of at least \(10\,\mathrm{s}\)

    Remember to convert the starting angle to radians before using np.sin().

  2. Plot \(\theta\) against time. What happens to the amplitude? Is this physically sensible?

  3. Fix the obvious problem by replacing forward Euler with the Euler–Cromer method. Plot the result again and compare the two simulations.

  4. Determine the period from your Euler–Cromer simulation. Measure across several oscillations rather than relying on only one.

  5. Repeat the simulation using progressively smaller values of \(\Delta t\). Continue until the period is stable to a precision that you consider adequate.

TipCheck for convergence

A numerical result is not trustworthy simply because the code runs. If changing \(\Delta t\) changes your answer appreciably, the calculation has not yet converged.

4.2 Investigate the small-angle approximation

Use your converged Euler–Cromer simulation to answer the following:

  1. Determine the period for a starting angle of \(10^\circ\).

  2. Compare it with the small-angle prediction

    \[ T_0=2\pi\sqrt{\frac{L}{g}}. \tag{10}\]

  3. Temporarily replace np.sin(angle) \(\longrightarrow\) angle and determine the period again.

  4. Compare the three results:

    • the nonlinear simulation using \(\sin\theta\)
    • the numerical small-angle simulation using \(\theta\)
    • the analytical small-angle prediction, \(T_0\)

    How far off is the small-angle approximation? Report both the absolute and percentage differences.

4.3 Compare \(T_{10^\circ}\) and \(T_{20^\circ}\)

Reset the code to the original nonlinear model containing \(\sin\theta\).

  1. Determine \(T_{10^\circ}\) and \(T_{20^\circ}\) using the same converged value of \(\Delta t\).

  2. Calculate \[ T_{20^\circ}-T_{10^\circ}. \]

  3. Express this difference as a percentage of \(T_{10^\circ}\).

  4. Compare your numerical results with the first analytical correction:

    \[ T(\theta_0) \approx T_0\left(1+\frac{\theta_0^2}{16}\right), \tag{11}\]

    Remember \(\theta_0\) must be expressed in radians.

  5. Finally, compare the predicted difference with the uncertainty in your experimental measurements. Was your experiment precise enough to detect it?

Back to top