Asking difficult questions with a numerical model

Activity 1 | Week 3 | IS

1 What to expect

So far, we have developed a model of the simple pendulum and solved it using both the small-angle approximation and numerical integration. Our investigative journey is summarised in Figure 1.

flowchart LR
    A["Measure<br/>T₁₀ and T₂₀"]
    B["Look for a<br/>difference"]
    C["Construct a<br/>physical model"]
    D["Small-angle<br/>approximation"]
    E["Numerical<br/>model"]
    F["Forward Euler<br/>(simple, but unstable)"]
    G["Euler–Cromer<br/>(better behaved)"]

    A --> B --> C
    C --> D
    C --> E
    E --> F --> G

    classDef experiment fill:#e8f1fb,stroke:#27628a,stroke-width:2px,color:#111;
    classDef theory fill:#f5ead7,stroke:#a26721,stroke-width:2px,color:#111;
    classDef computation fill:#e5f3e7,stroke:#347a3d,stroke-width:2px,color:#111;

    class A,B experiment;
    class C,D theory;
    class E,F,G computation;
Figure 1: Our journey from measurement to modelling

What we have done helps us understand the difference between \(T_{10^\circ}\) and \(T_{20^\circ}\). However, the model has not yet given us much insight into how we could improve the experiment.

That is what I want to focus on today. I want us to make the model more realistic, compare it with evidence, and then use it to investigate questions that would be difficult to isolate experimentally.

ImportantThe aim

Make the model realistic enough to give us useful insight; and then use it to ask questions that are difficult, slow, or impractical to answer experimentally.

This is one of the most powerful uses of computation in science. A model allows us to explore situations beyond our immediate experimental reach. For example:

  • AlphaFold predicts three-dimensional protein structures
  • SRIM models how ions travel through matter
  • CASINO simulates electron trajectories in solids
  • Spartan predicts molecular structures and properties
  • COMSOL

These models do not replace experiments. They help us decide what to expect, what to measure, and which experiments are worth doing.

1.1 The questions

Once we have a more realistic model, we can investigate questions such as:

  1. Does damping change the period?
    Is the change large enough to matter in our experiment?

  2. How many oscillations should we time?
    Does timing more oscillations always improve the result, or does amplitude decay eventually become important?

  3. What happens if we accidentally push the bob?
    How does a non-zero initial angular velocity affect the period we measure?

2 Improving our model

There are many ways to improve the pendulum model, with varying levels of sophistication. Let us start simply and take the model one step closer to reality by adding damping.

Damping allows the simulated pendulum to slow down, as a real pendulum does. It also gives us another opportunity to compare a model with data.

# Model Equation What it captures What it leaves out
1 Small-angle, undamped \(\daa{\theta}{t}=-(g/L)\theta\) Nearly sinusoidal motion at small amplitude Amplitude dependence and energy loss
2 Nonlinear, undamped \(\daa{\theta}{t}=-(g/L)\sin\theta\) Amplitude dependence of the period Air resistance and friction
3 Nonlinear with linear damping \(\daa{\theta}{t}=-(g/L)\sin\theta-\gamma\da{\theta}{t}\) Amplitude dependence and gradual energy loss Dry friction, nonlinear drag, and out-of-plane motion
4 Linear damping and dry friction \(\daa{\theta}{t}=-(g/L)\sin\theta-\gamma\da{\theta}{t}-\alpha_{\mathrm f}\operatorname{sgn}\left(\da{\theta}{t}\right)\) Velocity-dependent damping and constant-magnitude friction Static friction, nonlinear drag, and out-of-plane motion
Table 1: Model ladder for the simple pendulum.

We will focus on Model 3:

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

Here, \(\gamma\) controls how quickly the oscillations decay. Our first challenge is to find a value of \(\gamma\) that gives a reasonable description of your real pendulum.

One neat approach would be for every group to record a video, track the bob, and extract \(\gamma\) independently. As enjoyable1 as that might be, we do not have enough time today. Instead, we will use a common dataset so that we can quickly get the model working and focus on what it can teach us.

I will show you how to:

  1. add damping to the simulation
  2. compare the model with the common data
  3. estimate a value of \(\gamma\) that matches the data from your experiment
  4. use the calibrated model to investigate an experimental question

3 What you have to do

3.1 Step 1: Rewrite your old code

Rewrite your old simulation as a function that can be called using different parameter values.

Then modify it so that it:

  • detects downward zero crossings of angular speed \(\omega\) internally.
  • records the interpolated crossing times
  • stops after a specified number of oscillations

This mirrors what you did experimentally: instead of running for an arbitrary time, you measured a chosen number of complete oscillations.

A useful function structure is:

def simulate_pendulum(
    amplitude_degrees,
    number_of_oscillations,
    dt=0.01,
    length=1,
    gamma=0
):
   period_info = []
   oscillations = 0
   while oscillations < number_of_oscillations:
   ...
   ...
   ...

   period_info = np.diff(period_info)
   return time, angle, speed, periods

The basic strategy is:

  1. advance the simulation by one Euler–Cromer step
  2. check whether \(\theta\) has crossed zero in the downward direction—from positive to negative
  3. use linear interpolation to estimate the crossing time
  4. store the crossing time
  5. stop when you have recorded enough crossings

In pseudocode:

Create an empty list of crossing times

While fewer than N crossings have been recorded:

    Save the old time and angle

    Advance the simulation by one Euler–Cromer step

    If the angular speed changes from positive to negative:

        Interpolate the zero-crossing time

        Add the crossing time to the list

Calculate the individual periods from the differences
between successive crossing times

For two points \((t_\text{old},\omega_\text{old})\) and \((t_\text{new},\omega_\text{new})\) on opposite sides of zero, the interpolated crossing time is

\[ t_{\mathrm{cross}} = \frac{\omega_\text{old}t_\text{new}-\omega_\text{new}t_\text{old}} {\omega_\text{old}-\omega_\text{new}}. \tag{2}\]

After the simulation, calculate the individual periods using

periods = np.diff(crossing_times)

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

Here is what I obtained with my code.

Figure 2

Step 2: Add damping to the model

It is almost absurdly easy to incorporate damping: we only need to change the angular-acceleration term.

Our undamped model is

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

I will model the damping as being proportional to the angular velocity. The damped model is

\[ \frac{d^2\theta}{dt^2} = -\frac{g}{L}\sin\theta -\gamma\frac{d\theta}{dt}, \tag{3}\]

where \(\gamma\) is the damping rate and has units of \(\mathrm{s^{-1}}\). Setting \(\gamma=0\) recovers our original undamped model.

Letting

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

we can write Equation 3 as two first-order equations:

\[ \begin{aligned} \frac{d\omega}{dt} &= -\frac{g}{L}\sin\theta-\gamma\omega,\\ \frac{d\theta}{dt} &= \omega. \end{aligned} \]{#eq-damped-pendulum-first-order

For velocity-first Euler–Cromer,

\[ \begin{aligned} \omega_{n+1} &= \omega_n+ \left[ -\frac{g}{L}\sin(\theta_n) -\gamma\omega_n \right]\Delta t,\\ \theta_{n+1} &= \theta_n+\omega_{n+1}\Delta t. \end{aligned} \tag{4}\]

In Python, the only change is:

angular_acceleration = (-(g / length) * np.sin(angle)- gamma * angular_velocity)

The Euler–Cromer updates remain:

angular_velocity += angular_acceleration * dt
angle += angular_velocity * dt

Start with \(\gamma=0.2\,\mathrm{s^{-1}}\) and compare the result with \(\gamma=0\). Here is what I obtained with my code for 20 oscillations.

Figure 3: Comparison of damped and undamped pendulum motion over 20 oscillations.

3.2 Estimate \(\gamma\) for your experiment

Use trial and error to estimate \(\gamma\) from the amplitude decay in your \(10^\circ\) and \(20^\circ\) experiments.

For each starting angle:

  1. choose a value of \(\gamma\)
  2. run the simulation for the same number of oscillations as your experiment
  3. compare the simulated amplitude decay with your observations
  4. adjust \(\gamma\) until you obtain a reasonable match

There are more systematic approaches, such as bisection and optimisation, but I want to keep the complications to a minimum for now.

Ideally, one value of \(\gamma\) should describe both starting angles. If you obtain substantially different values, this may indicate that:

  • the measurements are uncertain
  • the damping is not proportional to angular velocity
  • our model is missing another effect, such as dry friction

This is not a failure but useful information for our next modelling iteration. However, lets leave this for another time.

3.3 Period over successive oscillations

Use your estimated value of \(\gamma\) to simulate 30 oscillations starting from \(10^\circ\) and \(20^\circ\).

Plot the period of each oscillation against the oscillation number. Then ask:

  1. Does the period remain constant?
  2. Does it approach a limiting value as the amplitude decreases?
  3. Is the change large enough to matter experimentally?
  4. How do the results differ between starting angles of \(10^\circ\) and \(20^\circ\)?

4 Last task: How many oscillations should we have timed?

4.1 Timing uncertainty

In Week 1, we asked how many oscillations should be included in one timing measurement. We can now return to that question with a better understanding of the pendulum.

Timing more oscillations reduces the effect of reaction time. Suppose the uncertainty associated with each press of the stopwatch is \(\delta t\). Since we press it once to start and once to stop, the uncertainty in the total measured time is approximately

\[ \delta t_{\mathrm{total}} \approx \sqrt{2}\,\delta t. \tag{5}\]

If we time \(N\) oscillations, the uncertainty in the period is therefore

\[ \delta T \approx \frac{\sqrt{2}\,\delta t}{N}. \tag{6}\]

However, timing more oscillations also means that:

  • the amplitude decreases
  • the nonlinear period changes as the amplitude decreases
  • out-of-plane motion and other imperfections may accumulate
  • the motion eventually becomes difficult to observe

There is a trade-off: more oscillations reduce timing uncertainty, but give the pendulum more time to change.

Comparing two periods

We measured \(T_{10^\circ}\) and \(T_{20^\circ}\) independently. The approximate uncertainty in their difference is therefore

\[ \delta\left(T_{20^\circ}-T_{10^\circ}\right) = \sqrt{ \left(\delta T_{20^\circ}\right)^2 + \left(\delta T_{10^\circ}\right)^2 }. \]

If the two period measurements have similar uncertainties, then

\[ \delta\left(T_{20^\circ}-T_{10^\circ}\right) \approx \frac{2\delta t}{N}. \tag{7}\]

4.2 What you have to do:

  1. Use \(\delta t=0.2\,\mathrm{s}\) in Equation 12 and stimate the uncertainty in one period measurement for \[ N=1,\ 5,\ 10,\ 20,\ \text{and}\ 30. \]

  2. Use Equation 7 to estimate the uncertainty in \(T_{20^\circ}-T_{10^\circ}\).

  3. Compare these uncertainties with the difference predicted by your model.

  4. Use the damped model to determine how much the amplitude and period change over the same numbers of oscillations.

  5. Recommend how many oscillations you would time if you repeated the Week 1 experiment.

Your recommendation should balance:

  • timing precision
  • amplitude decay
  • the changing period
  • the practical difficulty of observing the motion
Back to top

Footnotes

  1. Or not.↩︎