HP Forums

(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".

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)

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)*ε


>>> 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)

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.


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

(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!

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

(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" ?

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

(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.

(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.

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

(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

(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

(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

+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.


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


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