Perimeter of Ellipse
01-23-2020, 01:40 PM (This post was last modified: 01-24-2020 05:04 PM by Albert Chan.)
Post: #21
 Albert Chan Senior Member Posts: 972 Joined: Jul 2018
RE: Perimeter of Ellipse
(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".
06-05-2020, 03:28 AM (This post was last modified: 06-06-2020 11:33 AM by Albert Chan.)
Post: #22
 Albert Chan Senior Member Posts: 972 Joined: Jul 2018
RE: Perimeter of Ellipse
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)
06-05-2020, 08:09 PM (This post was last modified: 06-06-2020 05:14 PM by Albert Chan.)
Post: #23
 Albert Chan Senior Member Posts: 972 Joined: Jul 2018
RE: Perimeter of Ellipse
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)
06-06-2020, 05:12 PM (This post was last modified: 06-09-2020 12:55 AM by Albert Chan.)
Post: #24
 Albert Chan Senior Member Posts: 972 Joined: Jul 2018
RE: Perimeter of Ellipse
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
08-01-2020, 12:31 PM
Post: #25
 Albert Chan Senior Member Posts: 972 Joined: Jul 2018
RE: Perimeter of Ellipse
(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 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
 « Next Oldest | Next Newest »

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