Conference rusure::math

                     <<< Note 1502.0 by LIOS01::SAPIENZA >>>
                       -< Circumference of an ellipse? >-

>>       Given an ellipse described by the equation:
    
>>    		x^2	y^2
>>    		---  +  ---  = 1
>>    		a^2	b^2
    
>>       Does someone have a formula and/or program for determining the
>>    circumference of the ellipse?

From Handbook of Mathematics by I. N. Bronshtein and K. A. Semendyayev
(ISBN 0-442-21171-6).

Assume a > b (i.e., the ellipsis is "horizontal").

Then the distance between one focus and the origin is c = sqrt(a*a - b*b) and
the excentricity is e = c / a.

The circumference of the ellipse is

	L = 4 * pi * E(e)

			    1		 1*3		 1*3*5
	  = 2 * pi * (1 - (---)^2*e^2 - (---)^2*e^4/3 - (-----)^2*e^6/5 ...)
			    2		 2*4		 2*4*6

where E(e) = E(e, pi/2) is the complete elliptic integral of the second kind.

If we set x = (a-b)/(a+b), then

	L = pi*(a+b)*(1 + x^2/4 + x^4/64 + x^6/256 + 25*x^8/16384 ...)

The elliptic integrals are reasonably well described in Abramowitz and Stegun;
they also give numerical methods.

    Re .0, .1:
    
    Just to forestall any questions, I would like to point out that the
    circumference of an ellipse cannot generally be expressed with a closed
    form formula using the common operations (multiplication, sine,
    logarithm, et cetera).  The series in .1 is all the answer there is,
    except for changing the form in various ways.
    
    
    				-- edp

    
    Ok, while I'm trying to understand the series formula, perhaps you can
    explain why the circumference of a circle can degenerate to 2(pi)R. I'm
    assuming it follows somehow from the equation in .1s reply.
    
    
    Frank
    

RE: .3 (Frank)

    Here are the key parts of the equation in .1:

>	L = pi*(a+b)*(1 + x^2/4 + x^4/64 + x^6/256 + 25*x^8/16384 ...)

    where

>	x = (a-b)/(a+b)

    a=b=R for a circle, thus x = (R-R)/(R+R) = 0/2R = 0.  So L =
    pi*(R+R)*(1 + 0 + 0 + ...) = 2piR.

					Topher

re .3    
>    Ok, while I'm trying to understand the series formula, perhaps you can
>    explain why the circumference of a circle can degenerate to 2(pi)R. I'm
>    assuming it follows somehow from the equation in .1s reply.

For the circle, take a = b = r.

From .1, computing the eccentricity of a circle yields c =
sqrt(a*a - b*b) = sqrt(r*r - r*r) = 0, and e = c / a = 0 / r = 0
(of course).  From the series expansion E(0) = 1, and that yields
L = 2 * pi * r (that last factor, probably a factor of a as
opposed to b or r, was missing in .1).

Or using the second form, with again with a = b = r, you get
L = pi*(a+b)*(1 + x^2/4 + ...) where x = (a-b)/(a+b) = 0,
which also reduces to L = 2 pi r.

The area of an ellipse I believe is simply pi * a * b, but the
formula for arc length doesn't "stretch" so nicely.  Consider a
square of side r being stretched differently in the x and y
directions so that it now is a rectangle with height a and width
b.  The area becomes ab, but the diagonal doesn't scale smoothly;
from r sqrt(2) it changes to sqrt(a*a + b*b).  If you start with
a rectangle of 2r by r and scale as before (r -> a vertically and
r -> b horizontally), the area moves from 2r^2 to 2ab but the
diagonal goes from r sqrt(5) to sqrt(4*a*a + b*b).  Note that was
scaled by a different factor than for the square.  Well when you
stretch a circle to an ellipse, a tiny segment of the
circumference at one point gets stretched by a different factor
than a segment of the circumference elsewhere with a different
slope.  So the end result isn't easy to add up...the integral
doesn't have a closed form solution, so they just named it and
called that the answer.  (Of course, its usefulness in other
context is what probably led to its being given a name.)

Dan

    
    Re .4
    
       Yes, I just saw that and came back to clarify my posting. Seems you
    beat me to it. My original confusion was in looking at the first equation
    given. That is:
    
    	L = 4 * pi * E(e)
    
    	  = 2 * pi * (1 - ....)
    
      After I decided to skip that and looked at the rest of .1s reply I
    saw where 2(pi)r came from. But I still have a question about the last
    term in the example.
    
       The posting in .1 gives the following:
    
    	L = pi * (a+b) * (1 + x^2/4 + x^4/64 + x^6/256 + 25*x^8/16384 ...)
    
       The question is, where did -----------------------^ this come from?
    
       Could someone post this in unexpanded form so that I can write a
    subroutine to calculate it? With the mystery '25' in the term, I can't
    guess as to what follows in the remaining terms.
    
    
    Thanks again,
    Frank
    

                     <<< Note 1502.6 by LIOS01::SAPIENZA >>>

>>       The posting in .1 gives the following:
    
>>    	L = pi * (a+b) * (1 + x^2/4 + x^4/64 + x^6/256 + 25*x^8/16384 ...)
    
>>       The question is, where did -----------------------^ this come from?

I don't know where the 25 comes from.

I suspect that the book I quoted in .1 got the "L" formula from the "E" formula
by expressing x in terms of e and performing *a lot* of arithmetic.  Another
possible approach is to make a variable transformation in the original integral
for the arc length.

A general reflection: while a power series is easy to evaluate, the error in
the approximation (or, equivalently, the number of terms you must use to get
some given accuracy) increases drastically as the argument increases (in this
case, as the ellipsis gets narrower and narrower).  Also (not applicable in
this case), as the argument gets big (say 4 or 5 or so), one can easily get
into problems of floating-point overflow.  If one needs to evaluate some
special function, he normally does best in calculating a minimax approximation
(unfortunately, I can't provide a reference for how to do this; maybe somebody
else can).

If you need a function to calculate E(e), I think you'd be better off by using
the following algorithm (taken from An Atlas Of Functions, by Spanier and
Oldham (ISBN 0-89116-573-8), section 61.8, pp 613-4) that is claimed to give 24
bits of precision.

WARNING: Their definition of E is off by a constant factor compared to the one
I gave in .1.  Also, they use "p" where I've used "e"

	set g = sqrt(1-p*p)	! p is the argument to the function
	set a = t = e = 1
	set e = e + g * g
    2$:	set s = g * a
	set a = (a + g) / 2
	set g = sqrt(s)
	set t = t + t
	set e = e - t * (a*a - s)
	if t*(a*a-s) >= 1E-8 then go to 2$ endif
	set k = pi / (g + g)
	set e = e * k / 2
	return e (i.e., E(p))

They suggest using p = .5 as test data; it should return 1.46746221

   The general term in the series in x is

	1.3.5...(2n-3)
       (--------------)^2
           (2^n).n!

    For example (omitting the squares)

     1    1     5      7      7.9      9.11      9.11.13    11.13.15
    ---  ---  ------ ------ -------- --------- ----------- ----------- ...
    4.2  8.2  16.2.4 32.2.4 64.2.4.6 128.2.4.6 256.2.4.6.8 512.2.4.6.8

    The algorithm earlier is also given in AMS55 - it is basically the
    quadratically convergent arithmetic geometric mean iteration.

    - Jim

I have translated the algorithm in .6 to FORTRAN, if anyone wants it. Just 
ask.

>They suggest using p = .5 as test data; it should return 1.46746221

I get 1.4674618244 using either REAL*4 or REAL*8 data. I suspect the author 
actually used 1.E-6 for a cutoff.