Objectives

The simple pendulum

Equation of motion

We can deduce the equation of motion of the pendulum from the law of conservation of the total mechanical energy.

This energy \[E_m = E_c + E_p \]

is the sum of a kinetic energy and a potential energy due to the gravitational field exerted by the Earth on the mass of the pendulum.

If \(\theta\) is the angle of the pendulum with the vertical,

and \(\dot{\theta} = \frac{d\theta}{dt}\), the kinetic energy is

\[E_c = \frac{1}{2}ml^2 \dot{\theta}^2\]

This energy becomes null at the maximum amplitude reached \(\theta_{max}\). The potential energy is

\[E_p = -mgl \cos \theta\]

At any time, we have:

\[\frac{1}{2}ml^2 \dot{\theta}^2 -mgl\cos \theta = -mgl\cos \theta_{max} \] Hence

\[ \dot{\theta}^2 = \frac{2g}{l} (\cos \theta -\cos \theta_{max} ) \]

Harmonic pendulum limit

If we assume that \(\theta\) is small, we can use the second order approximation:

\[\cos(\theta ) = 1 - \frac{1}{2} \theta^2 + \ldots\]

This gives:

\[ \dot{\theta}^2 = - \frac{g}{l} \theta^2 + C \]

This shows that the time scale of the motion is \(\sqrt{\frac{g}{l}}\)

Indeed changing variables lead to:

\[\dot{x}^2 = 1 - x^2 \]

whose solution is the harmonic motion \(x(t)= \cos ( \omega(t-t_0))\)

Alternatively, writing the balance of forces (Newton’s third law) gives:

\[ \frac{d^2 \theta}{dt^2} + \frac{g}{l} \sin \theta = 0 \]

And the small angle approximation

\[ \sin \theta \approxeq \theta \]

gives the harmonic pendulum motion equation:

\[ \frac{d^2 \theta}{dt^2} + \frac{g}{l} \theta = 0 \]

whose solution is \[\theta(t) = \theta_{max} \cos( 2 \pi (t-t_{max}) / T_0)\]

This is a harmonic motion with period \(T_0=2 \pi \sqrt{ \frac{l}{g}}\)

The italian physicist Galileo Galilei was the first to note that this period is the constant regardless of the mass attached, and depends only on the length of the pendulum (and g).

Nonharmonic pendulum

If we do not want to make any approximation, and accept to displace the pendulum from any angle from the rest position, we must keep the equation:

\[ \dot{\theta}^2 = \frac{2g}{l} (\cos \theta -\cos \theta_{max} ) \]

We can separate time and angles

\[ dt = \sqrt{\frac{l}{2g}} \frac{d\theta}{ \sqrt{\cos \theta -\cos \theta_{max} }} \]

and integrate over one period to obtain:

\[T(\theta_{max}) = \frac{2T_0}{\pi} K(k) \]

with \[k = \sin (\frac{\theta_{max}}{2}) \]

and \(K(k)\) is the elliptic integral of the first kind:

\[K(k) = \int_0^{\frac{\pi}{2}} \frac{ds}{\sqrt{1-k^2 \sin^2(s) }}\]

Approximate period solution by Borda

This elliptic integral cannot be integrated exactly in terms of simple functions.

But we can again try to approximate this integral. Consider the behavior if \(k\) is small.

\[(1-k^2 \sin^2(s) )^{-1/2} \approxeq 1+ \frac{1}{2}k^2 \sin^2(s) + \ldots \]

Hence, by integration, we get the following approximation

\[T = T_0 (1+ \frac{1}{4} k^2 ) \]

As for \(k\) small, we have $k = /2 $, we get the formula found by Borda in :

\[T = T_0 (1+ \frac{1}{16} \theta_{max}^2 + \ldots ) \]

This formula is an improvement over the harmonic approximation as it adds one term that depends on the maximum amplitude. Still, this formula is an another approximation and will fail and become inaccurate for high amplitudes.

Exact period solution : calculation

Let us recall the elliptic integral:

\[I(a,b)=\frac{2}{\pi} \int_0^{\pi/2} \frac{d\phi}{\sqrt{a^2 \cos^2 \phi + b^2 \sin^2 \phi}} \]

The value of this integral is easily shown to be invariant by the transformation \[T(a,b) = ((a+b)/2,\sqrt{ab})\]. This gives:

\[ I(a,b) = I(T(a,b))=I(a_1,b_1) = \ldots \]

Hence

\[ I(a,b) = I(T^n(a,b)) = I(a_n,b_n) \ldots = I(a_{\infty},b_{\infty}) \]

and

\[ I(a,b) = \frac{1}{l(a,b)} \] where \[l(a,b)\] is the limit of the so called arithmetico-geometric suite \[(a_n,b_n)=T^n(a,b)\] when \(n\) tends to infinity.

Hence, to calculate \(I(a,b)\) we just need to find the limit \(l(a,b)\), a property that we will use in our R code below.

Note that \[K(k) = I(1,\sqrt{1-k^2}) \], therefore, we have found a way to calculate the anharmonic period up to any desired precision.

Now, we can Write a function in R to use this theory. It should become clear to you what this function does, even if you have not entirely followed the mathematical details above. All we do here is calculate the limit of the arithmetico-geometric sequence for a properly chosen set of initial values \((a,b)\), in relation to the elliptic integral given above.

# g is the local acceleration of gravity
# l is the length of the pendulum
# thetamax is the maximum amplitude

mypendulum <- function(g,l,thetamax){

# this is the harmonic period
T0=2*pi*sqrt(l/g)
T=T0
# we initialize the AG suite to (a,b)
a=1
b=cos(thetamax/2)
#print(b)

# we will iterate as long as we have not reached this preset precision
tole=0.00001
notpreciseenough = TRUE
while( notpreciseenough ){
  
# the next three lines apply the transform
# that defines the next pair (a,b) of the suite
c = (a+b)/2
b=sqrt(a*b)
a=c
Tprec = T
#print(c(a,b,c,Tprec,T0))
# we update the value of T
T = T0/a
# and calculate the change in value
e=abs(T-Tprec)
#print(e)
# here we reset the boolean if we have reached the 
# required precision
# this will allow us to exit the while loop
if(e < tole){notpreciseenough = FALSE}
}
#print(c(e,a,b))
# we return a list of two elements
# the harmonic and anharmonic period
return(list(T0=T0,T=T))
}

Try the function

g=9.81
l=1 
# a pendulum of 1 meter has a period of 2 seconds
# for 9.81 = g
# 2.0060666807106475
thetamax = pi/100
rm(T,T0)
## Warning in rm(T, T0): object 'T' not found
## Warning in rm(T, T0): object 'T0' not found
r   = mypendulum(g,l,thetamax)
T0  = r[[1]]
T   = r[[2]]
paste("harmonic period: ",T0)
## [1] "harmonic period:  2.00606668071065"
paste("anharmonic period: ",T)
## [1] "anharmonic period:  2.00619043198655"

Now plot the period as a function of max amplitude

g=9.81
l=1
# a pendulum of 1 meter has a period of 2 seconds
# for 9.81 = g
# 2.0060666807106475

ma = (pi/2) * 1

tv <- seq(0, ma, by = ma/100 )
out = c()

thetamax = pi/100

for(tmax in tv){
rm(T,T0)
thetamax = tmax
r   = mypendulum(g,l,thetamax)
T0  = r[[1]]
T   = r[[2]]
Trel = (T-T0)/T0
out = c(out,Trel)
#paste("harmonic period: ",T0)
#paste("anharmonic period: ",T)

}

plot(tv,out,main="Relative error on oscillation period")

#out

paste('maximum relative error is ',max(out), '%')
## [1] "maximum relative error is  0.180340599016096 %"

Exercise: write your own function to see how the approximation of Borda fails when the amplitude becomes large. Overlay the plot of the Borda period and the plot of the exact period to visualize their difference.

Try stopping after two iterations in the AG sequence. How accurate is the approximation:

\[T_2 = \frac{T}{a_2}\]

Can you think of other sources of errors in the expression of the period of the pendulum?

Richer in Cayenne (1672)

Richer, who had set the pendulum of his astronomical clock leaving Paris, to beat the second, found that to correctly follow the daytime movement in Cayenne, he had to shorten his pendulum by 1 line a quarter; which involved a reduction in gravity g under the formula giving the period, and proved that Cayenne was farther from the center of the Earth than Paris, under the law of attraction. Finally, another presumption, if we consider a fluid sphere in rotation, the centrifugal force none at the pole will not be at the equator and will tend to spread the molecules of the fluid to the outside, until the attraction of the fluid counterbalances this action, that’s why Jupiter and Saturn, whose diameter had been precisely observed by Cassini, are flattened.

A clock set at the equator advances 226 ″ per day to the pole. If we denote by w the angular speed of the Earth’s rotation, and by m the w2 R / g ratio of the equatorial axial acceleration to that of gravity, astronomical measurements combined with the result of Picard, made it possible to evaluate this number at 1/288.

Huygens, considering a spherical earth in rotation, wonders what shape would take a subject surface fluid: on the one hand to the central attraction of this earth, on the other hand by axifuge acceleration; composing the two vectors, it shows that their result is normal to a 1/576 flattening spheroid. In his Principia Mathematica Philosophia naturalis, Newton studies the combined effect of universal attraction and centrifugal force on the shape of the Earth. It shows that a fluid and homogeneous Earth of revolution with equatorial symmetry would take the form of a flattening spheroid equal to 1/230, or 5/4 m.

When we read the Principia, we are struck by Newton’s insistence on justifying its result through experience. He knew the values ​​measured in Paris, London, Copenhagen, by Picard and others, he knew that Richer had found for the pendulum to second a length of 439.25 lines and had shown that on a flattened ellipsoid the gravimetric flattening was in inverse proportion to the ratio of the axes: the result of Richer and some other physicists at low latitudes therefore confirmed his views.

Adopting the average value of the Picard radius and the pendulum measures then known, it gives a table of the lengths of the pendulum beating the second L = Lo (1 + b sin2 l) and also gives the arc of a degree of amplitude on its 1/230 flattening spheroid (l is latitude).

The universal attraction shed such light on planetary dynamics that it should have won the support of the geodesists on the shape of the earth. This was not the case: a violent fight opposed the supporters of the flattening and those of the oblong shape during the fifteen years which followed the publication in 1723 from the meridian of Jacques Cassini. Uncertainty will be lifted in 1737 by the mission of Lapland.

When Newton wanted to demonstrate, in 1666, that the Earth attracts the Moon by the same force of gravity which is exerted on the bodies in free fall, it needed the length of the terrestrial degree; he used Wright’s value and found a 15% difference between the calculation and the measurements, which was unacceptable. Resuming the calculation, in 1685, with the value of Picard for the terrestrial radius and those of Huygens and Richer for gravity, he obtained the experimental proof of universal gravity.

Role of latitude

Correction for latitude:

\[g(\phi) = g_0 *( 1 + a_0 \sin(\phi)^2 + a_2 \sin(2\phi)^2) \]

#in ms^-2

# the numbers used depend on a regression model

# model 1:
#g =9.780327 * (1+0.0053024sin2ϕ−0.0000058sin22ϕ)


lat =  pi/2

lat =  pi/4

sinlat = sin(lat)
sin2lat = sin(2*lat)

gbase = 9.780327

a0 =  1.0000000
a1 =  0.0053024
a2 = −0.0000058


g = gbase*(a0 + a1*sinlat^2 + a2 * sin2lat^2)
gmodel1 = g
gmodel1
## [1] 9.8062
# model 2: 1967. fitted to observed g data
# reference ellipsoid
gbase = 978.03185 # in gals
a0 = 1.
a1 = 0.005278895
a2 = 0.000023462


g = gbase*(a0 + a1*sinlat^2 + a2 * sinlat^4)
gmodel2 = g/100.

gmodel2
## [1] 9.806191
# north or south pole: 
#[1] 9.832186
#[1] 9.832177


# halifax NS : 9.80567
# see calculator https://www.wolframalpha.com/widgets/view.jsp?id=e856809e0d522d3153e2e7e8ec263bf2

-Exercise:

Build a data frame of cities with altitude and latitudes and acceleration of gravity

https://en.wikipedia.org/wiki/List_of_cities_by_elevation

Mexico: h=2660 m g = 9.779

-The value of gravity in Washington DC

https://stat.ethz.ch/R-manual/R-devel/library/boot/html/gravity.html

expressed as deviations from 980.000 in centimetres per second squared.

library(boot)
gravity

Now change the latitude and get the pendulum period anywhere on Earth.

g=9.81 # plug the local gravity value here
l=1
# a pendulum of 1 meter has a period of 2 seconds
# for 9.81 = g
# 2.0060666807106475
thetamax = pi/100
rm(T,T0)


r   = mypendulum(g,l,thetamax)
T0  = r[[1]]
T   = r[[2]]
paste("harmonic period: ",T0)
## [1] "harmonic period:  2.00606668071065"
paste("anharmonic period: ",T)
## [1] "anharmonic period:  2.00619043198655"

Role of height

Free air correction or Elevation correction

According to the law of Newton, any two masses \(m\) and \(m'\) attract each other with a gravitational force of the form

\[F = -G \frac{mm'}{r^2}\]

where G is a universal constant and \(r\) is the distance between the two masses.

Hence the Earth (mass \(M\)) attract any small mass \(m\) by a force that decreases quadratically ith the distance \(r\)$ between the center of the test mass \(m\) and the center of the Earth.

This implies that \(g=g(r)=g(r_0)\frac{r_0^2}{(r_0+h)^2}\), where \(r_0\) is the radius of the Earth (supposed to be spherical) and \(g_0\) is the reference gravity acceleration ‘at sea level’.

This can be approximated to first order by

\[g = g_0(1-\frac{2h}{r_0}) = g_0(1 + a_h h)\]

Hence, at elevation/altitude \(h\) with respect to sea level, the acceleration is modified by an amount \[\Delta_{fa}g = 2g_0/r_0\]

If we set the ‘sea level’ value at \(g_0 = 9.78\) and choose \(r_0=6378km\), this gives an elevation correction of \[ -3.07 * 10^{-7} m s^{-2} \] per meter of elevation

r0=6378000.
ah = -2/r0
ah # -3.10^{-7}
## [1] -3.135779e-07

The elevation correction needs to be added to the latitude correction:

h=0
#in ms^-2

# the numbers used depend on a regression model

# model 1:
#g =9.780327 * (1+0.0053024sin2ϕ−0.0000058sin22ϕ)
lat =  pi/2

sinlat = sin(lat)
sin2lat = sin(2*lat)

gbase = 9.780327

a0 =  1.0000000
a1 =  0.0053024      # 
a2 = −0.0000058
ah = -3.13 * 10^(-7) # h in meters

g = gbase*(a0 + a1*sinlat^2 + a2 * sin2lat^2 + ah * h)

gmodel3 = g
gmodel3
## [1] 9.832186

Role of Earth rotation

Role of the attraction of the Moon and the Sun

A look at the following figure

clearly shows that a distant planet (e.g. Moon or Sun) does exert the same gravitational potential at the point \(O\) (center of the Earth) and at a point located on the surface (\(P\)), due to the difference between the distance to the planet. This leads to tidal forces.

It can be shown that to second order, the tidal potential \(W\) (difference between the two potentials) is:

\[W = V_P - V_0 \approxeq W_2 = \frac{3}{4} f \mu \frac{r_0^2}{r^3} (cos(2z) + 1/3)\] where \(\mu = \frac{mM}{m+M}\) is the reduced mass of the planet,\(r_0=OP\), and \(f=GM\) is the geocentric gravitational constant, and that relation implies a variation of the vertical acceleration of gravity of magnitude given by:

\[\Delta_2 g = -\frac{\partial W_2}{\partial r_0 } = \frac{G m r_0 }{r^3}(3 \cos^2(z) - 1) \]

Can you calculate the magnitude of this anomaly in microgals, knowing that ?

massmoon  = 7.342 *  10^22  #  kg
# 1,23 % (123/10000) * massearth
massearth = 5.9722 * 10^24 # kg 
# mu = massmoon #  (massmoon*massearth)/(massmoon+massearth)
# mu
G = 6.67408 * 10^(-11)     #  m3 kg-1 s-2
#f=G*massearth # geocentric gravitational constant
# the distance Earth-Sun  is  20000 * r_0
# the distance Earth-Moon is     60 * r_0
#f           # 3.98 10^14  
r0=6378000.  # radius of Earth
r=60*r0      # Earth Moon distance
#r
c  = G*massmoon*r0*r^{-3}
c*3  #
## [1] 1.673032e-06
# 1 milligal is 10^{-5} ms^{-2} USI
dgm = c*3 * 10^5  # in mgals
dgsol = 10^3 * dgm
paste("the moon correction is ",dgsol,"microgals")
## [1] "the moon correction is  167.303205488353 microgals"

The formulas above are complicated. So let us calculate how the Moon modifies the attraction exerted on us via a very simple formula: we will just say that we are closer to the Moon.

Therefore

\[\Delta g = Gm / (r-r_0)^2 - Gm/r^2\]

dg = G*massmoon/(r-r0)^2 - G*massmoon/r^2
dg  # in ms^(-2)
## [1] 1.143871e-06

The solid pendulum

The formulas of the simple pendulum can become … vastly wrong!

Indeed the simple pendulum model assumes that the entire mass of the pendulum is concentrated in a point. This approximation is valid if the pendulum is composed of a small solid attached to a long thin rope/thread. However, when the length of the pendulum becomes comparable to the length of the massive body, the approximation does not hold. In fact, the formula for the period of the harmonic pendulum becomes very wrong when the length tends towards zero, because the period is predicted to be null. However, if the center of rotation of the pendulum and of the solid coincide, we expect the solid to not oscillate and stay at rest, i.e. an infinite period!

Christiaan Huyghens understood this perfectly and developed a fully functional theory of the solid pendulum. His theory led others, namely Prony and Kater, to a device that was the best tool for measuring gravity acceleration for 150 years, the Kater pendulum.

Kater pendulum

Up to the 1930s, the most precise measurements of earth gravity were obtained by a method designed by Henry Kater in 1818, which is itself based on an idea suggested by Prony, and, of course, Christiaan Huyghens.

The idea of this instrument is the concept of reversed pendulum, of which a typical example is shown in the following figure:

It can be shown, by an application of the laws discovered by Hyughens, that if the distances of the two masses are adjusted so that the two periods of oscillations of the pendulum and the pendulum ‘upside down’ coincide, the common period of oscillation should be the same as for a simple pendulum of equivalent length \(a+b\).

Applying the laws of Huyghens to the small oscillation amplitude limit of the Kater pendulum shown in figure gives:

\[\tau = 2 \pi \sqrt{\frac{a+b}{g}}\]

\[g = 4 \pi^2 \frac{a+b}{\tau^2}\]

  • Experiment 1

Instead of finding experimentally how to match the two periods, which involves tedious and time-consuming manipulation, one can fit a quadratic polynomial (we will see how later in this course) to each of the period-distance experimental points, and then, determine graphically the points of intersections of the two curves

\(T_1 = Q_1(x)\) and \(T_2 = Q_2(x)\).

https://www.phys.uconn.edu/~hamilton/phys258/N/k.pdf

Fit a quadratic polynomial

Student friendly Kater pendulum:

http://physics.mercer.edu/petepag/p_pend.html

-Experiment 2 : Kater cubic

It is of anecdotic interest to note that certain configurations of the Kater pendulum lead not to two but to three matching positions, as shown by this paper:

http://personalpages.to.infn.it/~zaninett/projects/pendolo.pdf

You can try to use their data and plot their experimental points in R:


x (cm) T1 (s) T2 (s)
10 2.3613 2.0615
20 2.1492 2.0337
30 2.0363 2.0089
35 2.0016 1.9999
40 1.9838 1.9931
45 1.9733 1.9911
50 1.9754 1.9894
55 1.9799 1.9908
58 1.9846 1.9924
65 2.0055 2.0002
68 2.0173 2.0064
75 2.0470 2.0273
85 2.0939 2.0678
90 2.1224 2.0969
92 2.1334 2.1071
106 2.2178 2.2174
110 2.2441 2.2589
120 2.3078 2.3776

Bouguer and the Chimborazo

36 pouces 7019 lignes

To find the period of his pendulum, Bouguer counts the number of oscillations in 24 hours. This will give the period in hours. Having the period, he can then deduct the ‘effective length’ of his pendulum at the location of the measure.

One of the measures takes place on top of the Pinchincha,

During the 1735-1745 French Academy of Sciences expedition to Peru, Pierre Bouguer conducted two experiments confirming Newtonian gravitational attraction and estimating the mean density of the Earth. One set of experiments determined the variation in gravity with altitude using a pendulum at sea level, Quito (2,860 m) and the summit of Pichincha (4,784 m). Qualitatively correct, Bouguer reported a smaller decrease in gravity than that predicted from altitude increase alone, but he calculated that the mean density of the Earth was nearly five times that of the near-surface rocks, an overestimate by a factor of at least two. The better reported experiment was an attempt to detect the deflection of the vertical near the mountain Chimborazo. There was a large difference between Bouguer’s predicted plumb-line deflection, 103“, and that which he and La Condamine measured, just 7”. I have investigated both experiments using a digital elevation model to compute the vertical and horizontal components of the gravity field caused by topography, and include the regional gravity signature of the Andean crustal root. The modelling indicates not only that Bouguer’s pendulum measurements were extremely accurate, but also that his observations allow a good determination of the significant isostatic effect. In contrast, on Chimborazo, contrary to recent suggestions, isostatic effects are negligible, and Bouguer’s deflection was, within error, in line with the modelled plumb-line deflection from topography. Both Bouguer’s pendulum and plumb-line measurements were reliable and therefore he should now be redeemed from any inference of failure.

de La Condamine C-M (1751) Journal du voyage fait par ordre du Roi a l’équateur, servant d’introduction historique a la mesure des trois premiers degres du méridien. Imprimerie Royale, Paris

Bouguer anomaly

Consider now the fingerprint of mountains and oceans on the acceleration of gravity.

On the map below, you see the anomaly of gravity (after removal of the average earth acceleration, the free air anomaly, the latitude effect, the effect of sun and moon). This residual gravity, called Bouguer anomaly shows an extraordinary strong signal.

The blue zones manifestly correspond to mountain ranges (Pyrenneans and Alps).

The red zones correspond to oceanic basins.

Negative Bouguer anomalies appear in blue where there is a ‘mass deficit’. This occurs in the mountain ranges (Pyrenneans and Alps) because the continental crust is composed of granitic rocks which are lighter than the rocks found deeper in the mantle. Mountains are like icebergs, and are ‘compensated’ by deep roots of light continental crust material. This creates a large mass deficit which shows in blue in the figure. However, note the yellow/red anomaly zone that cuts the Alps (inside the blue zone, down right). This is the Ivrea zone, which is denser and heavier, either because it is the result of rock transformation to another phase (eclogites) or because it is a piece of ioriginally deeper mantle caught into the alpine collision.

In the center of France, there is also a blue negative Bouguer anomaly. This is believed to be due to upwelling of asthenospheric material (lighter in density) or hot spot. Indeed the ‘Massif central’ is a location for volcanism (Chaine des Puys).

Anomalie négative au niveau du massif central : effet du point chaud/ remontée de matériau asthénosphérique moins dense. Sutures hercyniennes détectables.

Let us briefly outline how to calculate the order of magnitude of the Bouguer anomaly.

\[\Delta g = 2 \pi G \int_0^H \rho(y)dy\]

\[\Delta g = 2 \pi \rho_c G H \]

Let us suppose that we have a seamount

#LINKS

Supraconducing gravimeter Ballistic gravimeter

bootstrap table styles to pandoc tables function bootstrapStylePandocTables() { $('tr.header').parent('thead').parent('table').addClass('table table-condensed'); } $(document).ready(function () { bootstrapStylePandocTables(); });