Area and volume of hypersphere cap

A spherical cap is the portion of a sphere above some horizontal plane. For example, the polar ice cap of the earth is the region above some latitude. I mentioned in this post that the area above a latitude φ is

A = 2\pi R^2(1-\sin\varphi)

where R is the earth’s radius. Latitude is the angle up from the equator. If we use the angle θ down from the pole, we get

A = 2\pi R^2(1-\cos\theta)

I recently ran across a generalization of this formula to higher-dimensional spheres in [1]. This paper uses the polar angle θ rather than latitude φ. Throughout this post we assume 0 ≤ θ ≤ π/2.

The paper also includes a formula for the volume of a hypersphere cap which I will include here.

Definitions

Let S be the surface of a ball in n-dimensional space and let An(R) be its surface area.

A_n(R) = \frac{\pi^{n/2}}{\Gamma(n/2)} R^{n-1}

Let Ix(a, b) be the incomplete beta function with parameters a and b evaluated at x. (This notation is arcane but standard.)

I_x(a, b) = \int_0^x t^{a-1}\, (1-t)^{b-1}\, dt

This is, aside from a normalizing constant, the CDF function of a beta(a, b) random variable. To make it into the CDF, divide by B(a, b), the (complete) beta function.

B(a, b) = \int_0^1 t^{a-1}\, (1-t)^{b-1}\, dt

Area equation

Now we can state the equation for the area of a spherical cap of a hypersphere in n dimensions.

A_n^{\text{cap}}(R) = \frac{1}{2}A_n(R)\, I_{\sin^2\theta}\left(\frac{n-1}{2}, \frac{1}{2} \right )

Recall that we assume the polar angle θ satisfies 0 ≤ θ ≤ π/2.

It’s not obvious that this reduces to the equation at the top of the post when n = 3, but it does.

Volume equation

The equation for the volume of the spherical cap is very similar:

V_n^{\text{cap}}(R) = \frac{1}{2}V_n(R)\, I_{\sin^2\theta}\left(\frac{n+1}{2}, \frac{1}{2} \right )

where Vn(R) is the volume of a ball of radius R in n dimensions.

V_n(R) = \frac{\pi^{n/2}}{\Gamma\left(\frac{n}{2} + 1\right)} R^n

Related posts

[1] Shengqiao Li. Concise Formulas for the Area and Volume of a Hyperspherical Cap. Asian Journal of Mathematics and Statistics 4 (1): 66–70, 2011.

Random points in a high-dimensional orthant

In high dimensions, randomly chosen vectors are very likely nearly orthogonal. I’ll unpack this a little bit then demonstrate it by simulation. Then I’ll look at what happens when we restrict our attention to points with positive coordinates.

***

The lengths of vectors don’t contribute to the angles between them, so we may as well focus our attention on spheres.

Suppose we pick two points on a familiar sphere in 3 dimensions, our home planet. By symmetry, we can assume the first vector runs from the center of the earth to the north pole. Otherwise we can imagine picking the first point, then rotating the earth so that the chosen point becomes the north pole.

What vectors are orthogonal (perpendicular) to the north pole? The vectors from the center to a point on the equator. The vectors that are nearly orthogonal to the north pole are the vectors that are near the equator.

On a high-dimensional sphere, nearly all the area is near the equator (more on that here). This means that if we choose two vectors, they are very likely to be nearly orthogonal.

Simulating points on sphere

I’ll illustrate this by choosing pairs of points at random on a unit sphere in 100 dimensions and computing their dot product. This gives the cosine similarity of the points, the cosine of the angle between them. I’ll repeat this process 1,000 times and look at the average cosine similarity.

    from scipy.stats import norm
    import numpy as np

    np.random.seed(20230809)

    numruns = 1000
    avg = 0
    dim = 100
    for _ in range(numruns):
        x = norm(0,1).rvs(size = dim)
        y = norm(0,1).rvs(size = dim)
        avg += np.dot(x, y)/(np.dot(x,x)*np.dot(y,y))**0.5
     print(avg/numruns)

Here we use the standard method of generating random points on a sphere: generate each component as a normal (Gaussian) random variable, then divide by the length. The code above combines the normalization with finding the dot product.

The code above prints 0.00043. The average dot product between vectors is very small, i.e. the vectors are nearly orthogonal.

Simulating points in orthant

Next let’s look at one orthant of the sphere. (“Orthant” is the generalization to higher dimensions of quadrant in two dimensions and octant in three dimensions.) We’ll generate random samples from the first orthant by first generating random samples from the whole sphere, then taking the absolute value to fold all the points from the other orthants into the first orthant.

The only necessary change to the code is to put abs() around the code generating the random points.

    x = abs(norm(0,1).rvs(size = dim))
    y = abs(norm(0,1).rvs(size = dim))

When we run the code again it prints 0.639. This is the average value of the cosine of the angle between points, so the average angle is about 0.877 radians or about 50°.

This isn’t quite the right thing to do. We should convert to angles first, then take the average. But that leads to essentially the same result.

If my pencil-and-paper calculations are correct, the exact expected dot product value approaches 2/π = 0.6366 as dimension goes to infinity, and the value of 0.639 above would be consistent with this.

More high-dimensional geometry posts

Ancient estimate of π and modern numerical analysis

A very crude way to estimate π would be to find the perimeter of squares inside and outside a unit circle.

Circle with inscribed and circumscribed squares

The outside square has sides of length 2, so 2π < 8. The inside square has sides of length 2/√2, so 8/√2 < 2π. This tells us π is between 2.82 and 4. Not exciting, but it’s a start.

Finding the perimeter of an inscribed hexagon and a circumscribed hexagon would be better. Finding the perimeter of inscribed and circumscribed n-gons for larger values of n would be even better.

Archimedes and successors

Archimedes (287–212 BC) estimated π using inscribed and circumscribed 96-gons.

Aryabhata (476–450 AD) used Archimedes idea with n = 384 (3 × 27) and Vièta (1540–1603) used n = 39,3216 (3 × 217).

Accuracy

The perimeter of a regular n-gon inscribed in the unit circle is

Tn = 2n sin(π/n)

and the perimeter of a regular n-gon circumscribed on the unit circle is

Un = 2n tan(π/n).

We look at these formulas and see a way to calculate Tn and Un by plugging π/n into functions that calculate sine and cosine. But Archimedes et al had geometric ways of computing Tn and Un and used these to estimate π.

We can use the formulas above to evaluate the bounds above.

Archimedes: 3.1410 < π < 3.1427.

Aryabhata: 3.14155 < π < 3.14166.

Vièta: 3.14159265355 < π < 3.1415926536566.

Extrapolation

In 1654, Huygens discovered that algebraic combinations of the perimeter estimates for an n-gon and an (n/2)-gon could be better than either. So, for example, while Archimedes estimate based on a 96-gon is better than an estimate based on a 48-gon, you could combine the two to do better than both.

His series were

Sn = (4TnTn/2)/3

and

Vn = (2Un + Tn/2)/3

Huygens did not prove that Sn and Vn were lower and upper bounds on 2π, but they are, and are much better than Tn and Un. Using a 48-gon and a 96-gon, for example, Huygens could get twice as many correct digits as Archimedes did with a 96-gon alone.

From a modern perspective, Huygens’ expressions can be seen as a first step toward what Richardson called “the deferred approach to the limit” in 1927, or what is often called Richardson extrapolation. By expanding the expressions for Tn and Un as Taylor series in h = 1/n, we can take algebraic combinations of expressions for varying values of h that eliminate small powers of h, increasing the order of the approximation.

Huygens came up with expressions that canceled out terms of order h2 and leaving terms of order h4. (The odd-order terms were already missing from sin(h)/h and tan(h)/h.)

Milne (1896–1950) came up with a systematic approach that could, for example, use polygons with nn/2, and n/4 sides to come up with an approximation for π with order h6 and to form approximation with order equal to any even power of h. By making clever use of perimeter values known to Archimedes, Milne could obtain more accurate results than those that Aryabhata and Vièta obtained.

Richardson saw a general pattern in the work of Milne and others. We take expressions involving limits as h goes to 0, i.e. derivatives, and to do some manipulation before letting some h go to 0. Taking the derivative would be to take a limit immediately. But if we let some terms cancel each other out before taking the limit, we can get approximations that converge more quickly.

Richardson developed this idea used it as a way to solve a variety of problems more efficiently, such as numerical integration and numerically solving differential equations.

References

[1] L. F. Richardson. The deferred approach to the limit. I. Single lattice, Philos. Trans. Roy. Soc. London Ser. A, 226, pp. 299-349 (1927)

[2] D. C. Joyce. Survey of Extrapolation Process in Numerical Analysis by SIAM Review , Oct., 1971, Vol. 13, No. 4, pp. 435–490

 

Relating perimeter, inner radius, outer radius, and sides of a triangle

Suppose a triangle T has sides a, b, and c.

Let s be the semi-perimeter, i.e. half the perimeter.

Let r be the inner radius, the radius of the largest circle that can fit inside T.

Let R be the outer radius, the radius of the smallest circle that can enclose T.

Then three simple equations relate a, b, c, s, r, and R.

\begin{align*} a + b + c &= 2s \\ ab + bc + ac &= s^2 + r^2 +4rR \\ abc &= 4Rrs \end{align*}

Given a, b, and c, use the first equation to solve for s, then the third equation for Rr, then the second for r, then go back to the last equation to find R.

Given s, r, and R, you can calculate the right hand sides of the three equations above, which are the coefficients in a cubic equation for the sides a, b, and c.

x^3 - (2s)x^2 + (s^2 + r^2 + 4Rr)x -(4Rrs)= 0

Note that this last statement is not about triangles per se. It’s a consequence of

(x-a)(x-b)(x-c) = x^3 - (a+b+c)x^2 + (ab + bc + ac) -abc

which would be true even if a, b, and c were not the sides of a triangle. But since they are sides of a triangle here, the coefficients can be interpreted in terms of geometry, namely in terms of perimeter, inner radius, and outer radius.

Related posts

Multiplication via parabola

Here’s a graphical way to multiply two positive numbers a and b using the parabola y = x².

  1. Start at the origin, move a units to the left, then go up vertically to the parabola, and draw a point.
  2. Go back to the origin, move b units to the right, go up vertically to the parabola, and draw another point.
  3. Connect the points and see where they cross the y-axis. That point is ab.

Here’s an example multiplying 3 and 5.

Multiplying numbers geometrically via a parabola

Here’s why this works. The slope of the line is the change in y over the change in x which is

m = (b² − a²)/(b − (−a)) = ba.

Use the equation of a line

yy0 = m(xx0)

with x0 = b and y0 = b² to get

y − b² = (ba)(xb).

Stick in x = 0 and you get y = ab.

Equation of an ellipse or hyperbola through five points

Yesterday I wrote about the equation of a circle through three points. This post will be similar, looking at the equation of an ellipse or hyperbola through five points. As with the earlier post, we can write down an elegant equation right away. Making the equation practical will take more work.

The equation of a general conic section through five points

\begin{vmatrix} x^2 & y^2 & xy & x & y & 1 \\ x_1^2 & y_1^2 & x_1 y_1 & x_1 & y_1 & 1 \\ x_2^2 & y_2^2 & x_2 y_2 & x_2 & y_2 & 1 \\ x_3^2 & y_3^2 & x_3 y_3 & x_3 & y_3 & 1 \\ x_4^2 & y_4^2 & x_4 y_4 & x_4 & y_4 & 1 \\ x_5^2 & y_5^2 & x_5 y_5 & x_5 & y_5 & 1 \\ \end{vmatrix} = 0

Direct but inefficient solution

The most direct (though not the most efficient) way to turn this into an equation

Ax^2 + By^2 + Cxy + Dx + Ey + F = 0

is to expand the determinant by minors of the first row.

A = \begin{vmatrix} y_1^2 & x_1 y_1 & x_1 & y_1 & 1 \\ y_2^2 & x_2 y_2 & x_2 & y_2 & 1 \\ y_3^2 & x_3 y_3 & x_3 & y_3 & 1 \\ y_4^2 & x_4 y_4 & x_4 & y_4 & 1 \\ y_5^2 & x_5 y_5 & x_5 & y_5 & 1 \\ \end{vmatrix} ,\, B = - \begin{vmatrix} x_1^2 & x_1 y_1 & x_1 & y_1 & 1 \\ x_2^2 & x_2 y_2 & x_2 & y_2 & 1 \\ x_3^2 & x_3 y_3 & x_3 & y_3 & 1 \\ x_4^2 & x_4 y_4 & x_4 & y_4 & 1 \\ x_5^2 & x_5 y_5 & x_5 & y_5 & 1 \\ \end{vmatrix}

C = \begin{vmatrix} x_1^2 & y_1^2 & x_1 & y_1 & 1 \\ x_2^2 & y_2^2 & x_2 & y_2 & 1 \\ x_3^2 & y_3^2 & x_3 & y_3 & 1 \\ x_4^2 & y_4^2 & x_4 & y_4 & 1 \\ x_5^2 & y_5^2 & x_5 & y_5 & 1 \\ \end{vmatrix} ,\, D = - \begin{vmatrix} x_1^2 & y_1^2 & x_1 y_1 & y_1 & 1 \\ x_2^2 & y_2^2 & x_2 y_2 & y_2 & 1 \\ x_3^2 & y_3^2 & x_3 y_3 & y_3 & 1 \\ x_4^2 & y_4^2 & x_4 y_4 & y_4 & 1 \\ x_5^2 & y_5^2 & x_5 y_5 & y_5 & 1 \\ \end{vmatrix}

 E = \begin{vmatrix} x_1^2 & y_1^2 & x_1 y_1 & x_1 & 1 \\ x_2^2 & y_2^2 & x_2 y_2 & x_2 & 1 \\ x_3^2 & y_3^2 & x_3 y_3 & x_3 & 1 \\ x_4^2 & y_4^2 & x_4 y_4 & x_4 & 1 \\ x_5^2 & y_5^2 & x_5 y_5 & x_5 & 1 \\ \end{vmatrix} ,\, F = - \begin{vmatrix} x_1^2 & y_1^2 & x_1 y_1 & x_1 & y_1\\ x_2^2 & y_2^2 & x_2 y_2 & x_2 & y_2\\ x_3^2 & y_3^2 & x_3 y_3 & x_3 & y_3\\ x_4^2 & y_4^2 & x_4 y_4 & x_4 & y_4\\ x_5^2 & y_5^2 & x_5 y_5 & x_5 & y_5\\ \end{vmatrix}

More efficient solution

It would be much more efficient to solve a system of equations for the coefficients. Since the determinant at the top of the page is zero, the columns are linearly dependent, and some linear combination of those columns equals the zero vector. In fact, the coefficients of such a linear combination are the coefficients A through F that we’re looking for.

We have an unneeded degree of freedom since any multiple of the coefficients is valid. If we divide every coefficient by A, then the new leading coefficient is 1. This gives us five equations in five unknowns.

\begin{bmatrix} y_1^2 & x_1 y_1 & x_1 & y_1 & 1 \\ y_2^2 & x_2 y_2 & x_2 & y_2 & 1 \\ y_3^2 & x_3 y_3 & x_3 & y_3 & 1 \\ y_4^2 & x_4 y_4 & x_4 & y_4 & 1 \\ y_5^2 & x_5 y_5 & x_5 & y_5 & 1 \\ \end{bmatrix} \begin{bmatrix} B \\ C \\ D \\ E \\ F \\ \end{bmatrix} = - \begin{bmatrix} x_1^2 \\ x_2^2 \\ x_3^2 \\ x_4^2 \\ x_5^2 \\ \end{bmatrix}

Python example

The following code will generate five points and set up the system of equations Ms = b whose solution s will give our coefficients.

    import numpy as np

    np.random.seed(20230619)

    M = np.zeros((5,5))
    b = np.zeros(5)
    for i in range(5):
        x = np.random.random()
        y = np.random.random()    
        M[i, 0] = y**2
        M[i, 1] = x*y
        M[i, 2] = x
        M[i, 3] = y
        M[i, 4] = 1
        b[i] = -x*x

Now we solve for the coefficients.

    s = np.linalg.solve(M, b)
    A, B, C, D, E, F = 1, s[0], s[1], s[2], s[3], s[4]

Next we verify that the points we generated indeed lie on the conic section whose equation we just found.

    for i in range(5):
        x = M[i, 2]
        y = M[i, 3]
        assert( abs(A*x*x + B*y*y + C*x*y + D*x + E*y + F) < 1e-10 )

Ellipse or hyperbola?

It turns out that the five points generated above lie on a hyperbola.

If we had set the random number generator seed to 20230620 we would have gotten an ellipse.

When we generate points at random, as we did above, we will almost certainly get an ellipse or a hyperbola. Since we’re generating pseudorandom numbers, there is some extremely small chance the generated points would lie on a line or on a parabola. In theory, this almost never happens, in the technical sense of “almost never.”

If 4ABC² is positive we have a an ellipse. If it is negative we have a hyperbola. If it is zero, we have a parabola.

Related posts

A parabola, a triangle, and a circle

Let P be a parabola. Draw tangents to P at three different points. These three lines intersect to form a triangle.

Theorem: The circumcircle of this triangle passes through the focus of P. [1]

In the image above, the dashed lines are tangents and the black dot is the focus of the parabola.

(See this post for an explanation of what focus and directrix mean.)

By the way, creating the image above required finding the circle through the three intersection points. See the previous post for how to find such a circle.

[1] Ross Honsberger. Episodes in Nineteenth and Twentieth Century Euclidean Geometry. Mathematical Association of America. 1995. Page 47.

Equation of a circle through three points

Three given points on a circle with unknown center

A few months ago I wrote about several equations that come up in geometric calculations that all have the form of a determinant equal to 0, where one column of the determinant is all 1s.

The equation of a circle through three points (x1, y1), (x2, y2), and (x3, y3) has this form:

\begin{vmatrix} x^2 + y^2 & x & y & 1 \\ x_1^2 + y_1^2 & x_1 & y_1 & 1 \\ x_2^2 + y_2^2 & x_2 & y_2 & 1 \\ x_3^2 + y_3^2 & x_3 & y_3 & 1 \\ \end{vmatrix} = 0

This is elegant equation, but it’s not in a form that’s ready to program. This post will start with the equation above and conclude with code.

Let Mij be the determinant of the minor of the (i, j) element of the matrix above. That is, Mij is the determinant of the 3 × 3 matrix obtained by deleting the ith row and jth column. Then we can expand the determinant above by minors of the last column to get

M_{11}(x^2 + y^2) - M_{12}x + M_{13}y - M_{14} = 0

which we can rewrite as

x^2 + y^2 -\frac{M_{12}}{M_{11}} x + \frac{M_{13}}{M_{11}} y - \frac{M_{14}}{M_{11}} = 0.

A circle with center (x0, y0) and radius r0 has equation

(x - x_0)^2 + (y - y_0)^2 - r_0^2 = 0.

By equating the coefficients x and y in the two previous equations we have

\begin{align*} x_0 &= \phantom{-}\frac{1}{2} \frac{M_{12}}{M_{11}} \\ y_0 &= -\frac{1}{2} \frac{M_{13}}{M_{11}} \\ \end{align*}

Next we expand the determinant minors in more explicit form.

\begin{align*} M_{11} &= \begin{vmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \\ \end{vmatrix} \\ &= (x_2y_3 - x_3y_2) - (x_1y_3 - x_3y_1) + (x_1y_2 - x_2y_1) \\ &= (x1y_2 + x_2y_3 + x_3y_1) - (x_2y_1 + x_3y_2 + x_1y_3) \end{align*}

Having gone to the effort of finding a nice expression for M11 we can reuse it to find M12 by the substitution

x_i \to s_i \equiv x_i^2 + y_i^2

as follows:

\begin{align*} M_{12} &= \begin{vmatrix} x_1^2 + y_1^2 & y_1 & 1 \\ x_2^2 + y_2^2 & y_2 & 1 \\ x_3^2 + y_3^2 & y_3 & 1 \\ \end{vmatrix} = \begin{vmatrix} s_1^2 & y_1 & 1 \\ s_2^2 & y_2 & 1 \\ s_3^2 & y_3 & 1 \\ \end{vmatrix} \\ &= (s_1y_2 + s_2y_3 + s_3y_1) - (s_2y_1 + s_3y_2 + s_1y_3) \end{align*}

And we can find M13 from M12 by reversing the role of the xs and ys.

\begin{align*} M_{13} &= \begin{vmatrix} x_1^2 + y_1^2 & x_1 & 1 \\ x_2^2 + y_2^2 & x_2 & 1 \\ x_3^2 + y_3^2 & x_3 & 1 \\ \end{vmatrix} = \begin{vmatrix} s_1^2 & x_1 & 1 \\ s_2^2 & x_2 & 1 \\ s_3^2 & x_3 & 1 \\ \end{vmatrix} \\ &= (s_1x_2 + s_2x_3 + s_3x_1) - (s_2x_1 + s_3x_2 + s_1x_3) \end{align*}

Now we have concrete expressions for M11, M12, and M13, and these give us concrete expressions for x0 and y0.

To find the radius of the circle we simply find the distance from the center (x0, y0) to any of the points on the circle, say (x1, y1).

The following Python code incorporates the derivation above.

    def circle_thru_pts(x1, y1, x2, y2, x3, y3):
        s1 = x1**2 + y1**2
        s2 = x2**2 + y2**2
        s3 = x3**2 + y3**2
        M11 = x1*y2 + x2*y3 + x3*y1 - (x2*y1 + x3*y2 + x1*y3)
        M12 = s1*y2 + s2*y3 + s3*y1 - (s2*y1 + s3*y2 + s1*y3)
        M13 = s1*x2 + s2*x3 + s3*x1 - (s2*x1 + s3*x2 + s1*x3)
        x0 =  0.5*M12/M11
        y0 = -0.5*M13/M11
        r0 = ((x1 - x0)**2 + (y1 - y0)**2)**0.5
        return (x0, y0, r0)

The code doesn’t use any loops or arrays, so it ought to be easy to port to any programming language.

For some inputs the value of M11 can be (nearly) zero, in which case the three points (nearly) lie on a line. In practice M11 can be nonzero and yet zero for practical purposes and so that’s something to watch out for when using the code.

The determinant at the top of the post is zero when the points are colinear. There’s a way to make rigorous the notion that in this case the points line on a circle centered at infinity.

Update: Equation of an ellipse or hyperbola through five points

Similar triangles and complex numbers

Suppose the vertices of two triangles are given by complex numbers a, b, c and x, y, z. The two triangles are similar if

\left| \begin{matrix} 1 & 1 & 1 \\ a & b & c \\ x & y & z \\ \end{matrix} \right| = 0

This can be found in G. H. Hardy’s classic A Course of Pure Mathematics. It’s on page 93 in the 10th edition.

Corollary

The theorem above generalizes a result from an earlier post giving a test for whether points determine an equilateral triangle. A triangle ABC is equilateral if it is similar to BCA. For triangles to be similar, corresponding vertex angles must be equal. If you can rotate the order of the angles in one triangle and the two triangles remain similar, the angles must all be equal.

Proof

Hardy’s proof of the equation above is direct and terse. I’ll quote it here in its entirety.

The condition required is that

\overline{AB}\,/\,\overline{AC} = \overline{XY}\,/\,\overline{XZ}

(large letters denoting the points whose arguments are the corresponding small letters), or

(b-a)(c-a) = (y-x)/(z-x)

which is the same as the condition given.

Commentary

To expand on this a little bit, Hardy is saying that we need to show that if we take two corresponding vertices, A and X, the ratio of the lengths of corresponding sides is the same in both triangles. If memory serves, this is called the “SAS postulate.”

The determinant condition implies the last equation in Hardy’s proof. This is not obvious, but it’s easy to work out with pencil and paper. The reader may be left thinking “OK, I can see it’s true, but why would anyone think of that?” A more motivated derivation would require more background than the reader of the book is presumed to have.

You might object that Hardy’s last equation is a ratio of complex numbers, not a ratio of lengths. But you can take the absolute value of both sides. And the absolute value of a fraction is the absolute value of the numerator divided by the absolute value of the denominator. And now you’re taking a ratio of side lengths. In symbols,

\frac{b-a}{c-a} = \frac{y-x}{z-x} \implies \frac{|b-a|}{|c-a|} = \frac{|y-x|}{|z-x|}

Related posts

Comparing approximations for ellipse perimeter

This post will compare the accuracy of approximations for the perimeter of an ellipse.

The exact perimeter is given in terms of an elliptic integral. (That’s where the elliptic integrals  gets their name.) And so an obvious way to approximate the perimeter would be to expand the elliptic integral in a power series. Unfortunately this series converges slowly for highly eccentric ellipses.

Ernst Kummer found a series that converges more quickly [1]. Let a and b be the semi-major and semi-minor axes of an ellipse and let

 x = \left(\frac{a-b}{a+b} \right )^2

Then Kummer’s series is

\pi(a+b)\left(1 + \frac{x}{4} + \frac{x^2}{64} + \frac{x^3}{256} + \frac{25x^4}{16384} + \cdots \right )

The following plot compares the accuracy of truncating Kummer’s series after powers of x equal to 1 through 4.

The noise at the bottom of the chart comes from the limits of machine precision. For moderate values of eccentricity, three terms of Kummer’s series gives an approximation so accurate that all the error is coming from the limitations of floating point resolution.

I’ve written about approximations for the perimeter of an ellipse before. The approximation based on the 3/2 norm is about as accurate as Kummer’s series truncated after the x term. The error in Ramanujan’s approximation is between that of Kummer’s series with the cubic and quartic terms.

For the orbit of Pluto (e = 0.2488) the first-order Kummer approximation is good to 9 figures and the higher-order approximations are essentially good to machine precision.

For the orbit of Haley’s comet (e = 0.9671) the successive Kummer approximations are good to 2, 3, 4, and 5 figures.

[1] Carl E. Linderholm and Arthur C. Segal. An Overlooked Series for the Elliptic Perimeter. Mathematics Magazine, Vol. 68, No. 3 (Jun., 1995), pp. 216–220.