Perimeter of Ellipse

[Joe Horn]

Although simple formulae for the perimeter of an ellipse exist, they are only approximations. Exact formulae are complicated. The following program uses a converging iteration technique instead. Ported from a QBASIC program by Gérard P. Michon.

Syntax: EllipsePerimeter(a,b), where a & b are the lengths of the axes (order doesn't matter)
Output: perimeter of ellipse

Examples:
EllipsePerimeter(0.5, 0.5) --> 3.14159265359
EllipsePerimeter(3, 4) --> 22.1034921607

Code:

gk(h);
cayley(x);

EXPORT EllipsePerimeter(a,b)
BEGIN
LOCAL x,h,P;
a:=ABS(a);
b:=ABS(b);
IF a < b THEN x:=a; a:=b; b:=x; END;
IF b < 0.28*a THEN
 P := 4*a*cayley((b/a)^2);
ELSE
 h := ((a-b)/(a+b))^2;
 P := pI*(a+b)*gk(h);
END;
RETURN P;
END
;

gk(h)
BEGIN
LOCAL z,x,n;
z := 0; x := 1; n := 0;
WHILE z + x <> z
 DO
 n := n + 1;
 x := h * x * ((n-1.5)/n)^2;
 z := z + x;
END;
RETURN 1 + z;
END;

cayley(x)
BEGIN
LOCAL y,t,n,z,u,v,w;
y := LOG(16 / x) - 1;
t := x / 4;
n := 1;
z := 0;
u := t * y;
v := (n - .5) / n;
w := .5 / ((n - .5) * n);
WHILE z <> z + u
 DO
 z := z + u;
 n := n + 1;
 t := x * t * v;
 v := (n - .5) / n;
 t := t * v;
 y := y - w;
 w := .5 / ((n - .5) * n);
 y := y - w;
 u := t * y;
END;
RETURN 1 + z;
END;

[Wes Loewer]

I was going to suggest the following:

Code:

EXPORT EllipsePerimeter(a,b)
BEGIN
 RETURN ∫(√(((a*COS(T))^2+(b*SIN(T))^2)),T,0,2*π);
END;

but it seems that PPL doesn't like local variables inside of an integral.

Here's a work-around:

Code:

EXPORT EllipsePerimeter(a,b)
BEGIN
 A:=a;
 B:=b;
 RETURN ∫(√(((A*COS(T))^2+(B*SIN(T))^2)),T,0,2*π);
END;

[Wes Loewer]

or this work-around:

Code:

EXPORT EllipsePerimeter(a,b)
BEGIN
 RETURN approx(int(√(((a*COS(T))^2+(b*SIN(T))^2)),T,0,2*π));
END;

[TASP]

Did not know that about perimeters of ellipses.

Learn something here everyday.

[parisse]

I guess the perimeter is computed using the arithmetic geometric mean, correct?

[SlideRule]

(03-05-2016 04:19 PM)Joe Horn Wrote:  Although simple formulae for the perimeter of an ellipse exist, they are only approximations. Exact formulae are complicated. The following program uses a converging iteration technique instead.

When I was employed by an Engineering firm specializing in Bridge Design, we used elliptical substructure rather than rectangular. The design was easy enough BUT the insitu replication of that design was elusive, so we substituted cirular arcs (tangent) for the desired elipse. The desired substitute was reflective about both inplane axis, so only 1/4 of the form need be designed. Seldom were more the 5 contiguous arcs needed while 3 arcs satisfied most requirements. The circle(s) approximations functioned well within the desired strutural requirements while easing the construction layout.

BEST!
SlideRule

[Joe Horn]

(03-06-2016 02:16 PM)Wes Loewer Wrote:  or this work-around:

Code:

EXPORT EllipsePerimeter(a,b)
BEGIN
 RETURN approx(int(√(((a*COS(T))^2+(b*SIN(T))^2)),T,0,2*π));
END;

Your program above (unlike mine) always agrees with Wolfram Alpha's perimeter of an ellipse, even for very eccentric ellipses. Congratulations!

[parisse]

I can add this functionnality to perimeter, using the approximate value of the integral
http://www-fourier.ujf-grenoble.fr/%7epa...(-1,1,2))&
(make sure to clear your browser cache if you have already tried Xcas offline).

[Albert Chan]

From HP Forum Archive 21, message 5, and some optimizations, I get:

Quote:local sqrt, pi = math.sqrt, math.pi

function ellipsePerimeter(a, b) -- agm2 method
    local s = 0.5 * (a*a+b*b)
    local t = 0.25
    while true do
        local k = b - a
        local s1 = s - t*k*k
        if s == s1 then return 2*pi*s/(a+0.5*k) end
        b = sqrt(a*b)
        a = a + 0.5*k
        s = s1
        t = t + t
    end
end

Example:

lua> ellipsePerimeter(1, 2) --> 9.688448220547675
lua> ellipsePerimeter(2, 3) --> 15.86543958929059
lua> ellipsePerimeter(3, 4) --> 22.103492160709504

[ggauny@live.fr]

(03-07-2016 03:34 PM)Joe Horn Wrote:  

(03-06-2016 02:16 PM)Wes Loewer Wrote:  or this work-around:

Code:

EXPORT EllipsePerimeter(a,b)
BEGIN
 RETURN approx(int(√(((a*COS(T))^2+(b*SIN(T))^2)),T,0,2*π));
END;

Your program above (unlike mine) always agrees with Wolfram Alpha's perimeter of an ellipse, even for very eccentric ellipses. Congratulations!

[Gerson W. Barbosa]

(03-24-2019 12:42 PM)Albert Chan Wrote:  From HP Forum Archive 21, message 5, and some optimizations, I get:

Quote:local sqrt, pi = math.sqrt, math.pi

function ellipsePerimeter(a, b) -- agm2 method
    local s = 0.5 * (a*a+b*b)
    local t = 0.25
    while true do
        local k = b - a
        local s1 = s - t*k*k
        if s == s1 then return 2*pi*s/(a+0.5*k) end
        b = sqrt(a*b)
        a = a + 0.5*k
        s = s1
        t = t + t
    end
end

Example:

lua> ellipsePerimeter(1, 2) --> 9.688448220547675
lua> ellipsePerimeter(2, 3) --> 15.86543958929059
lua> ellipsePerimeter(3, 4) --> 22.103492160709504

On Free42 using Hugh Steers’s original algorithm I get

1 ENTER 2 XEQ PEL -> 9.688448220547676198428503196391833

2 ENTER 3 XEQ PEL -> 15.86543958929058979133166302778308

3 ENTER 4 XEQ PEL -> 22.10349216070950504528558646387247

Code:



00 { 75-Byte Prgm }
01▸LBL "PEL"
02 RCL÷ ST Y
03 RCL ST X
04 X↑2
05 1
06 STO 01
07 STO 02
08 +
09 2
10 ÷
11 STO 03
12 SIGN
13▸LBL 00
14 RCL+ ST Y
15 2
16 ÷
17 X<> 02
18 RCL 02
19 X=Y?
20 GTO 01
21 R↓
22 RCL× ST Y
23 LASTX
24 X<>Y
25 SQRT
26 X<> ST Z
27 -
28 2
29 ÷
30 X↑2
31 RCL× 01
32 STO- 03
33 R↓
34 2
35 STO× 01
36 CLX
37 RCL 02
38 GTO 00
39▸LBL 01
40 STO÷ 03
41 R↑
42 PI
43 ×
44 STO+ ST X
45 RCL× 03
46 END

On the wp34s, for low-eccentricity ellipses, the following approximation might be handy, since no programming is required:

2\(\pi\)[3*agm(a,b) - 2√(a*b)]

The error in the lenght of the orbit of Pluto should be less than 5 meters when using that approximation, if my calculation is correct.

[Valentin Albillo]

(12-04-2019 05:30 PM)Gerson W. Barbosa Wrote:  The error in the lenght of the orbit of Pluto should be less than 5 meters when using that approximation, if my calculation is correct.

Your calculation may be correct but I seriously doubt the orbital parameters for Pluto are actually known to such precision, most especially since the orbit's semi-major axis is some 5.90638 billion km, i.e.: 5.90638 trillion meters.

An error "less than 5 meters" in such huge value seems unrealistic from a purely physical point of view. From a purely mathematical point of view, that's another matter and your error estimation may be correct, though again the orbit's eccentricity is not known to high accuracy either (~0.2488, quite large when compared with the much more circular orbit of Earth, say).

Regards.
V.

[Gerson W. Barbosa]

(12-04-2019 06:28 PM)Valentin Albillo Wrote:  An error "less than 5 meters" in such huge value seems unrealistic from a purely physical point of view. From a purely mathematical point of view, that's another matter and your error estimation may be correct, though again the orbit's eccentricity is not known to high accuracy either (~0.2488, quite large when compared with the much more circular orbit of Earth, say).

Hello, Valentin,

This has been done only to assess the quality of the approximation, if any. Anyway, despite the uncertainties involved, the order of magnitude of the difference should be correct.

Some doublechecking:

semi-major axis: a = 5906376272 km
eccentricity: e = 0.24880766 ->
semi-minor axis: b = 5720637952.8 km,

for which the Free42 program returns the exact result

36529672878.01583840603514193230844 km

Now, on the wp34s,

5906376272 ENTER

5720637952.8

AGM ->

5813136193.07

3 * 5906376272 ENTER

5720637952.8 * √ 2 * -

2 * π *

-> 36529672878.01109446182 km

Difference:

0.00474394421514193230844 km

or

4.74 m

———-

This first Ramanujan approximation brings the error down to under one meter:

π[3(a + b) - √((3a + b)(a + 3b))]

His second approximation is even better, with an error less than one micrometer:

https://www.johndcook.com/blog/2013/05/0...-ellipse/

Best regards,

Gerson.

[Gerson W. Barbosa]

This is an equivalent Decimal BASIC program of a C-program by Hugh Steers here. The port of this one or the original C-program to the Prime should be straightforward. The computational effort is a fraction of that involved in the computation of the elliptic integral, but that may not make a significant difference in fast calculators like the Prime.

OPTION ARITHMETIC DECIMAL_HIGH
INPUT  PROMPT "a, b: ":ae,be
LET nd = 1000
LET tm = TIME 
LET b = be/ae
LET s = (1 + b*b)/2
LET a = 1 
LET t = 1
DO
   LET a1 = (a + b)/2
   IF a = a1 THEN EXIT DO
   LET c = (a - b)/2
   LET b = SQR(a*b)
   LET a = a1
   LET s = s - t*c*c
   LET t = t + t
LOOP
LET p = 2*PI*ae*s/a
LET r = TIME - tm
LET r$ = STR$(p - INT(p))                           
PRINT
PRINT STR$(INT(p));
PRINT r$(0:1);
FOR i = 2 TO nd + 1
   PRINT r$(i:i);
   IF MOD((i - 1),10) = 0 THEN PRINT " ";
   IF MOD((i - 1),50) = 0 THEN 
      PRINT
      PRINT REPEAT$(" ",1 + LEN(STR$(INT(p))));
   END IF
NEXT i
IF MOD (i - 2,50) <> 0  OR nd = 0 THEN PRINT
PRINT 
PRINT "Runtime: ";
PRINT  USING "0.##": r;
PRINT " seconds"
END


---------------


a, b: 3,4

22.1034921607 0950504528 5586463872 4607782782 8910007059 
   3992632612 5910396946 4959177286 6218309324 0390352083 
   5815953162 7467734138 2129000335 3896736622 8442994174 
   6192287204 5014970302 5317033312 3929977604 0095275908 
   3220930141 4469558250 9419563787 3762863238 9667370976 
   2571981363 4931419230 9352895844 1906042118 2907229022 
   9686418117 2635887950 2052296943 0126681600 8830044762 
   7954774028 4076289798 3279618691 9046286983 5594397451 
   2356208773 7408645009 3309824766 6784863302 2726732463 
   0516083153 4972329406 0801844377 6345370790 6974256179 
   0156140823 0342896177 6736816934 6443281912 6179659392 
   2331855348 5127216985 8740979075 6632664289 7373983969 
   0524635361 8083080355 2201668145 2949128953 8470111771 
   1967401421 1827316426 8227257828 0215582400 2244838587 
   9322889687 7768410992 9363265305 7955533884 5140145875 
   1572214188 4766228802 8889719384 4198561575 6084702523 
   4710513780 8869298904 8167827143 0007524734 3709238837 
   9951007032 1404352133 5658034012 3750989012 2020953649 
   3147176616 2140936028 0496095218 3637964950 6346311041 
   3538836238 7660229689 7814732208 7377779791 016495137 
   
Runtime: 0.05 seconds

[Gerson W. Barbosa]

(12-04-2019 10:27 PM)Gerson W. Barbosa Wrote:  ———-

This first Ramanujan approximation brings the error down to under one meter:

π[3(a + b) - √((3a + b)(a + 3b))]

His second approximation is even better, with an error less than one micrometer²:

https://www.johndcook.com/blog/2013/05/0...-ellipse/

Instead of

p ~ 2π[3*agm(a,b) - 2√(a*b)]

I will suggest a better approximation involving the arithmetic-geometric mean and the geometric mean:

p ~ 2π{a + b - a*b/agm(a,b) - 2[agm(a,b) - √(a*b)]}

Let’s use it to compute the length of the orbit of Pluto¹:

——

wp34s program

001 LBL  A
002 ©ENTER
003 STO+  T
004 RCL×  Y
005 √
006 STO  I
007 x⇆  L
008 ⇆  ZYXT
009 AGM
010 STO/  Y
011 RCL- I
012 STO+  X
013 +
014 -
015 #  π
016 ×
017 STO+ X
018 END

——

5906376272 ENTER 5720637952.8 A ->

36529672878.01583848170946635744574

Exact result:

36529672878.01583840603514193230844 km

Difference:

0.000000075674324425137 km,

or

75.674 μm

——


¹ Assuming the parameters are exact and the orbit is perfectly elliptical.

——-
01-03-2020 11:57 PM

PS:

I present another formula involving AGM:

p ~ 2π[agm(a,b)(192(1 - h) - h²) - 128(1 - h)√(ab)]/[64(1 - h) - h²]

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

Error: 32.59 nm (orbit of Pluto)

Further improvement is still possible .



² The actual error produced by Ramanujan’s second formula is slightly less than one nanometer, not one micrometer.

——-
01-05-2020 01:03 AM

PSS:

Just a small refinement – more are still possible – and the overall error in the length of the orbit of Pluto is only 19.34 fm (femtometers!).

p ~ 2π{agm(a,b)[77h² - 768(h - 1)] - [54h² - 512(h - 1)]√(ab)}/[23h² - 256(h - 1)]

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

Example:

a = 5906376272 km

b = 5720637952.8 km

p ~ 36529672878.01583840603514191296851 km

Exact:

p = 36529672878.01583840603514193230844 km

[Gerson W. Barbosa]

(12-07-2019 07:57 PM)Gerson W. Barbosa Wrote:  

(12-04-2019 10:27 PM)Gerson W. Barbosa Wrote:  ———-

This first Ramanujan approximation brings the error down to under one meter:

π[3(a + b) - √((3a + b)(a + 3b))]

His second approximation is even better, with an error less than one micrometer:

https://www.johndcook.com/blog/2013/05/0...-ellipse/

Instead of

p ~ 2π[3*agm(a,b) - 2√(a*b)]

I will suggest a better approximation involving the arithmetic-geometric mean and the geometric mean:

p ~ 2π{a + b - a*b/agm(a,b) - 2[agm(a,b) - √(a*b)]}

Let’s use it to compute the length of the orbit of Pluto¹:

.
.
.



36529672878.01583848170946635744574

Exact result:

36529672878.01583840603514193230844 km

Difference:

0.000000075674324425137 km,

or

75.674 μm

——


¹ Assuming the parameters are exact and the orbit is perfectly elliptical.

We can do better:

\(p\approx \frac{\pi \left [ \left ( 15h\sqrt{1+\frac{3h}{8-3h\sqrt{2}}}-80\right )\left ( a^{2}+\frac{6}{5}ab+b^{2} \right ) +4h\left ( a^{2}+2ab+b^{2} \right ) \right ]}{\left ( 12h\sqrt{1+\frac{3h}{8-3h\sqrt{2}}}+4h-64\right )\left ( a+b \right )}\)


where

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

This one wasn't particularly difficult to find (might explain later).

p ~ 36529672878.01583840603557428740268 km

p = 36529672878.01583840603514193230844 km

Difference: 0.432 nm

Contrary to what is stated above, the difference obtained with Ramanujan's second approximation is less than one nanometer, not one micrometer (actually 0.905 nm), about the double of the above result, but his formula is way more simple.

[Albert Chan]

(12-29-2019 07:58 PM)Gerson W. Barbosa Wrote:  We can do better:

\(p\approx \frac{\pi \left [ \left ( 15h\sqrt{1+\frac{3h}{8-3h\sqrt{2}}}-80\right )\left ( a^{2}+\frac{6}{5}ab+b^{2} \right ) +4h\left ( a^{2}+2ab+b^{2} \right ) \right ]}{\left ( 12h\sqrt{1+\frac{3h}{8-3h\sqrt{2}}}+4h-64\right )\left ( a+b \right )}\)


where

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

This one wasn't particularly difficult to find (might explain later).

We can simplify this a bit.

Let \(\large c = \sqrt{1+\frac{3h}{8-3h\sqrt{2}}}\quad → p ≈ \pi(a+b)\Large\left({3ch^2 + 12(c-1)h - 64 \over 4(3c+1)h-64}\,\right) \)

If \(\large c = 1\), above simplified to: \(\quad\large p ≈ \pi(a+b)\Large\left({3h^2 - 64 \over 16h-64}\,\right) \)

This matched ellipse perimeter Pade [2,1] approximation.

[Albert Chan]

(12-07-2019 07:57 PM)Gerson W. Barbosa Wrote:  Just a small refinement – more are still possible – and the overall error in the length of the orbit of Pluto is only 19.34 fm (femtometers!).

p ~ 2π{agm(a,b)[77h² - 768(h - 1)] - [54h² - 512(h - 1)]√(ab)}/[23h² - 256(h - 1)]

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

Example:

a = 5906376272 km

b = 5720637952.8 km

p ~ 36529672878.01583840603514191296851 km

Exact:

p = 36529672878.01583840603514193230844 km

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.

[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)]²

[Albert Chan]

(01-19-2020 11:00 PM)Albert Chan Wrote:  \(\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.

Amazingly, the less eccentric ellipse have the eccentricity of h

lua> sqrt = math.sqrt
lua> a, b = 2667950000, 678281900 -- numbers from here
lua> = sqrt(1-(b/a)^2)     → e = 0.967142904276921
lua> = (a-b)/(a+b)           → h = 0.5945995852827773 = e'
lua>
lua> a, b = (a+b)/2, sqrt(a*b)
lua> = sqrt(1-(b/a)^2)     → e' = 0.5945995852827773
lua> = (a-b)/(a+b)           → h' = 0.10863394673347321 = e''

Prove: Assume a ≥ b, let a', b' = (a+b)/2, √(ab)

e' = √(1-(b'/a')²) = √(1 - 4ab/(a+b)²) = (a-b)/(a+b) = h