The circumference of an ellipse

Apr 10, 2021 mathematics computation

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_library/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:

 1  from mpmath import mp
 2  mp.dps = 50
 3  a,b = mp.mpf('3'), mp.mpf('2')
 4  for i in range(10):
 5      a,b = (a+b)/2, mp.sqrt(a*b)
 6      print('{0:52s} {1:52s}'.format(str(a),str(b)))
 7
 82.5                                                  2.4494897427831780981972840747058913919659474806567
 92.4747448713915890490986420373529456959829737403283  2.4746160019198827700554874647235766528956885806854
102.4746804366557359095770647510382611744393311605069  2.474680435816873015671798992747485357556612143505
112.4746804362363044626244318718928732659979716520059  2.4746804362363044625888873352899297125054403531594
122.4746804362363044626066596035914014892517060025827  2.4746804362363044626066596035914014892516421855508
132.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
142.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
152.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
162.4746804362363044626066596035914014892516740940667  2.4746804362363044626066596035914014892516740940667
172.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$:

 1  mp.dps = 2000
 2  a,b = mp.mpf('3'), mp.mpf('2')
 3  for i in range(10):
 4      a,b = (a+b)/2, mp.sqrt(a*b)
 5      mp.nprint(a-b,10)
 6
 70.05051025722
 80.0001288694717
 98.388628939e-10
103.55445366e-20
116.381703188e-41
122.057141145e-82
132.137563719e-165
142.307963982e-331
152.690598786e-663
163.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:

 1  mp.dps = 100
 2  a = mp.mpf(3)
 3  b = mp.mpf(2)
 4  p = 1
 5  s = a**2-b**2
 6  a0 = a
 7  e = 2*mp.pi/a*(a0*a0-1/2*s)
 8  for i in range(11):
 9      a, b, p = (a+b)/2, np.sqrt(a*b), p*2
10      s += p*(a**2-b**2)
11      print(2*mp.pi/a*(a0*a0-1/2*s))
12
1315.70796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058534
1415.86502654157467714117218441586859215641446367453641805450510958042565274839788072998452444741042018
1515.86543958660157016021207778974106592713197464956329054113604818310136964968401110893643543265751673
1615.86543958929058979121772312468572564968472056328417484804713380278017532700185422580941088121856672
1715.86543958929058979133166302778307249672987827943566144626073834574026314097232621543195947183210587
1815.86543958929058979133166302778307249673008284832650068966726311774248223910968899525488214058871233
1915.86543958929058979133166302778307249673008284832650068966726311774248223910968899591430967903912196
2015.865439589290589791331663027783072496730082848326500689667263117742482239109688995914309679039122
2115.86543958929058979133166302778307249673008284832650068966726311774248223910968899591430967903912207
2215.86543958929058979133166302778307249673008284832650068966726311774248223910968899591430967903912222
2315.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:

 1mp.dps=2000
 2a = mp.mpf(3)
 3b = mp.mpf(2)
 4p = 1
 5s = a**2-b**2
 6a0 = a
 7e = 2*mp.pi/a*(a0*a0-1/2*s)
 8for i in range(11):
 9    a, b, p = (a+b)/2, np.sqrt(a*b), p*2
10    s += p*(a**2-b**2)
11    e1 = 2*mp.pi/a*(a0*a0-1/2*s)
12    mp.nprint(e1-e,10)
13    e = e1
14
152.094395102
160.1570632736
170.0004130450269
182.689019631e-9
191.139399031e-19
202.045688908e-40
216.594275385e-82
226.852074221e-165
237.398301331e-331
248.624857552e-663
251.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)`	`=2PI()/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.