The circumference of an ellipse

Share on:

Note: This blog post is mainly computational, with a hint of proof-oriented mathematics here and there. For a more in-depth analysis, read the excellent article "Gauss, Landen, Ramanujan, the Arithmetic-Geometric Mean, Ellipses, pi, and the Ladies Diary" by Gert Akmkvist and Bruce Berndt, in The American Mathematical Monthly, vol 95 no. 7 (August-September 1988), pages 585-608, and happily made available online for free at https://www.maa.org/sites/default/files/pdf/upload%5Flibrary/22/Ford/Almkvist-Berndt585-608.pdf

This article was the 1989 winner of the MAA Halmos-Ford award "For outstanding expository papers in The American Mathematical Monthly" and you should certainly read it.

Introduction

It's one of the delights or annoyances of mathematics, which ever way you look at it, that there's no simple formula for the cicrcumference of an ellipse comparable to \(2\pi r\) for a circle.

Indeed, for an ellipse with semi axes \(a\) and \(b\), the circumference can be expressed as the integral

\[ 4\int_0^{\pi/2}\sqrt{a^2\cos^2\theta + b^2\sin^2\theta},d\theta \]

which is one of a class of integrals called elliptic integrals, and which cannot be expressed using algebraic or standard transcendental functions.

However, it turns out that there are ways of very quickly computing the circumference on an ellipse to any desired accuracy, using methods which originate with Carl Friedrich Gauss.

The arithmetic-geometric mean

Given \(a>b>0\), we can define two sequences by:

\[ a_0 = a,\qquad a_{k+1}=(a_k+b_k)/2 \]

and

\[ b_0=b,\qquad b_{k+1}=\sqrt{a_kb_k}. \]

Thus the \(a\) values are the arithmetic means of the previous pair; the \(b\) values the geometric mean of the pair.

Since

\[ b<\sqrt{ab}<\frac{a+b}{2}<a \]

and since \(a_{k+1}<a_k\) and \(b_{k+1}>b_k\), it follows that the \(a\) values are decreasing and bounded below, and the \(b\) values are increasing and bounded above, so they both converge. Also, if \(c_k=\sqrt{a_k^2-b_k^2}\) then (see Almqvist \& Berndt pp 587-588):

\begin{align*} c_{k+1}&=\sqrt{a_{k+1}^2-b_{k+1}^2}\\\
&=\sqrt{\frac{1}{4}(a_k+b_k)^2-a_kb_k}\\\
&=\frac{1}{2}(a_k-b_k)\\\
&=\frac{c_k^2}{4a_{k+1}}\\\
&<\frac{c_k^2}{4M(a,b)}. \end{align*}

This shows that not only do the sequences converge to the same limit, but that the sequences converge quadratically; each iteration being double the precision of the previous.

The common limit is called the arithmetic-geometric mean of \(a\) and \(b\) and will be denoted here as \(M(a,b)\).

To give an indication of this speed, in Python:

  from mpmath import mp
  mp.dps = 50
  a,b = mp.mpf('3'), mp.mpf('2')
  for i in range(10):
      a,b = (a+b)/2, mp.sqrt(a*b)
      print('{0:52s} {1:52s}'.format(str(a),str(b)))

2.5                                                  2.4494897427831780981972840747058913919659474806567
2.4747448713915890490986420373529456959829737403283  2.4746160019198827700554874647235766528956885806854
2.4746804366557359095770647510382611744393311605069  2.474680435816873015671798992747485357556612143505
2.4746804362363044626244318718928732659979716520059  2.4746804362363044625888873352899297125054403531594
2.4746804362363044626066596035914014892517060025827  2.4746804362363044626066596035914014892516421855508
2.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
2.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
2.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
2.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
2.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667

You see that we have reached the limit of precision in six steps. To better demonstrate the speed of convergence, use greater precision and display the difference \(a_k-b_k\):

  mp.dps = 2000
  a,b = mp.mpf('3'), mp.mpf('2')
  for i in range(10):
      a,b = (a+b)/2, mp.sqrt(a*b)
      mp.nprint(a-b,10)

0.05051025722
0.0001288694717
8.388628939e-10
3.55445366e-20
6.381703188e-41
2.057141145e-82
2.137563719e-165
2.307963982e-331
2.690598786e-663
3.656695286e-1327

You can see that the precision does indeed roughly double at each step.

Elliptic integrals and the AGM

This integral:

\[ K(k) = \int_0^{\pi/2}\frac{1}{\sqrt{1-k^2\sin^2\theta}}d\theta \]

is called a complete elliptic integral of the first kind. If

\[ k^2=1-\frac{b^2}{a^2} \]

then

\[ \int_0^{\pi/2}\frac{1}{\sqrt{1-k^2\sin^2\theta}} d\theta = a\int_0^{\pi/2} \frac{1}{\sqrt{a^2\cos^2\theta +b^2\sin^2\theta}} d\theta \]

and the integral on the right is denoted \(I(a,b)\). Gauss observed (and proved) that

\[ I(a,b)=\frac{\pi}{2M(a,b)}. \]

This is in fact equivalent to the assertion that

\[ K(x) = \frac{\pi}{2M(1-x,1+x)}. \]

Gauss started with assuming that \(M(1-x,1+x)\) was an even function, so could be expressed as a power series in even powers of \(x\). He then compared the power series representing \(K(x)\) with the integral computed by expanding the integrand as a power series using the binomial theorem and integrating term by term, to obtain

\[ \frac{1}{M(1-x,1+x)}=\frac{2}{\pi}K(x)=1+\left(\frac{1}{2}\right)^{2}x^2+ \left(\frac{1\cdot 3}{2\cdot 4}\right)^{2}x^4+ \left(\frac{1\cdot 3\cdot 5}{2\cdot 4\cdot 6}\right)^{2}x^6+\cdots \]

and this power series can be written as

\[ \sum_{n=0}^{\infty}\frac{(2n)!^2}{2^{4n}(n!)^4}x^{2n}. \]

Another proof, apparently originally due to Newman, uses first the substition \(t=a\tan\theta\) to rewrite the integral:

\[ I(a,b)=\int_0^{\pi/2}\frac{1}{\sqrt{a^2\cos^2\theta+b^2\sin^2\theta}} d\theta = \int^{\infty}_{-\infty}\frac{1}{\sqrt{(a^2+t^2)(b^2+t^2)}} dt \]

and make the substitution

\[ u = \frac{1}{2}\left(t-\frac{ab}{t}\right) \]

which produces (after "some care", according to the Borwein brothers in their book Pi and the AGM):

\[ I(a,b) = \int^{\infty}_{-\infty}\frac{1}{\sqrt{(\left(\frac{a+b}{2}\right)^2+t^2)(ab+t^2)}} dt =I(\frac{a+b}{2},\sqrt{ab}) = I(a_1,b_1). \]

Continuing this process we see that \(I(a,b) = I(a_k,b_k)\) for any \(k\), and taking the limit, that \(I(a,b)=I(M,M)\). Finally

\[ I(M,M) = \int_0^{\pi/2}\frac{1}{\sqrt{M^2\cos^2\theta+M^2\sin^2\theta}} d\theta = \int_0^{\pi/2}\frac{1}{M} d\theta = \frac{\pi}{2M}. \]

So we now know that the complete elliptic integral of the first kind is computable by means of the AGM. But what about the complete elliptic integral of the second kind?

Complete elliptic integrals of the second kind

These are defined as

\begin{align*} E(k)&=\int_0^{\pi/2}\sqrt{1-k^2\sin^2\theta} d\theta\\\
&=\int_0^1\sqrt{\frac{1-k^2t^2}{1-t^2}},dt. \end{align*}

Note: An alternative convention is to write:

\begin{align*} E(m)&=\int_0^{\pi/2}\sqrt{1-m\sin^2\theta} d\theta\\\
&=\int_0^1\sqrt{\frac{1-mt^2}{1-t^2}},dt. \end{align*}

This elliptic integral is the one we want, since the perimeter of an ellipse given in cartesion form as

\[ \frac{x^2}{a^2}+\frac{y^2}{b^2}=1 \]

is equal to

\[ 4aE\left(\sqrt{1-\frac{b^2}{a^2}}\right) \]

(using the first definition). This can be written as

\[ 4aE(e) \]

where \(e\) is the eccentricity of the ellipse.

The computation of elliptic integrals goes back as far as Euler, taking in Gauss and then Legendre along the way. Adrien-Marie Legendre (1752-1833) is better known now than he was in his lifetime; his contributions to elliptic integrals are now seen as fundamental, and in paving the way for the greater work of Abel and Jacobi.

In particular, Legendre showed that

\[ \frac{K(k)-E(k)}{K(k)}=\frac{1}{2}(c_0^2+2c_1^2+4c_2^2+\cdots + 2^nc_n^2+\cdots ) \]

where

\[ c_n^2=a_n^2-b_n^2. \]

The above equation can be rewritten to provide a fast-converging series for \(E(k)\):

\begin{align*} E(k)&=K(k)(1-\frac{1}{2}(c_0^2+2c_1^2+4c_2^2+\cdots + 2^nc_n^2+\cdots))\\\
&=\frac{\pi}{2M(a,b)}(1-\frac{1}{2}(c_0^2+2c_1^2+4c_2^2+\cdots + 2^nc_n^2+\cdots)). \end{align*}

In the sequences for the AGM, let \(M(a,b)\) be approximated by \(a_n\). This produces a sequence that converges very quickly to \(E(k)\):

\begin{align*} e_0&=\frac{\pi}{2a}(1-\frac{1}{2}c_0^2)\\\
e_1&=\frac{\pi}{2a_1}(1-\frac{1}{2}(c_0^2+2c_1^2))\\\
e_2&=\frac{\pi}{2a_2}(1-\frac{1}{2}(c_0^2+2c_1^2+4c_2^2)) \end{align*}

This can easily be managed recursively, given \(a\) and \(b\), by setting

\[ a_0=a,\qquad,b_0=b,\qquad p_0=1,\qquad s_0=a^2-b^2 \]

and then iterating by

\begin{align*} a_{k+1}&=(a_k+b_k)/2\\\
b_{k+1}&=\sqrt{a_kb_k}\\\
p_{k+1}&=2p_k\\\
s_{k+1}&=s_k+p_k(a_{k+1}^2-b_{k+1}^2) \end{align*}

Then the values

\[ e_k = \frac{2\pi}{a_k}(a_0^2-\frac{1}{2}s_k) \]

approach the perimeter of the ellipse.

We can demonstrate this in Python for \(a,b=3,2\), again using the multiprecision library mpmath:

  mp.dps = 100
  a = mp.mpf(3)
  b = mp.mpf(2)
  p = 1
  s = a**2-b**2
  a0 = a
  e = 2*mp.pi/a*(a0*a0-1/2*s)
  for i in range(11):
      a, b, p = (a+b)/2, np.sqrt(a*b), p*2
      s += p*(a**2-b**2)
      print(2*mp.pi/a*(a0*a0-1/2*s))

15.70796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534
15.86502654157467714117218441586859215641446367453641805450510958042565274839788072998452444741042018
15.86543958660157016021207778974106592713197464956329054113604818310136964968401110893643543265751673
15.86543958929058979121772312468572564968472056328417484804713380278017532700185422580941088121856672
15.86543958929058979133166302778307249672987827943566144626073834574026314097232621543195947183210587
15.86543958929058979133166302778307249673008284832650068966726311774248223910968899525488214058871233
15.86543958929058979133166302778307249673008284832650068966726311774248223910968899591430967903912196
15.865439589290589791331663027783072496730082848326500689667263117742482239109688995914309679039122
15.86543958929058979133166302778307249673008284832650068966726311774248223910968899591430967903912207
15.86543958929058979133166302778307249673008284832650068966726311774248223910968899591430967903912222
15.86543958929058979133166302778307249673008284832650068966726311774248223910968899591430967903912252

We see that we have reached the limits of this precision very quickly. And as above, we can look instead at the differences between succesive approximations:

mp.dps=2000
a = mp.mpf(3)
b = mp.mpf(2)
p = 1
s = a**2-b**2
a0 = a
e = 2*mp.pi/a*(a0*a0-1/2*s)
for i in range(11):
    a, b, p = (a+b)/2, np.sqrt(a*b), p*2
    s += p*(a**2-b**2)
    e1 = 2*mp.pi/a*(a0*a0-1/2*s)
    mp.nprint(e1-e,10)
    e = e1

2.094395102
0.1570632736
0.0004130450269
2.689019631e-9
1.139399031e-19
2.045688908e-40
6.594275385e-82
6.852074221e-165
7.398301331e-331
8.624857552e-663
1.172173128e-1326

The convergence is seen to be quadratic.

Using Excel

It might seem utterly regressive to use Excel - or indeed, any spreadsheet - for computations such as this, but in fact Excel can be used very easily for recursive computations, as long as you're happy to use only IEEE double precision.

In Cells A1-E2, enter:

A B C D E
1 3 2 1 =A1^2-B1^2
2 =(A1+B1)/2 =SQRT(A1*B1) =2*C1 =D1+C2*(A2^2-B2^2) =2*PI()/A2*($A$1^2-1/2*D2)

and then copy A2-E2 down a few rows. Because of overflow errors, if you show 14 decimal places in column E, they'll never quite settle down. But in fact you reach the limits of double precision by row 5.