Elliptic integrals of 1st and 2nd kind calculated by complex agm Gjermund Skailand Junior Member Posts: 15 Threads: 2 Joined: Oct 2008 06-27-2013, 06:00 PM 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) Member Posts: 167 Threads: 33 Joined: Jul 2011 06-27-2013, 09:26 PM A wonderful combination of two amazing mathematical insights! Thanks. Paul Dale Posting Freak Posts: 3,229 Threads: 42 Joined: Jul 2006 06-27-2013, 11:10 PM This should be fairly straightforward to implement on the 34S. It has the complex AGM function built in :) - Pauli Gjermund Skailand Junior Member Posts: 15 Threads: 2 Joined: Oct 2008 06-29-2013, 03:39 PM ```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 ``` « Next Oldest | Next Newest »

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

Forum Jump: