Perimeter of Ellipse
03-05-2016, 04:19 PM
Post: #1
 Joe Horn Senior Member Posts: 1,774 Joined: Dec 2013
Perimeter of Ellipse
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;

<0|ɸ|0>
-Joe-
03-06-2016, 11:55 AM
Post: #2
 Wes Loewer Senior Member Posts: 333 Joined: Jan 2014
RE: Perimeter of Ellipse
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;
03-06-2016, 02:16 PM
Post: #3
 Wes Loewer Senior Member Posts: 333 Joined: Jan 2014
RE: Perimeter of Ellipse
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;
03-06-2016, 02:40 PM
Post: #4
 TASP Senior Member Posts: 401 Joined: Mar 2015
RE: Perimeter of Ellipse
Did not know that about perimeters of ellipses.

Learn something here everyday.

2speed HP41CX,int2XMEM+ZEN, HPIL+DEVEL, HPIL+X/IO, I/R, 82143, 82163, 82162 -25,35,45,55,65,67,70,80
03-06-2016, 06:42 PM
Post: #5
 parisse Senior Member Posts: 1,186 Joined: Dec 2013
RE: Perimeter of Ellipse
I guess the perimeter is computed using the arithmetic geometric mean, correct?
03-07-2016, 01:16 PM (This post was last modified: 03-07-2016 01:16 PM by SlideRule.)
Post: #6
 SlideRule Senior Member Posts: 1,313 Joined: Dec 2013
RE: Perimeter of Ellipse
(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
03-07-2016, 03:34 PM
Post: #7
 Joe Horn Senior Member Posts: 1,774 Joined: Dec 2013
RE: Perimeter of Ellipse
(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!

<0|ɸ|0>
-Joe-
03-09-2016, 08:39 AM (This post was last modified: 03-09-2016 08:42 AM by parisse.)
Post: #8
 parisse Senior Member Posts: 1,186 Joined: Dec 2013
RE: Perimeter of Ellipse
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).
03-24-2019, 12:42 PM
Post: #9
 Albert Chan Senior Member Posts: 1,658 Joined: Jul 2018
RE: Perimeter of Ellipse
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
07-11-2019, 05:02 PM
Post: #10
 ggauny@live.fr Senior Member Posts: 528 Joined: Nov 2014
RE: Perimeter of Ellipse
(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!

Gérard.
12-04-2019, 05:30 PM
Post: #11
 Gerson W. Barbosa Senior Member Posts: 1,419 Joined: Dec 2013
RE: Perimeter of Ellipse
(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.
12-04-2019, 06:28 PM
Post: #12
 Valentin Albillo Senior Member Posts: 776 Joined: Feb 2015
RE: Perimeter of Ellipse
(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.

All My Articles & other Materials here:  Valentin Albillo's HP Collection

12-04-2019, 10:27 PM
Post: #13
 Gerson W. Barbosa Senior Member Posts: 1,419 Joined: Dec 2013
RE: Perimeter of Ellipse
(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.
12-05-2019, 03:17 PM
Post: #14
 Gerson W. Barbosa Senior Member Posts: 1,419 Joined: Dec 2013
RE: Perimeter of Ellipse
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
12-07-2019, 07:57 PM (This post was last modified: 01-05-2020 01:03 AM by Gerson W. Barbosa.)
Post: #15
 Gerson W. Barbosa Senior Member Posts: 1,419 Joined: Dec 2013
RE: Perimeter of Ellipse
(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
12-29-2019, 07:58 PM
Post: #16
 Gerson W. Barbosa Senior Member Posts: 1,419 Joined: Dec 2013
RE: Perimeter of Ellipse
(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.
01-19-2020, 03:56 AM
Post: #17
 Albert Chan Senior Member Posts: 1,658 Joined: Jul 2018
RE: Perimeter of Ellipse
(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.
01-19-2020, 11:00 PM (This post was last modified: 01-20-2020 03:26 PM by Albert Chan.)
Post: #18
 Albert Chan Senior Member Posts: 1,658 Joined: Jul 2018
RE: Perimeter of Ellipse
(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.
01-21-2020, 02:37 AM (This post was last modified: 05-18-2021 07:07 PM by Gerson W. Barbosa.)
Post: #19
 Gerson W. Barbosa Senior Member Posts: 1,419 Joined: Dec 2013
RE: Perimeter of Ellipse
(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)]²
01-21-2020, 05:16 PM (This post was last modified: 01-23-2020 01:56 PM by Albert Chan.)
Post: #20
 Albert Chan Senior Member Posts: 1,658 Joined: Jul 2018
RE: Perimeter of Ellipse
(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
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)