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



#2

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


#3

A wonderful combination of two amazing mathematical insights! Thanks.

#4

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


- Pauli


#5

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


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP Prime: complex numbers in CAS. Alberto Candel 1 243 12-06-2013, 02:36 PM
Last Post: parisse
  [HP Prime] Plots containing complex numbers bug? Chris Pem10 7 465 12-05-2013, 07:40 AM
Last Post: cyrille de Brébisson
  Complex Number Entry on Prime Jeff O. 19 894 11-16-2013, 12:34 PM
Last Post: Jeff O.
  My 2nd and 3rd Prime Programs Jeff O. 7 420 10-31-2013, 08:18 AM
Last Post: Jeff O.
  HP Prime complex results Javier Goizueta 0 141 10-06-2013, 12:59 PM
Last Post: Javier Goizueta
  HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 659 09-30-2013, 03:44 AM
Last Post: From Hong Kong
  [HP-Prime xcas] operations with complex numbers + BUGs + Request CompSystems 9 552 09-08-2013, 10:40 PM
Last Post: CompSystems
  Virtual HP-IL Video Interface ILVideo the 2nd! Christoph Giesselink 3 262 08-15-2013, 06:49 PM
Last Post: Sylvain Cote
  HP-11C and complex numbers Antlab 5 412 06-28-2013, 08:59 AM
Last Post: Antlab
  71B Complex User Defined Functions in Calc Mode? Michael Burr 0 147 03-18-2013, 09:38 PM
Last Post: Michael Burr

Forum Jump: