Perimeter of Ellipse

[Albert Chan]

(01-21-2020 05:16 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.

We could even do it in a novel way, from the *more* eccentric ellipse.
Below, LHS is the more eccentric ellipse, with \({a_0 + b_0 \over 2} = a, \sqrt{a_0 b_0} = b\)

\( p(a_0,b_0) = 2× \left( p(a,b) - {\pi b^2 \over AGM(a,b)}\right)\)

Using my Casio FX-115MS, AGM2 method for Halley's comet, numbers from here

A = 2667950000
B = 678281900
C = A²
D = B²
E = 0.5

F=B-A : C=C-EF² : B=√AB : A=A+F/2 : E=2E

Quote:\(p(a,b) = {p(a_0,b_0) \over 2} + {\pi b^2 \over AGM(a,b)}\)

Keep pressing "=" until F=0 (so that A,B,C all converged)

\(\large\pi\)(C+D)/B       → ellipse perimeter = 11,464,318,984.1 km

This was similar to my lua code, with C=(A²+B²)/2, E=0.25, so that converged C/B = "radius".

[Albert Chan]

From previous 2 posts, we can get some nice Elliptic Integral Identity:
Below, all Elliptic functions uses modulus argument.

Eccentricity of ellipse, e^2 = 1 - (b/a)^2
We proved "next" ellipse (a' = (a+b)/2 , b' = sqrt(a*b)), has eccentrity of h = (a-b)/(a+b)

K(e) = pi/2 / agm(1,b/a) = pi/2 / agm(1+e,1-e)
K(h) = pi/2 / agm(1+h,1-h) = pi/2 / agm(1,b/a) / (1+h)

→ K(e) = (1+h) K(h)

(01-21-2020 05:16 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)\)

With ellipse_perimeter = 4 a E(e), we have:

4 a E(e) = 4 (a+b) E(h) - 4 b K(e) = 4 (a+b) E(h) - 4 b (1+h) K(h)

→ E(e) = 2/(h+1) * E(h) - (1-h) * K(h)

[Albert Chan]

We can get E(e),K(e) by recursively going for E(h),K(h), like this.
Note: code use parameter m, not modulus k, m = k^2
Note: code quit when m is small enough for taylor linear term estimate.

E(m=ε) ≈ pi/2 - (pi/8)*ε ,     K(m=ε) ≈ pi/2 + (pi/8)*ε

Code:

from cmath import sqrt, pi
    
def EK(m, verbal=False, c=pi/2):
    e = abs(m)
    if 1+e*e == 1: m *= 0.25*c; return c-m, c+m
    b = sqrt(1-m)
    if not b: return 1, 99  # K(1) assumed 99 (should be inf)
    h = m / ((b+2)*b+1)     # = (1-b)/(1+b)
    e, k = EK(h*h, verbal)
    hk = h*k
    e = e-k + hk + b*e      # = 2/(1+h)*e - (1-h)*k
    k = hk + k              # = (1+h)*k
    if verbal: print 'm = %s\n\tE(m) = %s\n\tK(m) = %s' % (m, e, k)
    return e, k

>>> e, k = EK(0.96, verbal=True)
m = (2.89332317725e-05+0j)
      E(m) = (1.57078496468+0j)
      K(m) = (1.57080768903+0j)
m = (0.0212862362522+0j)
      E(m) = (1.56240357945+0j)
      K(m) = (1.57925700384+0j)
m = (0.444444444444+0j)
      E(m) = (1.3781039379+0j)
      K(m) = (1.80966749549+0j)
m = 0.96
      E(m) = (1.05050222698+0j)
      K(m) = (3.01611249248+0j)
>>> 4 * 50 * abs(e)                       # ellipse_perimeter(10,50), error = -3 ULP
210.10044539689011

>>> e, k = EK(2, verbal=True)        # see https://www.hpmuseum.org/forum/thread-15...#pid132745
m = (5.57959210499e-05-0j)
      E(m) = (1.57077441556+0j)
      K(m) = (1.57081823849+0j)
m = (0.0294372515229-0j)
      E(m) = (1.55917174457+0j)
      K(m) = (1.58255172722+0j)
m = (-1-0j)
      E(m) = (1.91009889451+0j)
      K(m) = (1.31102877715+0j)
m = 2
      E(m) = (0.599070117368+0.599070117368j)
      K(m) = (1.31102877715-1.31102877715j)

[Albert Chan]

Improvement to my previous EK(m)
When m is small, the term E(m) - K(m) is losing significant digits.
Worse, catastrophic cancellation occurs at the base of recursion.

EKmc(m) returns E(m)-c, K(m)-c, where c=pi/2, similarly to expm1(x) = exp(x) - 1
EK(m) is now a simple wrapper for EKmc(m), adding back the c's.

Code:

from cmath import sqrt, pi
c = pi/2

def EKmc(m):
    'return E(m)-c, K(m)-c, where c=pi/2, m!=1'
    if 1+abs(m) == 1: m *= 0.25*c; return -m, m
    b = sqrt(1-m)
    h = m / ((b+2)*b+1)     # = (1-b)/(1+b)
    e, k = EKmc(h*h)
    hk, hc = h*k, h*c
    return (hk-k) + (e+b*e) - b*hc, (hk+k) + hc

def EK(m):
    'return E(m), K(m), assumed K(1)=99 (should be inf)'
    if abs(m-1)==0: return 1, 99
    e, k = EKmc(m)
    return e+c, k+c

>>> e, k = EK(0.96)
>>> 4 * 50 * abs(e)                       # ellipse_perimeter(10,50), error = 0 ULP
210.10044539689002

[Albert Chan]

(01-21-2020 05:16 PM)Albert Chan Wrote:  

(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

If we are at focus F2, and we let distance a0 = AF2, b0 = BF2, eccentricity of ellipse is also h0

ellipse major-axis a = OB = (a0+b0)/2
ellipse minor-axis b = OC = √((CF2)² - (OF2)²) = √(a² - (a-b0)²) = √((2a-b0)*b0) = √(a0*b0)

This matched ellipse perimeter recursive AGM formula!

Example, lets try the earth-moon orbit, average orbital distance 384748 km

a0 = 406731 km
b0 = 364397 km

e = h0 = (a0-b0)/(a0+b0) = 42334/771128 ≈ 0.0549

Update, Jan 26, 2023
In 2 years, California State University, San Bernardino picture link is dead.
It is ironic that CSUSB motto is "We Define The Future".
Above picture were from WayBack Machine, Thanks!

[hazem]

My name is Hazem Albadry from Iraq born in 1956

Finally I discovered the exact formula of Perimeter of an Ellipse name (Hazem C-Ellipse) as following:

C = a * ((π-2)/45 * arcsine(b/a) + 4)

Where a is the major radius of ellipse
Where b is the minor radius of ellipse

When a=b then the formula will be C = 2aπ (As circle)
And when b=0 then the formula will be C = 4a (As lines)

I hope this help anyone who interested in Perimeter of an Ellipse.

Best regards.

Hazem Hameed Rashid Albadry
Hazemhameed90@gmail.com
Baghdad / Iraq

[rprosperi]

(04-11-2023 09:43 PM)hazem Wrote:  My name is Hazem Albadry from Iraq born in 1956

Finally I discovered the exact formula of Perimeter of an Ellipse name (Hazem C-Ellipse) as following:

Thank you for sharing this here.

In the equations, what is "n" ?

[Werner]

Look closely, Bob, that's not an 'n' ;-)
And I'm afraid to say that any mathematical formula that requires degrees mode for the arcsin is suspicious, to say the least.
Let's try it with A=10 and B=50 (it requires B>=A), Albert's example from above.
Then the result is 214.6..
So not even a very good approximation..

Cheers, Werner

[rprosperi]

(04-12-2023 05:43 AM)Werner Wrote:  Look closely, Bob, that's not an 'n' ;-)

Good eyes Werner, thanks. The ironic part is I paused and wondered why there was no Pi in the formula, which led me to notice the "n", which isn't an "n". You know the rest.

[Valentin Albillo]

(04-12-2023 05:43 AM)Werner Wrote:  Look closely, Bob, [...]

Correct. The proper way to unambiguously show Pi would be like this: \(\pi\)

Quote:And I'm afraid to say that any mathematical formula that requires degrees mode for the arcsin is suspicious, to say the least. [...] So not even a very good approximation..

Indeed. The OP says his formula is exact but nothing could be further from the truth. Computing the perimeter of a generic ellipse exactly requires, well, elliptic functions, not elementary functions such as inverse trigonometric ones.

V.

[floppy]

Exact formula for ellipse perimeter calculation? "Yes".

Via division of 2 infinite iterated functions MAGM and AGM (I dont call them elliptic functions; other does).
http://www.ams.org/notices/201208/rtx120801094p.pdf

On your calc, you just stop the AGM and MAGM routines when the delta of the iteration result is good enough for you (I stop when it is E-8; one time I tried E-9 but it went undefinitively and suppose this is due to the HP41 calculation approximation making a convergence unstable. A 128bit PC would give an outstanding precision).

Code for HP41 standard below (no SandMath ROM or other modules required) which can be overtaken for other computer.
The documentation should make possible a full understanding of the "iterated" behaviour.

Since everything can be improved.. any suggestion to improve the code is welcome.

Use case:
10
ENTER
6
XEQ ALPHA ELPER ALPHA
51.05399775

Update1:
- code revisited/simplified since the convergence is stable on HP41 and an exit criteria like E-8 not mandatory
- file upload redone

Code:

; Ellipse perimeter calculation according MAGM and AGM
;
; under CC BY SA CreativeCommons 4.0 floppy @ https://www.hpmuseum.org
;
; change log
;  2023 04 16 creation
;
; verification https://www.mathsisfun.com/geometry/ellipse-perimeter.html
; 
; Register overwritten
;   00  a then later a**2..
;   01  b then later b**2..
;   02  zn   parameter
;   03  AGM  temporary result
;   
LBL "ELPER" 
;                    X   Y   Z   T
STO 00            ;  a   b
                  ;  ... a in R00  
X=0?
GOTO 00
X<>Y
X=0?
GOTO 00
STO 01            ; ... b in R01
CF 00             ;  clear the flag 00 used in AGM
XEQ 01            ; "AGM"
STO 03            ; ... agm in R03
0
STO 02
RCL 00
X^2 
STO 00            ; a^2
                  ; R00 content modified
RCL 01
X^2
STO 01            ; b^2
                  ; R01 content modified
CF 00
XEQ 02            ; "MAGM"
RCL 03
/
2
*
PI
*
RTN
;
; MAGM calculation with recursive function
;   http://www.ams.org/notices/201208/rtx120801094p.pdf
;
; under CC BY SA CreativeCommons 4.0 floppy @ https://www.hpmuseum.org
;
LBL 02  
;                   X            Y           Z              T
;                   an           bn  
+
2
/                ; (an+bn)/2
RCL 00
RCL 01
RCL 02           ; zn            bn          an          (an+bn)/2
ST- Y
ST- Z           ;  zn           an-zn      bn-zn         (an+bn)/2
X<> Z           ;  bn-zn        an-zn        zn          (an+bn)/2
*               ;  (an-zn)*(bn-zn)    zn    (an+bn)/2    (an+bn)/2
SQRT            ;  SQRT((an-zn)*(bn-zn))    zn    (an+bn)/2    (an+bn)/2
+               ;  (SQRT(an-zn)*(bn-zn))+zn  (an+bn)/2     (an+bn)/2   (an+bn)/2
STO 01          ;  bN         (an+bn)/2     (an+bn)/2   (an+bn)/2
                ;  ... bN in R01
LASTX           ;  (SQRT(an-zn)*(bn-zn))   bN   (an+bn)/2   (an+bn)/2
RCL 02          ;  zn    (SQRT(an-zn)*(bn-zn))   bN      (a+b)/2
-
CHS             ;  zn-(SQRT(an-zn)*(bn-zn))  bN  (an+bn)/2   (an+bn)/2
STO 02          ;  zN           bN          (an+bn)/2      (an+bn)/2
RDN             ;  bN          (an+bn)/2      (an+bn)/2     zN
ENTER           ;  bN           bN          (an+bn)/2   (an+bn)/2
X<> Z           ;  (an+bn)/2    bN          bN          (an+bn)/2
STO 00          ;  aN           bN          bN          aN
X=Y?
SF 00           ; convergence according cal precision achieved
FC?C 00
GOTO 02
RTN
;
; AGM calculation with recursive function
;   https://de.wikipedia.org/wiki/Arithmetisch-geometrisches_Mittel
;
LBL 01           ;  X   Y    Z    T
;                   a   b
STO Z            ;  a   b    a
X<>Y             ;  b   a    a
ST+ Z            ;  b   a    a+b
STO T            ;  b   a    a+b  b
*                ;  ba  a+b  b    b
SQRT             ;  sqrt(ba) a+b  b   b
R^               ;  b  sqrt(ba) a+b  b 
X=Y?
SF 00            ; convergence according cal precision achieved
RDN              ;  sqrt(ba) a+b  b   b
X<>Y
2
/                ; (a+b)/2  SQRT(ab)  b   b
;                    =a        =b   below..
FC?C 00
GOTO 01
RTN
;
LBL 00
+
4
*
RTN
END

[hazem]

(04-12-2023 01:53 AM)rprosperi Wrote:  

(04-11-2023 09:43 PM)hazem Wrote:  My name is Hazem Albadry from Iraq born in 1956

Finally I discovered the exact formula of Perimeter of an Ellipse name (Hazem C-Ellipse) as following:

Thank you for sharing this here.

In the equations, what is "n" ?

It is Pi

[floppy]

(04-13-2023 02:06 PM)hazem Wrote:  

(04-12-2023 01:53 AM)rprosperi Wrote:  Thank you for sharing this here.

In the equations, what is "n" ?

It is Pi

Can you give the result with a=10 and b=6? (I dont get it)
Then I will compare with approximations available here https://www.mathsisfun.com/geometry/elli...meter.html

[Albert Chan]

(04-12-2023 07:22 PM)floppy Wrote:  On your calc, you just stop the AGM and MAGM routines when the delta of the iteration result is good enough for you
(I stop when it is E-8; one time I tried E-9 but it went [i]ndefinitively ...

If you compare sucessive GM's or AM's (but not AM vs GM), it will not get into infinite loop.
see https://www.hpmuseum.org/forum/thread-15...#pid134773

FNP(A,B) for perimeter of ellipse, by MAGM, until full convergence of GM's

(06-04-2021 07:14 PM)Albert Chan Wrote:  10 DEF FNP(A,B) @ S=A*A+B*B @ T=1
20 REPEAT @ G=B @ K=(A-B)/2 @ B=SQRT(A*B) @ A=A-K @ T=T+T @ S=S-T*K*K @ UNTIL G=B
30 FNP = S/B*PI @ END DEF

>fnp(10, 6)
51.0539977270

Exact ellipse perimeter, rounded to 12 digits = 51.0539977268

[floppy]

+1.
Now leaving the "simplicity" of having MAGM divided by AGM, at the cost of calculation time, for calculating ellipse perimeter, and lets show what I call "GaussKummer-like" with "AGM functional turbo boosting". If I remember it was already posted in this forum for HP15c. Now for HP41. NO APPROXIMATION: only limitation to the available precision with HP41.

Code:

; calcul ellipse perimeter GaussKummer like with AGM exact function for 
;   "convergence boosting"
;
;Inputs
;  a ENTER b ENTER XEQ PERE12  
;  a and b: half parameter of the ellipse
;
;Outputs
;  results in stack X:perimeter
;
;Modules used
; None
;
;use R 07..11
;keep free Reg 00-06 if use with SOL from MATH from other programs
; in case this program would be used for example for finding the "a" if
; the "perimeter" and the "b" are given
;
;create raw file with "rpncomp --raw-output PERE12.TXT", upload
;  in V41 emulator and store into virtual drive pyILPER
;
;under CC BY SA CreativeCommons 4.0 floppy @ https://www.hpmuseum.org/forum/
;
;idea from here
;        https://www.hpmuseum.org/forum/thread-5820-post-171326.html#pid171326
;
;change log
; date 2023 04 15 creation and release (based on PEREL6 simplified program)
;
;comment
; all an bn R0xn are from the "n" loop
; all with "N" are updated parameter in the new loop
;
LBL "PERE12"              
X=0?         ; in case a = 0 this is a flat ellipse
GTO 01
STO 07       ; b in 07
X<>Y
X=0?         ; in case b = 0 this is a flat ellipse
GTO 01
STO 08       ; a in 08
X^2
X<>Y
X^2
+
STO 09       ; an^2 + bn^2 in R09 
1
STO 10       ; 1 in R10n
;
LBL 00       ;     X             Y          Z          T
  RCL 08     ;     an
  RCL X      ;     an           an
  RCL 07     ;     bn           an          an
  *          ;     bn*an          an
  SQRT       ;  SQRT(an*bn)=bN    an              
  X<> 07     ;     bn             an                
             ;  ... R07 = bN
  RCL 07     ;     bN             bn        an
  X=Y?
  SF 00
  RDN        ;     bn             an       
  -          ;    an-bn           
  2          ;    2            an-bn         
  ST* 10     ;    2            an-bn              
             ;    ... R10N= 2*R10n
  /          ;  (an-bn)/2          
  ST- 08     ;  (an-bn)/2                       
             ;  ... R08N=aN=an-(an-bn)/2=(an+bn)/2
  X^2        ;  ((an-bn)/2)^2     
  RCL 10     ;  R10N      ((an-bn)/2)^2      
  *          ;  R10N*((an-bn)/2)^2           
  ST- 09     ;  R10N*((an-bn)/2)^2             
             ;  ...  R09N = R09n- R10N*((an-bn)/2)^2
  FC?C 00
  GOTO 00
RCL 09
RCL 07 
/
PI
*          
RTN          
LBL 01       ;in case b or a = 0 flat ellipse it goes here out
+
4
*
RTN      
END

Just tested on an HP41CX. Printing w/o comment (in order to see the length).

Code:

 01*LBL "PERE12"
 02 X=0?
 03 GTO 01
 04 STO 07
 05 X<>Y
 06 X=0?
 07 GTO 01
 08 STO 08
 09 X^2
 10 X<>Y
 11 X^2
 12 +
 13 STO 09
 14 1
 15 STO 10
 16*LBL 00
 17 RCL 08
 18 RCL X
 19 RCL 07
 20 STO T
 21 *
 22 SQRT
 23 STO 07
 24 RDN
 25 -
 26 CHS
 27 2
 28 ST* 10
 29 /
 30 ST- 08
 31 X^2
 32 RCL 10
 33 *
 34 ST- 09
 35 RDN
 36 X#Y?
 37 GTO 00
 38 1/X
 39 RCL 09
 40 *
 41 PI
 42 *
 43 RTN
 44*LBL 01
 45 +
 46 4
 47 *
 48 RTN
 49 .END.

Now latest shortest version on HP41 (no check of flat ellipse; no modules required)

Code:

 01*LBL "PERE12"
 02 STO 07
 03 X^2
 04 X<>Y
 05 STO 08
 06 X^2
 07 +
 08 STO 09
 09 1
 10 STO 10
 11*LBL 00
 12 RCL 08
 13 RCL X
 14 RCL 07
 15 STO T
 16 *
 17 SQRT
 18 STO 07
 19 RDN
 20 -
 21 CHS
 22 2
 23 ST* 10
 24 /
 25 ST- 08
 26 X^2
 27 RCL 10
 28 *
 29 ST- 09
 30 RDN
 31 X#Y?
 32 GTO 00
 33 1/X
 34 RCL 09
 35 *
 36 PI
 37 *
 38 RTN
 39 END