Perimeter of Ellipse

[Gerson W. Barbosa]

(01-19-2020 11:00 PM)Albert Chan Wrote:  If you have AGM, this recursive formula is better (based on AGM2 method):

\(\large p(a,b) = 2× \left( p({a+b \over 2}, \sqrt{ab}) - {\pi a b \over AGM(a,b)}\right)\)

In other words, calculate ellipse perimeter from a much less eccentric ellipse.

You can do this recursively, but then you might as well use AGM2.
Note: AGM only needed to calculate once, since AGM(a,b) = AGM((a+b)/2, √(ab)) = ...

Code:

from gmpy2 import *
get_context().precision = 160
pi = const_pi()

def pade22(a,b):
    h = ((a-b)/(a+b))**2
    return pi*(a+b) * ((-21*h-48)*h+256) / ((3*h-112)*h+256)

>>> a=mpfr('5906376272')
>>> b=mpfr('5720637952.8')
>>> print pade22(a,b)
36529672878.015838406030163739762066543303728710139

>>> p = pade22((a+b)/2, sqrt(a*b))
>>> print 2*(p - pi*a*b / agm(a,b))
36529672878.015838406035141932308423167594568357066

AGM improved pade22 doubled accuracy, from 22 to 45 digits.

That’s it! Thanks!

I’ve used it on my other example, the much more eccentric orbit of Halley’s comet ( a = 2667950000 km; b = 678282900 km ). Just one call to the function, as in your example

p(a, b) = π{(a + b)[h²(3h² - 136) + 320] + [a b/(a + b)](96h² - 256)}/[h²(3h² - 112) + 256] (*)
where h = [(a - b)/(a + b)]

and the error was reduced from 12315.162 km to 0.555 meters!

(*) When comparing the arithmetic and harmonic means with the equivalent radius, I came up with Peano approximation, of which that’s a refinement.

P.S.:

The above approximation can be simplified to

p(a, b) = π(a + b)[1 - (24h - 64)/(3h + 256/h - 112)]

where h = [(a - b)/(a + b)]²