Elliptic integrals of 1st and 2nd kind calculated by complex agm

[Gjermund Skailand]

Hi.
I do not know is this known, there was a recent thread about calculating perimeter of ellipses and elliptic integrals of first and second kind. There was presented an iteration with modified agm formula, but no reasoning why it worked.

In general refer to http://en.wikipedia.org/wiki/Elliptic_integral

Which states that elliptic integrals of 1st kind K(k) can most efficiently be computed in terms of the arithmetic-geometric mean:

K(k)=(pi/2)/agm(1-k, 1+k) (1)

For elliptic integrals of 1st and 2nd kind we have following relationship

dK(k) /dk = E(k)/(k(1-k^2)) – K(k)/k (2)

which can be re-arranged into

E(k)=k(1-k*k) dK(k)/dk + (1-k*k) K(k) (3)

Now take a second look on the agm,

K(zk)=(pi/2)/agm(1-zk, 1+zk) (4)

and allow zk be complex and define such iteration to define complex K.

This is useful because if we remember what Valentine taught us about numeric derivation a while ago: numeric derivation rule, let dk be in complex plane instead of real , we then have two advantages, we can “freely” choose the step dk without numerical problems and compute derivate for real values of k using s side-step in imaginary plane:

dK(k)/dk = im(K(k,dk))/ dk (5)
thus E(k) can be computed directly by computing K(zk) and dK(zk) with complex agm .

Further, provided dk is small enough, to avoid “interference” between real and imaginary values by choosing a suitable value for dk, say for the calculator use dk= (0, 1E-12) then
you can compute both dK(k)/dk and K(k) in same complex agm iteration.

thus if we want to calculate E(k) we use initial values

k0=k+i1E-12
A0=1-k0
B0=1+k0
A0 = (1-k, -1E-12)
B0 = (1+k, 1E-12)
Do complex AGM
{ Ai= (A0+B0)/2; Bi=SQRT(A*B) }
until Ai==Bi then
Zn=0.5*pi/An
K(k)= REAL(Zn)
dK(k)/dk= IM(Zn)/1E-12
E(k)= k(1-k^2) dK(k)/dk + (1-k^2) K(k)

In user rpl
AGMc << 1. 12. START DUP 2 + 2 / UNROT * SQRT NEXT DROP >>

Kdk = << 1E-12 DUP2 R->C -> k dk z << 1 z – 1 z + AGMc pi 2. / SWAP / C->R dk / >> >>

Ek = << -> k << k Kdk 1. k k * - DUP UNROT k * * UNROT * + >> >>

Notice the almost one-liner programs!

This can be made even more efficient considering that we in principle could use extended (infinite) precision with small dk. Since what we actually want to do is keep the real and imaginary values separately and we are free to choose dk as we like so we can split and simplify (de-complexify) as shown in part 2.

Any comments are welcome
Best regards
Gjermund

[Peter Murphy (Livermore)]

A wonderful combination of two amazing mathematical insights! Thanks.

[Paul Dale]

This should be fairly straightforward to implement on the 34S. It has the complex AGM function built in :)

- Pauli

[Gjermund Skailand]

Part 2 Making the calculations more efficient

The complex AGM can be simplied.

1 In the numerical derivation above we used complex numbers, 

chose dt = 1E-12 and abused floating point calulations because 

we wanted to keep real and imaginary parts separated. 

We were then using properties of complex algebra and floating point numbers, 

but we can do better with some manual work.

2) Consider what happens with a complex numbers where the real values are similar size, 

and the imaginary values are much smaller than the real values.

Then any multiplication between real and an imaginary value will be

of same small size as imaginary number. 

Any multiplication with two small imaginary numbers will be very

 very small and be of second order and disregarded.  

3) Therefore split all complex operations into real operations and simplify according to the above.

4) Since formulas are separated we now use idk =1 instead of 1E-12

We end up with:  

E(k): 

k0=k+i

A0=1-k0

B0=1+k0 

ar0=1-k

ai0=-1

br0=1+k

bri=1

Do modified "complex" AGM 

{ 

ar1 = 0.5(ar0 + br0)

ai1 = 0.5(ai0 + bi0)

br1 = sqrt(ar0*br0)  

bi1 =( 0.5(ar0*bi0 + ai0*br0))/ br1

}

iterate until (arn~brn ,  ain~bin)

then

K=K(k)=0.5pi/sqrt(arn)

dK=dK(k)/dk=  -0.5pi*ain/arn^2  

E(k)= k(1-k*k) dK + (1-k*k) K

This is similar to S. Adlaj's  modified agm iteration formula.

In user rpl

Kdk2 

<< 1. OVER -  -1.  1. 4 PICK + 1.  

  -> k ar ai br bi  

  << 1. 12. 

    START 

      .5 ar br + * 

      .5 ai bi + * 

      ar br * sqrt 

      .5 ar bi * ai br * + * OVER / 

      'bi' STO 'br' STO 'ai' STO 'ar' STO 

    NEXT 

    pi 2 / ar  / 

    DUP ai ar / NEG * 

  >>  

>>

Thus with this algorithm dK/dk is proportional to K * AI/AR which I think is a nice feature.

 Since all starting values are of similar size then the ar and ai will converge at similar speeds.

To show that this actually works we can compute pi for any k between 0 and 1 by Legendres equation

K(k) E(c) + E(k)K(c) - K(k) K(c) = pi/2,       

c=sqrt(1-k^2)

KE 

<< -> k  << k Kdk2 OVER UNROT 1. k k * - DUP UNROT k * * UNROT * + >> >>

Test

<< 

1 OVER DUP * - sqrt -> k c << k KE c KE >>  

-> Kk Ek Kc Ec << Kk Ec * Ek Kc * + Kk Kc * - 2 * >> 

>>

It is easy to check the accuracy since any value between 0 and 1 should result in pi .

If there is any interest I can post source code for  hp48 and hp50g of K and E in sys-rpl in extended precision 

Any comments are welcome

Best regards

Gjermund