Post Reply 
HP42s first major program (Double Integral) Best way to approach?
06-08-2020, 03:35 PM
Post: #56
RE: HP42s first major program (Double Integral) Best way to approach?
Here is EKmc() non-recursive Python code.
Note: code preceded by "from math import sqrt, pi", for complex number support, use cmath.

Code:
def EKmc(m, verbal=False, c=pi/2):
    H = []          # non-recursive version, by collecting h's
    while abs(1+m)-1: h = m/(1+sqrt(1-m))**2; H.append(h); m = h*h
    k = 0.25*c*m
    e = -k
    for h in reversed(H):
        m = h*h
        if verbal: print 'm =',m,'\n\tE(m)-c =',e,'\n\tK(m)-c =',k
        e = (e+e-k + m*k + (m-h)*c) / (1+h)        
        k = k + h*k + h*c
    return e, k

Try HV(1,1000) again. With d=1, we can simplify the HV formula.

>>> D = 1000
>>> e, k = EKmc(1/(D*D), verbal=True)
m = 2.44141113282e-28
      E(m)-c = -9.58739909907e-29
      K(m)-c = 9.58739909907e-29
m = 6.25000625001e-14
      E(m)-c = -2.45437171499e-14
      K(m)-c = 2.45437171499e-14
>>> e, k      # m = 1e-6
(-3.9269915532983253e-07, 3.9269930259211089e-07)

>>> D/3 * ((e+k) + D*D*(e-k) + pi)       # = HV(1,D)
785.39806522266554

EKmc has the speed of AGM2 method, but avoided catastrophic cancellation issue.
It also avoided AGM2 (complex number) convergence issue.

With K(m), we can recover AGM's: AGM(1, b) * K(1-b*b) = pi/2
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
RE: HP42s first major program (Double Integral) Best way to approach? - Albert Chan - 06-08-2020 03:35 PM



User(s) browsing this thread: 20 Guest(s)