The circumference of an ellipse
Note: This blog post is mainly computational, with a hint of prooforiented mathematics here and there. For a more indepth analysis, read the excellent article "Gauss, Landen, Ramanujan, the ArithmeticGeometric Mean, Ellipses, pi, and the Ladies Diary" by Gert Akmkvist and Bruce Berndt, in The American Mathematical Monthly, vol 95 no. 7 (AugustSeptember 1988), pages 585608, and happily made available online for free at https://www.maa.org/sites/default/files/pdf/upload_library/22/Ford/AlmkvistBerndt585608.pdf
This article was the 1989 winner of the MAA HalmosFord 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 arithmeticgeometric 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^2b_k^2}\) then (see Almqvist \& Berndt pp 587588):
\begin{align*} c_{k+1}&=\sqrt{a_{k+1}^2b_{k+1}^2}\\ &=\sqrt{\frac{1}{4}(a_k+b_k)^2a_kb_k}\\ &=\frac{1}{2}(a_kb_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 arithmeticgeometric 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_kb_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(ab,10)
6
70.05051025722
80.0001288694717
98.388628939e10
103.55445366e20
116.381703188e41
122.057141145e82
132.137563719e165
142.307963982e331
152.690598786e663
163.656695286e1327
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{1k^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{1k^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(1x,1+x)}. \]
Gauss started with assuming that \(M(1x,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(1x,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{1k^2\sin^2\theta} d\theta\\ &=\int_0^1\sqrt{\frac{1k^2t^2}{1t^2}}\,dt. \end{align*}
Note: An alternative convention is to write:
\begin{align*} E(m)&=\int_0^{\pi/2}\sqrt{1m\sin^2\theta} d\theta\\ &=\int_0^1\sqrt{\frac{1mt^2}{1t^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. AdrienMarie Legendre (17521833) 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^2b_n^2. \]
The above equation can be rewritten to provide a fastconverging 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^2b^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}^2b_{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**2b**2
6 a0 = a
7 e = 2*mp.pi/a*(a0*a01/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**2b**2)
11 print(2*mp.pi/a*(a0*a01/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**2b**2
6a0 = a
7e = 2*mp.pi/a*(a0*a01/2*s)
8for i in range(11):
9 a, b, p = (a+b)/2, np.sqrt(a*b), p*2
10 s += p*(a**2b**2)
11 e1 = 2*mp.pi/a*(a0*a01/2*s)
12 mp.nprint(e1e,10)
13 e = e1
14
152.094395102
160.1570632736
170.0004130450269
182.689019631e9
191.139399031e19
202.045688908e40
216.594275385e82
226.852074221e165
237.398301331e331
248.624857552e663
251.172173128e1326
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 A1E2, enter:
A  B  C  D  E  

1  3 
2 
1 
=A1^2B1^2 

2  =(A1+B1)/2 
=SQRT(A1*B1) 
=2*C1 
=D1+C2*(A2^2B2^2) 
=2*PI()/A2*($A$1^21/2*D2) 
and then copy A2E2 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.