Case Study: Numerical Analysis — Pendulum Clock¶
This case study illustrates how to mix gvar.GVar
s with numerical routines
for integrating differential equations (gvar.ode
) and for
finding roots of functions (gvar.root
). It also gives a
simple example of a simulation that uses gvar.GVar
s.
The Problem¶
The precision of a particular pendulum clock is limited by two dominant factors: 1) the length of the pendulum (0.25m) can be adjusted with a precision of at best ±0.5mm; and 2) irregularities in the drive mechanism mean that the maximum angle of swing (π/6) is uncertain by ±0.025 radians. The challenge is to determine how these uncertainties affect time-keeping over a day.
The angle theta(t)
of the pendulum satisfies a differential equation
d/dt d/dt theta(t) = -(g/l) sin(theta(t))
where g
is the acceleration due to gravity and the l
is the length
of the pendulum.
Pendulum Dynamics; Finding the Period¶
We start by designing code to integrate the differential equation:
import numpy as np
import gvar as gv
def make_pendulum(theta0, l):
""" Create pendulum solution y(t) = [theta(t), d/dt theta(t)].
Initial conditions are y(0) = [theta0, 0]. Parameter l is the
length of the pendulum.
"""
g_l = 9.8 / l
def deriv(t, y):
""" Calculate d/dt [theta(t), d/dt theta(t)]. """
theta, dtheta_dt = y
return np.array([dtheta_dt, - g_l * gv.sin(theta)])
y0 = np.array([theta0, 0.0])
return gv.ode.Integrator(deriv=deriv).solution(0.0, y0)
Given a solution y(t)
of the differential equation from this method,
we find the period of oscillation using gvar.root
: the period
is the time at which the pendulum returns to its starting point and its
velocity (y(t)[1]
) vanishes:
def find_period(y, Tapprox):
""" Find oscillation period of pendulum solution y(t).
Parameter Tapprox is the approximate period. The code finds the time
between 0.7 * Tapprox and 1.3 * Tapprox where y(t)[1] = d/dt theta(t)
vanishes. This is the period, provided Tapprox is correctly chosen.
"""
def dtheta_dt(t):
""" vanishes when dtheta/dt = 0 """
return y(t)[1]
return gv.root.refine(dtheta_dt, (0.7 * Tapprox, 1.3 * Tapprox))
Analysis¶
The last piece of the code does the analysis:
def main():
l = gv.gvar(0.25, 0.0005) # length of pendulum
theta_max = gv.gvar(np.pi / 6, 0.025) # max angle of swing
y = make_pendulum(theta_max, l) # y(t) = [theta(t), d/dt theta(t)]
# period in sec
T = find_period(y, Tapprox=1.0)
print('period T = {} sec'.format(T))
# uncertainty in minutes per day
fmt = 'uncertainty = {:.2f} min/day\n'
print(fmt.format((T.sdev / T.mean) * 60. * 24.))
# error budget for T
inputs = dict(l=l, theta_max=theta_max)
outputs = dict(T=T)
print(gv.fmt_errorbudget(outputs=outputs, inputs=inputs))
if __name__ == '__main__':
main()
Here both the length of the pendulum and the maximum angle of swing
have uncertainties and are represented by gvar.GVar
objects. These uncertainties
work their way through both the integration and root finding to give
a final result for the period that is also a gvar.GVar
. Running the code
results in the following output:
period T = 1.0210(20) sec
uncertainty = 2.79 min/day
Partial % Errors:
T
--------------------
l: 0.10
theta_max: 0.17
--------------------
total: 0.19
The period is T = 1.0210(20) sec
, which has an uncertainty of
about ±0.2%. This corresponds to an uncertainty of ±2.8 min/day
for the clock.
The uncertainty in the period is caused by the uncertainties in the
length l
and the angle of maximum swing theta_max
.
The error budget at the end of the output shows how much error comes
from each source: 0.17% comes from the angle, and 0.10% comes from
the length. (The two errors added in quadrature give the total.)
We could have estimated the error due to the length from the
standard formula 2π sqrt(l/g) for the period, which is approximately
true here. Estimating the uncertainty due to the angle is trickier, since it
comes from nonlinearities in the differential equation.
The error budget tells us how to improve the clock. For example, we can reduce the error due to the angle by redesigning the clock so that the maximum angle of swing is π/36 ± 0.025 rather than π/6 ± 0.025. The period becomes independent of the maximum angle as that angle vanishes, and so becomes less sensitive to uncertainties in it. Taking the smaller angle reduces that part of the period’s error from 0.17% to 0.03%, thereby cutting the total error almost in half, to ±0.10% or about ±1.5 min/day. Further improvement requires tighter control over the length of the pendulum.
Simulation¶
We can check the error propagation analysis above using
a simulation. Adding the following code at the end of main()
above
# check errors in T using a simulation
Tlist = []
for i in range(100):
y = make_pendulum(theta_max(), l())
T = find_period(y, Tapprox=1.0)
Tlist.append(T)
print('period T = {:.4f} +- {:.4f}'.format(np.mean(Tlist), np.std(Tlist)))
gives the following additional output:
period T = 1.0209 +- 0.0020
The new code generates 100 different values for the period T
, corresponding
to randomly chosen values for theta_max
and l
drawn from the
Gaussian distributions corresponding to their gvar.GVar
s. (In general, each
call x()
for gvar.GVar
x
is a new random number drawn from x
’s
Gaussian distribution.)
The mean and
standard deviation of the list of periods give us our final result.
Results fluctuate with only 100 samples; taking 10,000 samples shows that
the result is 1.0210(20), as we obtained
in the previous section above (using a tiny fraction of the computer time).
Note that the gvar.GVar
s in this simulation are uncorrelated and so their random
values can be generated independently. gvar.raniter()
should be used to
generate random values from correlated gvar.GVar
s.