*   Contents

• 1 Simple Pendulum : Period Etc.
• 2 Evaluating the Arithmetic-Geometric Mean
• 3 Connecting the Physics to the agM
• 4 Approximations for the Asymptotic Cases
  • 4.1 Small Amplitudes
  • 4.2 Large Amplitudes
• 5 Deriving The Equation of Motion
  • 5.1 Basic Physics
  • 5.2 Putting the Integral into Standard Form
  • 5.3 Gory Details: Standard Forms of the Elliptic Integral
• 6 The Sequence of Integrals
  • 6.1 Overview
  • 6.2 Gory Details: Cosine Factor
  • 6.3 Gory Details: Radius Factor
  • 6.4 Gory Details: Numerator
• 7 Ramifications
  • 7.1 Numerical Robustness
  • 7.2 Backwards Sequencing
  • 7.3 Eponymous Constants
• 8 Assumptions and Details
• 9 Software Implementation : C++ Code
• 10 Notes and References

1 Simple Pendulum : Period Etc.

Consider the simple pendulum shown in figure 1. The center-of-mass moves along a circle of radius L. The instantaneous position is given by the angle φ.

Figure 1: Simple Pendulum

We are interested in the pendulum’s period of oscillation, denoted T. Pendulums have been used for timekeeping since Day One of modern science, and have been used to symbolize physics itself.

Figure 2: Old AIP Logo

The period depends on the length of the pendulum (L) and on the local acceleration of gravity (g), both of which we assume are well-known and constant.¹

We need to look closely at the dependence of the period on the amplitude of the oscillation, θ. On the one hand, the period is “nearly” independent of amplitude when the amplitude is small. On the other hand, a very good pendulum clock keeps time to about 0.1 ppm, so even a weak dependence is significant to horologists. Even in the introductory physics lab, it is easy to get into the regime where this dependence affects the measurements.

Here is a very useful formula for the period T. It is exact, easy to remember, and easy to use:

As derived in section 5,

where agM() is the arithmetic-geometric mean, as discussed in detail below. It splits the difference between the familiar arithmetic mean and geometric mean. It can be calculated very efficiently. Forsooth, it is easier to calculate agM() from scratch that it is to calculate cos() from scratch.

Figure 3 shows how the period depends on amplitude.

Figure 3: Pendulum Period and Inverse Period versus Amplitude

Table 1 provides similar information in tabular form, namely the depencence of period on amplitude:

(See equation 47 for a version of equation 1 with better numerical stability in extreme cases.)

(Beware: In the literature there are dozens of formulas for approximating the period, using power series or whatever, but these are almost never useful, because the exact expression in equation 1 is simpler and better.)

Important remark: There area hundreds, maybe thousands of references that say there is no closed-form solution for the period of a pendulum, except in the small-angle approximation. They analyze small-angle case and then don’t even try to analyze larger angles. They just give up.

The point is, I consider equation 1 to be a perfectly respectable closed-form solution. The agM is perfectly well defined and easy to evaluate. Forsooth, it is easier to evaluate the agM from scratch than evaluate the cosine from scratch.

The equation of motion leads to a complete elliptic integral of the first kind. Such things crop up fairly often in math and physics. They have an undeserved reputation for being hard to deal with. If you encounter one, don’t give up. Evaluate it in terms of the agM.

There’s nothing particularly new about this. Gauss figured out how to evaluate elliptic integrals via the agM about 200 years ago. The work wasn’t published until some decades later, but still, there’s been plenty of time for the idea to spread. It’s mentioned in Abramowitz and Stegun (1964), including in the table of contents.² Back before everything was available online, generations of physicists relied on that book. Also, the Wikipedia page for “elliptic integral” has mentioned the agM since the second edit, more than 20 years ago (March, 2002).³ There’s also a nice AJP article from 2008, discussing the agM and its application to pendulums.⁴ So it’s not a big secret.

2 Evaluating the Arithmetic-Geometric Mean

Suppose we are given two numbers, a₀ and b₀, and we wish to find their arithmetic-geometric mean, denoted agM(a₀, b₀). Then we define a sequence of values, for all i≥0:

This is interesting because the sequence converges very quickly. After a few interations, a_i is equal to b_i for all practical purposes, and their shared value defines the agM.

Equivalently, we can define the agM recursively:

We remark that the agM has the following properties, which you would expect any sort of average to have:

Furthermore, the agM sequence has the wonderful property that for all i (except maybe i=0),⁵ a_i is an upper bound on the final answer, and b_i is a lower bound.

This means you always know how well converged the sequence is. This is unlike most sequences and series, where you have to do extra work to figure out how well converged things are.

Some c++ code to evaluate the agM can be found in section 9.

3 Connecting the Physics to the agM

So far we have merely asserted equation 1. Let’s see if we can derive it.

Starting from the basic equations of motion, we can derive a formula that expresses the period of oscillation (T) in terms of a complete elliptic integral of the first kind. The derivation is spelled out in section 5.

where for a pendulum a+0=1 and b₀=cos(θ/2). By symmetry, we can also write this as:

which doesn’t tell us much we didn’t already know, but makes the geometric interpretation more clear. The RHS is the average of 1/r for an ellipse centered at the origin, averaged over all angles, where the semi-major axis is a₀ and the semi-minor axis is b₀.

To make progress, we need a way of evaluating the integral. You could do it numerically (using Runge-Kutta or whatever), but there is a cleverer way.

Key idea: Suppose we have a sequence of ellipses, all with the same average 1/r, and suppose the sequence rapidly converges to an ellipse that is circular. Such a sequence is shown in figure 4.

Figure 4: The agM as a Sequence of Ellipses

The average 1/r for a circle is obvious. Game over, you win.

The only remaining task is to find a suitable sequence of ellipses. It turns out that the a_i and b_i defined by equation 2 do the job very nicely.

Specifically: In equation 7, we make the following substitution:^6,7

when we turn the crank,^7,8 we find that

where the a_i and b_i have the same meaning as in equation 2. By iterating this process, we find that T/T₀ converges to 1/agM(a, b), as advertised (i.e. equation 1).

In figure 4 the red ellipse is of particular relevance. It represents agM(1, √½), which is the starting point for evaluating T/T₀ when the amplitude is θ=90^∘. The dashed blue circle is an accurate representation of the final answer, T₀/T=0.8472, obtained by two iterations of the agM algorithm.

4 Approximations for the Asymptotic Cases

The literature is overflowing with complicated approximations to the elliptic integral. Most of these are worthless, because the exact expression equation 1 is easier to understand, and easier to use.

However, sometimes you can improve your intuition by looking at a low-order expansion in the neighborhood of the limiting cases: small amplitude (theta near zero) and large amplitudes (η near zero).

4.1 Small Amplitudes

The small-amplitude approximation is:

Here’s a much more accurate approximation, sometimes useful if you’re using a spreadsheet that can do trig functions and square roots, but can’t execute imperative iterative code (as in section 9). It is equivalent to doing one 1.5 rounds by hand, algebraically (as in equation 14):

Keep in mind that even with a spreadsheet, even if you can’t do arbitrarily many iterations, you can still do enough interations to obtain very high accuracy, directly from the definition of agM. If necessary, add a couple of columns to the spreadsheet, to hold the intermediate a_i and b_i.

Let’s conpare this to the exact value when θ is 45^∘:

So equation 10 is in the right ballpark. For smaller amplitudes, its accuracy is better.

For larger amplitudes, or if you need more accuracy, don’t bother using a power series with more terms. Instead, use more iterations of the agM algorithm. Equation 11 is accurate within 20 ppm all the way out to θ=90^∘.

Details: To derive equation 10, start with equation 1 and expand everything to lowest order:

To derive equation 11, the key idea is this: Directly using the definition of agM, peform the iterations using algebraic variables (rather than numbers). That gives us:

It’s rarely necessary to continue the process beyond this point, but you can if you want. Explicit expressions are given in reference 4.

4.2 Large Amplitudes

We now consider the case of large amplitudes, where θ is large and the supplementary angle η approaches zero.

In this case, the period approaches infinity. However, the approach is remarkably slow, as shown in figure 5. In fact, for η near zero, the period is proportional to the logarithm of η.

Figure 5: The Period Diverges Very Slowly

In the figure, note that each tick mark in the horizontal direction represents two orders of magnitude of change in η. Here’s a supplement to table 1, quantifying the large-angle behavior.

Here is an interesting expression for agM(1,b), valid when b is small:⁹

This is useful insofar as it illuminates the theory. It shows that the logarithmic divergence in figure 5 is not a fluke. (On the other hand, it is not of much use in numerical calculations, since the exact expression in equation 1 is simpler and better.)

In particular, here’s an amusing exercise that requires more theorizing than computing: Ideally, if you balance a pencil so that it is vertical, it should stay there for a very long time before falling over. Using equation 15 (or otherwise), estimate how hard it would be to balance a pencil so precisely that it takes at least 20 seconds to fall over. (Assume the tip is sharpened so that its radius of curvature is negligible.)

5 Deriving The Equation of Motion

Let’s derive equation 1.

5.1 Basic Physics

We start by considering the total energy, which we take to be a constant,¹ denoted E.

The kinetic part of the energy is:

At this point we invoke the oscillator assumption.¹ That is, we assume that it’s an oscillator, not a rotor. In other words, it doesn’t have enough energy to flip over the top and keep going in the same direction. The turnaround occurs at an angle θ, as shown in the diagram, at which point the kinetic energy is zero. At other points in the swing, the potential energy differs in proportion to the difference in height, so the kinetic energy must do the same.

Note that we did not directly use Newton’s laws of motion. The force law would have given us an equation involving the second derivative. Here we used energy concepts (not known in Newton’s day) to obtain an equation involving the square of the first derivative, which saves us a couple of steps and makes the physical significance of equation 17 more clear.

Let’s integrate both sides of equation 17. Integrating from bottom dead center to the turning point gives us a time equal to one quarter of the oscillation period (T):

5.2 Putting the Integral into Standard Form

We are about to show that equation 18 can be written as:

where the integral is the textbook standard form of a complete elliptic integral of the first kind with modulus k. For our pendulum, k is given by:

We can also rewrite the period as:

or equivalentely

where the integral in equation 21 (not including the 2/π out front) is the standard two-argument Cayley form of the complete elliptic integral of the first kind. For the pendulum we have:

Equation 22 has a nifty geometric interpretation: The RHS is just the average of 1/r, averaged over all angles, as we go around an ellipse with semi-major axis a and semi-minor axis b. We have already used this idea in section 3.

5.3 Gory Details: Standard Forms of the Elliptic Integral

This subsection is mostly routine algebra (although there is some trickery in equation 26). Feel free to skim, or skip ahead to section 6.

We can rewrite the trig functions in terms of the squares of of trig functions, using the double-angle identities, which are valid for any angle α:

This makes the integral look a little more “elliptical”:

It is slightly annoying to have θ appear both in the integrand and in the limit of integration. We can fix that, and bring the integral into a more conventional form, using the following substitution:

The idea is that at the turning point, where φ=θ, we have sin(u)=1 and u=π/2, which is a constant, independent of θ. The substitution is also motivated by the appearance of the square of sin(u) inside the square root. Most importantly, the substitution is motivated by 20/20 hindsight. Clever people have been down this road before, so we know where it leads. To make use of the substitution, we need its derivative:

So the whole integral becomes:

where the integral in the last line is a standard representation of the complete elliptic integral of the first kind.^10,8 In our case the modulus is k=sin(θ/2) or equivalently the parmeter is m=k²=sin²(θ/2).

When the amplitude goes to zero, the integral is trivial to evaluate. This gives us:

which allows us to write an important result. The RHS of equation 30 is the standard form of the complete elliptic integral of the first kind.

To obtain the two-argument Cayley form, we replace the 1 in equation 30 with sin² + cos².

These results are discussed and applied back in section 5.2

6 The Sequence of Integrals

6.1 Overview

We wish to prove equation 9, which we restate here:

using the notation for a_i and b_i as defined in equation 2.

It’s not obvious how to proceed, so we rely on Gauss to help us get started. We make the following substitution:

In some sense, that’s all you need to know. The rest is just algebra. Pages and pages of algebra.

Remark: When u is zero, u′ is also zero. Also, when u is π/2, u′ is also π/2. It pretty much has to be this way, on account of the symmetry of the ellipses in figure 4.

We restate equation 31 with a trivial modification, and then organize the rest of the work in three stages.

Stage 1: In section 6.2, we will show that the cosine factor is:

where D is the denominator of equation 33, namely:

In equation 35 we are happy to see a₁ and b₁. The square root can be visualized as the radius of an ellipse with semi-major axis a₁ and semi-minor axis b₁, which is concept we have seen a couple of times already.

Stage 2: In section 6.3, we will show that the radius factor is:

Stage 3: In section 6.4, we will show the numerator of equation 34 is:

We plug these subsidiary results into equation 34, and almost everything cancels out:

just as advertised.

The following three subsections are just algebra. They’re not very interesting. They are included mainly to reassure you that you could do the calculations if you wanted. No wizardly required. The only serious wizardry was back in equation 33, along with the underlying idea of a sequence of integrals.

Feel free to skip ahead to section 7

6.2 Gory Details: Cosine Factor

To obtain a formula for cos(u), we start with equation 33, then multiply both sides by D and turn the crank:

We would like to pull out a factor of cos²(u′) on the right, to go with the cos²(u) on the left. And once again we are motivated by 20/20 hindsight. The calculation continues:

6.3 Gory Details: Radius Factor

Next, let’s work out the other factor in the denominator of equation 34, namely the square root, i.e. the “old” radius. In that factor, we use equation 33 to plug in for for sin(u), and use equation 42 to plug in for cos(u). Then turn the crank:

The RHS is a perfect square, so we are left with an expression that is a lot simpler than the steps that led to it:

6.4 Gory Details: Numerator

Last but not least, we need to deal with the numerator in equation 34. In particular, we need to express it in terms of u′. We do that by restating equation 33 and differentiating both sides:

So:

7 Ramifications

Here’s some stuff that non-experts don’t need to know, but may find amusing.

7.1 Numerical Robustness

Equation 1 is mathematically correct and is directly usable in all but the most extreme situations. However, when doing floating-point calculations, roundoff becomes a problem when the amplitude θ is very large, specifically when the supplementary angle η is very small compared to 1. In such situations you should use the last line of equation 47, which is mathematically equivalent but numerically more robust. It remains usable down to η=10⁻¹⁰⁰ and beyond.

This was useful for preparing figure 5 and table 2.

7.2 Backwards Sequencing

In the sequence that defines the agM, you can move backwards using the expressions:¹¹

Note: In such a series, a_i ≥ b_i for all i (except possibly i=0),⁵ as a direct consequence of the definition.

7.3 Eponymous Constants

Gauss’s Constant: The reciprocal of agM(1,√2) is known as Gauss’s constant, G≈0.834627. It shows up in various contexts in analysis.¹² The period for a pendulum when the amplitude is 90^∘ is T/T₀=√2 G.

Galileo’s Constant: The quarter-period is the time it takes a pendulum to fall to the middle, starting from rest at some angle θ. This is proportional to the time it takes for an object to a distance L, starting from rest, where L is the length of the pendulum. For small amplitudes, this is a universal constant, namely π / √8. Stillman Drake calls this Galileo’s constant,¹³ although I’m not sure anybody else does. Galileo had to determine this quantity experimentally, since the necessary theory (Newtonian mechanics) hadn’t been invented yet.

8 Assumptions and Details

Our calculations are based on the following assumptions.

1.

Oscillator: We assume that the pendulum is actually pendulous, i.e. it hangs down.

In other words, we assume it is an oscillator, not a rotor. It doesn’t have enough energy to flip over the top and keep going in the same direction.

2.

Stable conditions: We assume the pivot is well supported and not wobbling. We assume the length L is well-known and not changing. Similarly we assume the local acceleration of gravity g is well-known and not changing.

Beware that g varies from place to place, by an amount that is easily measurable with a good pendulum clock.

3.

Moving in a single vertical plane.

It is certainly possible to have a spherical pendulum. When the amplitude is small, this is not too bad. In the special case of a conical pendulum, where the bob moves in a horizonal circle, this is not to bad. However, the general case is very complicated.

4.

Lossless: We neglect friction of all kinds. In other words, the system energy E is constant.

5.

Airless: We neglect buoyancy, breezes, and all other static and dynamic contributions from the air.

6.

Simple pendulum.

That is, not a compound pendulum, i.e. a pendulum attached to another pendulum.

7.

Rigid: In cases where the amplitude θ is larger than 90∘, we assume the bob is connected to the pivot by a rigid rod.

You can use a string for small angles, but not large angles, because it would go slack.

8.

Bob Only: We assume the pendulum has the conventional design, so that the mass of the connecting rod is negligible compared to the mass of the bob.

For an arbitrary distribution of mass, deriving the equation of motion would be slightly more complicated. The potential energy would still be proportional to −cos(φ), and kinetic energy would still be proportional to (dφ/dt)2, so the calculation runs parallel the one in section 5; you just have to figure out the proportionality factors.

9.

Non-Thermal: We ignore thermal fluctuations.

In the classical approximation, if the mass is large enough and you cool things down enough, the thermal fluctuations are negligible.

10.

Classical: In rough terms, we assume the system is “classical” in the sense that we don’t need to worry about relativistic effects.

In more precise terms, everything is subject to the laws of relativity. Momentum is first order in v/c, and kinetic energy is second order in v/c, so what we are really saying is the third-order (and higher) terms are neglgible. To put it another way, phenomena that were well known prior to 1887, prior to Michelson-Morley, are considered classical.

Similarly, we neglect general relativity. This is an excellent approximation.

Again, everything is subject to the laws of general relativity. In practice, we don’t need to worry about it when the distances are not enormously large, and the gravitational fields are not enormously strong.

We also neglect quantum effects. This is not an entirely trivial assumption. QM places some rather strict bounds on how long you can make a pencil balance on end.

Again, everything is subject to the laws of QM. Quantum mechanics contains, predicts, and explains classical mechanics — but not vice versa. In practice, if the bob is massive enough and you don’t look too closely, the equations of classical physics are good enough.

11.

In the definition of agM(a, b) we assume that a and b are positive real numbers. For the pendulum we are considering, this is automatically true, as a consequence of the other assumptions.

The agM function can by analytically continued to arbitrary complex numbers. This has some interesting properties, but is of no help in analyzing the pendulum.

9 Software Implementation : C++ Code

// <*><*><*><*>  arithmetic-geometric-mean.c  <*><*><*><*>
#include <cmath>
#include <stdexcept>
#include "arithmetic-geometric-mean.h"
using namespace std;
static double const qnan = numeric_limits<double>::quiet_NaN();

// Simple version of the arithmetic-geometric mean.
// A first step, for edification.
// Not industrial strength.
double agM_simple(double ai, double bi) {
  int maxit(20);
  for (int iter=0;;) {
    if (ai == bi) break;
    if (iter >= maxit) throw runtime_error("agM didn't converge");
    double anew = (ai+bi)/2.;
    double bnew = sqrt(ai*bi);
    ++iter;               // count # of iterations needed
    if (fabs(anew-bnew) >= fabs(ai-bi)) {
// previous values (ai,bi) were just as good, maybe better
      break;
    }
    ai = anew;
    bi = bnew;
  }
  return ai;
}

// Fancy version.
// More defensive about extreme cases and out-of-bounds cases.
//
// Some typical usages, among others:
// Simplest case:
//   mean = agM(a,b);
//   mean = agM(a);     // same as agM(a,1.)
// More information is available from the structured object:
//   AGM foo(a,b); mean = foo.ai;  iter = foo.iter;
//   AGM foo; mean = foo.agM(a,b); iter = foo.iter;

double AGM::agM(double const a0, double const b0) {
  ai = a0;
  bi = b0;
  iter = 0;
  if (bi == 0) return (ai=bi);
  if (ai == 0) return (bi=ai);
  if (ai-ai != 0) return (bi=ai);      // defend against nan and inf
  if (bi-bi != 0) return (ai=bi);
  if (bi < 0) return ai=bi=qnan;
  if (ai < 0) return ai=bi=qnan;
  for (;;) {
    if (ai == bi) break;
// check for badness /after/ checking for goodness:
    if (iter >= maxit) throw runtime_error("agM didn't converge");
// Two separate square roots, to defend against the possibility
// that the product of ai*bi would cause overflow or underflow.
// Separation slows down the typical (non-nasty) case, and I do
// know how to optimize it, but overall the routine is so efficient
// that optimizing it would be silly.
// FWIW: In non-nasty cases the root of the product may differ from
// the product of the roots by 1 ULP (unit in the last place).
    bnew = sqrt(ai) * sqrt(bi);
    anew = (ai+bi)/2.;
    ++iter;               // count # of iterations needed
    if (fabs(anew-bnew) >= fabs(ai-bi)) {
// previous values (ai,bi) were just as good, maybe better
      break;
    }
    ai = anew;
    bi = bnew;
  }
  return ai;
}

// maxit is 20 because 15 iterations are needed
// when a=1e300 and b=1e-300
//
// except for maxit, these are all out-of-bounds values,
// to make sure that anybody who uses them initializes them.
AGM::AGM()
 : maxit(20), iter(10*maxit), ai(-1.), bi(-1), anew(-1.), bnew(-1.) {}

AGM::AGM(double const _a, double const _b) : AGM() {
  agM(_a, _b);
}

// not a class member:
double agM(double const a0, double const b0) {
  return AGM(a0, b0).ai;
}

// <*><*><*><*>  arithmetic-geometric-mean.h  <*><*><*><*>
#pragma once
double agM_simple(double a0, double b0);
struct AGM {
  int maxit;
  int iter;
  double ai, bi, anew, bnew;
  double agM(double const _a, double const _b=1);
// constructors:
  AGM();        // doesn't do any work
  AGM(double const _a, double const _b=1);
};

// not a class member:
double agM(double const a0, double const b0=1);

10 Notes and References

1.

In section~8 there is a list of the assumptions we are making.

2.

Abramowitz and Stegun, Handbook of Mathematical Functions (1964)
https://personal.math.ubc.ca/~cbm/aands/

3.

Wikipedia article, “Elliptic Integral” https://en.wikipedia.org/wiki/Elliptic_integral

4.

Claudio Carvalhaes and Patrick Suppes, “Approximation for the period of the simple pendulum based on the arithmetic-geometric mean”
Am. J. Phys., Vol. 76, No. 12, December 2008 http://dx.doi.org/10.1119/1.2968864

5.

By symmetry, without loss of generality, we can assume a_0≥b_0, in which case equation~5 and equation~48 apply for all i including i=0.

6.

David A. Cox, “Gauss and the Arithmetic-Geometric Mean”
https://ctnt-summer.math.uconn.edu/wp-content/uploads/sites/1632/2016/02/coxctnt.pdf#page=9

7.

David A. Cox, “The arithmetic-geometric mean of Gauss” https://webspace.science.uu.nl/~wepst101/elliptic/cox_agm.pdf

8.

Djalil Chafaï, “Proof of the invariance of Cayley elliptic integrals”
https://djalil.chafai.net/blog/2021/06/16/landen-transformation-of-complete-elliptic-integrals/

9.

See equation 16 in Eric W. Weisstein, “Arithmetic-Geometric Mean”
https://mathworld.wolfram.com/Arithmetic-GeometricMean.html

10.

B.C. Carlson, “Elliptic Integrals”
Chapter 19 in the NIST Digital Library of Mathematical Functions
https://dlmf.nist.gov/19

11.

P. Van Mieghem, “The Arithmetic-Geometric Mean: A Pearl of Gauss” https://www.nas.ewi.tudelft.nl/people/Piet/papers/TUD20230605_arithmetic_geometric_mean.pdf

12.

John D. Cook, “Gauss’s Constant”
https://www.johndcook.com/blog/2021/10/17/gauss-constant/

13.

Stillman Drake, “Galileo’s Constant”
https://www.degruyter.com/document/doi/10.3138/9781487572037-029/html