Round the ellipse and back again

Here’s my best approximation of the perimeter of an ellipse.

One thing you might notice is most of the approximations Matt presents get progressively worse as a/b increases. Matt’s best one’s percent error does bounce around for a bit but then it takes off too for a/b > 3.6.

But that’s sort of the nature of approximations, when what you’re trying to approximate increases without bound: If your approximation doesn’t capture the exact behavior (and if it does, it’s not really an approximation, is it?) then it’ll necessarily diverge from the correct value by more and more. Really the best you can hope for is to get a good approximation within some bounds, knowing it’ll go increasingly bad once you’re far enough outside those bounds.

So I thought I’d see if I could get the divergence to occur much further out than a/b > 3.6, and it seemed to me if I wanted a better (in the sense of more accurate, not more beautiful) formula than Ramanujan’s, I should start with Ramanujan’s and fix it up.

Ramanujan’s approximation is:

p_{ram} = \pi(3(a+b)-\sqrt{(3a+b)(a+3b)})

I cracked open a spreadsheet and found that if \delta = p_{true}-p_{ram} then for b = 1 and 3 < a < 20, \delta^{0.62} is approximately linear in a. After a bunch of manipulation I got the corrected approximation

p_{corr} = \pi(3(a+b)-\sqrt{(3a+b)(a+3b)} + qb(a/b-3)^r)

where q = 0.0004714 and r = 1.6156.

This formula doesn’t work for a/b < 3, but if my numbers for the “exact” perimeter (series sum to 20 terms) are correct, then for 3 \le a/b < 17.5, the average absolute percent error is about 0.0012%. Beyond that it starts to diverge. Here’s the percent error, in red, along with the percent error for Ramanujan’s formula, in blue.

And here’s a Python function to calculate the corrected approximation (for a/b > 3, otherwise it returns Ramanujan’s approximation):


from sys import argv
from math import pi

def ellper (a, b):
    q = 0.0004714
    r = 1.6156

    pram = pi*(3*(a+b)-((3*a+b)*(a+3*b))**0.5)
    return pram + pi*(q*b*(a/b-3)**r) if a/b >= 3 else pram

a = float(argv[1])
b = float(argv[2])

print ellper (a, b)

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s