The following warnings occurred:
Warning [2] Undefined array key 245627 - Line: 275 - File: inc/plugins/threaded_mode.php PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php 275 errorHandler->error_callback
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Undefined array key 245628 - Line: 275 - File: inc/plugins/threaded_mode.php PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php 275 errorHandler->error_callback
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Undefined array key 245630 - Line: 275 - File: inc/plugins/threaded_mode.php PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php 275 errorHandler->error_callback
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Undefined array key 245736 - Line: 275 - File: inc/plugins/threaded_mode.php PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php 275 errorHandler->error_callback
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Undefined variable $thread - Line: 295 - File: inc/plugins/threaded_mode.php PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php 295 errorHandler->error_callback
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Trying to access array offset on value of type null - Line: 295 - File: inc/plugins/threaded_mode.php PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php 295 errorHandler->error_callback
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Undefined variable $fid - Line: 295 - File: inc/plugins/threaded_mode.php PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php 295 errorHandler->error_callback
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Undefined array key 245628 - Line: 331 - File: inc/plugins/threaded_mode.php PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php 331 errorHandler->error_callback
/inc/plugins/threaded_mode.php 332 ThreadedMode::buildtree
/inc/plugins/threaded_mode.php 304 ThreadedMode::buildtree
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Undefined array key 245736 - Line: 331 - File: inc/plugins/threaded_mode.php PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php 331 errorHandler->error_callback
/inc/plugins/threaded_mode.php 332 ThreadedMode::buildtree
/inc/plugins/threaded_mode.php 332 ThreadedMode::buildtree
/inc/plugins/threaded_mode.php 304 ThreadedMode::buildtree
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Undefined variable $theme - Line: 3 - File: inc/plugins/threaded_mode.php(305) : eval()'d code PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php(305) : eval()'d code 3 errorHandler->error_callback
/inc/plugins/threaded_mode.php 305 eval
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Trying to access array offset on value of type null - Line: 3 - File: inc/plugins/threaded_mode.php(305) : eval()'d code PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php(305) : eval()'d code 3 errorHandler->error_callback
/inc/plugins/threaded_mode.php 305 eval
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Undefined variable $theme - Line: 3 - File: inc/plugins/threaded_mode.php(305) : eval()'d code PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php(305) : eval()'d code 3 errorHandler->error_callback
/inc/plugins/threaded_mode.php 305 eval
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Trying to access array offset on value of type null - Line: 3 - File: inc/plugins/threaded_mode.php(305) : eval()'d code PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php(305) : eval()'d code 3 errorHandler->error_callback
/inc/plugins/threaded_mode.php 305 eval
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Undefined variable $lang - Line: 5 - File: inc/plugins/threaded_mode.php(305) : eval()'d code PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php(305) : eval()'d code 5 errorHandler->error_callback
/inc/plugins/threaded_mode.php 305 eval
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks
Warning [2] Attempt to read property "messages_in_thread" on null - Line: 5 - File: inc/plugins/threaded_mode.php(305) : eval()'d code PHP 8.1.2-1ubuntu2.14 (Linux)
File Line Function
/inc/class_error.php 153 errorHandler->error
/inc/plugins/threaded_mode.php(305) : eval()'d code 5 errorHandler->error_callback
/inc/plugins/threaded_mode.php 305 eval
/inc/plugins/threaded_mode.php 23 ThreadedMode::showthread_threaded
/inc/class_plugins.php 142 threaded_mode_showthread_threaded
/showthread.php 918 pluginSystem->run_hooks





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



#5

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


#6

A wonderful combination of two amazing mathematical insights! Thanks.

#7

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


- Pauli


#8

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 1,856 12-06-2013, 02:36 PM
Last Post: parisse
  [HP Prime] Plots containing complex numbers bug? Chris Pem10 7 3,551 12-05-2013, 07:40 AM
Last Post: cyrille de Brébisson
  Complex Number Entry on Prime Jeff O. 19 5,024 11-16-2013, 12:34 PM
Last Post: Jeff O.
  My 2nd and 3rd Prime Programs Jeff O. 7 2,451 10-31-2013, 08:18 AM
Last Post: Jeff O.
  HP Prime complex results Javier Goizueta 0 956 10-06-2013, 12:59 PM
Last Post: Javier Goizueta
  HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 4,078 09-30-2013, 03:44 AM
Last Post: From Hong Kong
  [HP-Prime xcas] operations with complex numbers + BUGs + Request CompSystems 9 3,361 09-08-2013, 10:40 PM
Last Post: CompSystems
  Virtual HP-IL Video Interface ILVideo the 2nd! Christoph Giesselink 3 1,542 08-15-2013, 06:49 PM
Last Post: Sylvain Cote
  HP-11C and complex numbers Antlab 5 2,148 06-28-2013, 08:59 AM
Last Post: Antlab
  71B Complex User Defined Functions in Calc Mode? Michael Burr 0 905 03-18-2013, 09:38 PM
Last Post: Michael Burr

Forum Jump: