Post Reply 
(34S) Integration (using the double-exponential method)
03-24-2017, 02:21 PM (This post was last modified: 06-15-2017 01:21 PM by Gene.)
Post: #1
(34S) Integration (using the double-exponential method)
This is my implementation of the double exponential integration method for the wp34s.

As stated in other threads, this method copes with discontinuities of the integrand or its derivatives at the integration interval ends, accepts infinities as interval ends and is, usually, much faster than the Romberg method.

There are here 3 versions of this program (the 3 versions are not identical in behaviour):
  1. The 1st one (v1.2r) is suited to be entered by hand (uses local labels for long jumps)
  2. The 2nd one (v1.0l) must be assembled (contains JMPs and GSBs) and is suited to be stored in RAM or LIB
  3. The 3rd one (v1.3x) is aimed to replace the built in integrator, must be build with the whole firmware


This is the version suited to be keyed in by hand.
Code:

// Double Exponential Integration for the wp34s calculator
//
// v1.2r-393 (20170327) by M. César Rodríguez GGPL
//
// USAGE:
//  Inputs:
//    - integration limits in X (upper limit) and Y (lower limit)
//    - name of integrand program in the alpha register
//  Outputs:
//    - approximate value of integral in X
//    - estimated error in Y
//    - upper integration limit in Z
//    - lower integration limit in T
//    - updates L (with upper limit)
//    - may corrupt ABCD
//
//  Remarks:
//    - accepts +infinite and -infinite as integration limits. If none
//      of the limits is infinite, the method applied is the tanh-sinh
//      one. If both are infinite then the sinh-sinh method is selected.
//      Otherwise the exp-sinh method is used
//    - the program tries to get all displayed digits correct, so
//      increasing the number of displayed digits will (hopefully) give
//      closer approximations (needing more time). But, in DBLON mode, the
//      program ignores the number of displayed digits and does its best
//      to get as much correct digits as possible
//    - during program execution the succesive aproximations are
//      displayed (as with the built-in integration program). Pressing
//      [<-] will interrupt the program & the last approximation
//      will be returned. The estimated error will be 0 in this case
//    - if the computed error estimation is not much smaller than the computed
//      result, it is assumed that all digits of the result are corrupted by
//      roundoff. In such cases, the reported result is 0 and the reported
//      error estimation equals the sum of the absolute values of the computed
//      result and error. This usually happens when the integral evaluates to 0
//    - all internal computations are carried out with flag D set (no
//      matter the value the user set to it), thus many errors when
//      evaluating the integrand (say 0/0 or 1/0, overflows or ...) will
//      not raise an error, but return infinite or NaN instead. Such
//      conditions will be detected by this program that will continue
//      the computation simply ignoring the offending point
//    - keep in mind that, when both integration limits are infinite (or
//      one is infinite and the other 0) the integrand may be evaluated
//      at really bigs/small points (up to 1e±201 in DBLOFF mode and
//      1e±3088 in DBLON mode). Ensure that your integrand program
//      behaves as expected for such arguments (remember that this program
//      works with flag D set, as stated previously)
//    - the double exponential method relies on the function to be
//      integrated being analytic over the integration interval, except,
//      perhaps, at the interval ends. Thus, if it is known that the
//      integrand, _or any of its derivatives_, has a discontinuity at some
//      point inside the integration interval, it is advised to break the
//      interval at such point and compute two separate integrals, one
//      at each side of the problematic point. Better results will be
//      get this way. Thus, beware of abs, fp, ip and such non analytic
//      functions inside the integrand program. Other methods (Romberg)
//      may behave better in this respect, but, on the other side, the
//      double exponential method manages discontinuities (of the
//      integrand or its derivatives) at the interval edges better and
//      is usually way faster that the Romberg method
//    - as many other quadrature methods (Romberg included), the
//      double exponential algorithm does not manage well highly
//      oscillatory integrands. Here, highly oscillatory means that the
//      integrand changes sign many (or infinite) times along the
//      integration interval
//
//  Method used:
//    The algorith used here is adapted from that described in:
//      H. Takahasi and M. Mori. "Double Exponential Formulas for Numerical
//      Integration." Publications of RIMS, Kyoto University 9 (1974),
//      721–741.
//    Another useful reference may be:
//      David H. Bailey, Xiaoye S. Li and Karthik Jeyabalan, "A comparison
//      of three high-precision quadrature schemes," Experimental
//      Mathematics, vol. 14 (2005), no. 3, pg 317-329.
//
// 275 steps
// 16 local registers

        LBL'DEI'          // Double-Exponential Integration
                INTM?       // check for invalid mode
                  ERR 13    // ERR_BAD_MODE
                LocR 16     // Local registers and flags:

                // local registers
                // .00    XB      backup of X (b)
                // .01    YB      backup of Y (a)
                // .02    bma2    (b - a)/2, a & b are the integration limits
                // .03    bpa2    (b + 2)/2
                // .04    eps     epsilon
                // .05    thr     convergence threshold
                // .06    lvl     level counter
                // .07    tm      maximum abscissa in the t domain
                // .08    h       space between abscissas in the t domain
                // .09    ssp     partial sum at each level
                // .10    j       abscissa counter
                // .11    ch      cosh(t)
                // .12    w       weight
                // .13    rp      abscissa (r) or fplus + fminus (p)
                // .14    ss      value of integral
                // .15    ss1     previous ss

                // local flags
                // F.00   DF      backup of D flag
                // F.01   rev     b is < a
                // F.02   ES      expsinh mode
                // F.03   bpa2z   bpa2 is zero in sinhsinh or expsinh modes
                // F.04   left    expsinh case with -infinite limit
                // F.05   TS      tanhsinh mode
                // F.06   lg0     true if level > 0

                // save stack **************************************************
                FS?S D      // backup & set flag D, specials are not errors here
                  SF .00
                cSTO .00    // save XY (ba)
                // A..D[->]    // save ABCD
                // check arguments  ************************************
                NaN?        // check for invalid limits
                  GTO 95    // b is NaN, exit
                x<>y
                NaN?
                  GTO 95    // a is NaN, exit
                x=? Y       // check for equal limits
                  GTO 90    // a == b, return 0
                // set mode flags, bma2 and bpa2  **********************
                x>? Y       // a > b?
                  SF .01    // yes, flag it
                INF?        // a == ±inf?
                  SF .02    // yes, flag it (temporarily in ES)
                x<>y
                INF?        // b == ±inf?
                  SF .03    // yes, flag it (temporarily in bpa2z)
                // set the ES, bpa2z, left & TS flags according to the
                // a & b values ++++++++++++++++++++++++++++++++++++++++
                FS? .02     // was a == ±inf?
                  FC? .03   // was b != ±inf?
                SKIP 004    // skip if a != ±inf or b != ±inf
                # 000       // here in the sinhsinh case  --------------
                # 002       // 2*bpa2 = 0, 2*bma2 = 2
                CF .02      // as this flag was dirty
                GTO 10      // done with sinhsinh
                FC? .02     // was a != ±inf?
                  FS? .03   // was b == ±inf?
                SKIP 007    // skip if a == ±inf or b == ±inf
                cENTER      // here in the tanhsinh case  --------------
                +           // 2*bpa2 ready
                x<> Z
                -
                ABS         // 2*bma2 ready
                SF .05      // flag the tanhsinh case
                GTO 10      // done with tanhsinh
                INF?        // here in the expsinh case ----------------
                  x<>y      // now X is the finite limit
                x>? Y       // finite limit > infinite one?
                  SF .04    // yes, left case
                RCL+ X      // 2*bpa2 ready
                CF .03      // as this flag was dirty
                x=0?        // finite limit == 0?
                  SF .03    // yes, flag it
                SF .02      // flag the expsinh case
                # 002       // 2*bma2  ready
        LBL 10
                # 002       // here with 2*bpa & 2*bma2 on stack
                /           // and mode flags ready
                STO .02     // save bma2
                x<>y
                RCL/ L
                STO .03     // save bpa2
                // done with mode flags, bpa2 and bma2  ****************
                // compute precision parameters ************************
                # 034       // 34 digits in DBLON mode
                +/-
                DBL?        // in DBL mode use all digits
                  SKIP 007
                # 009       // otherwise use displayed digits + 3,
                1/x         // let's get the number of displayed digits
                ROUND
                RCL- L
                EXPT        // now X is -(number of displayed digits + 1)
                # 002
                -           // -(the number of displayed digits + 3)
                ENTER       // here X = -(number of displayed + 3) (or -34)
                10^x
                STO .04     // save epsilon
                SQRT        // compute the convergence threshold, .05,
                SDL 001     // = 10*sqrt(epsilon)
                STO .05     // save the convergence threshold
                DROP
                ABS         // compute the maximum level
                LB
                ROUNDI
                # 002
                +           // maxlevel = round(log2(digits)) + 2
                STO .06     // save level
                FC? .03     // the sinhsinh case and the expsinh case
                  SKIP 004  // when the finite limit == 0 can use a
                # 001       // smaller epsilon
                # 000
                NEIGHB
                STO .04     // save epsilon for such cases
                // precision parameters (epsilon, thr, level) ready ****
                // compute maximum t value (tm) ************************
                // tm is computed in such a way that all abscissas are
                // representable at the given precision and that the
                // integrand is never evaluated at the limits
                FC? .05     // tanhsinh mode?
                  SKIP 004  // no, check other modes
                # 001       // start tm computation for tanh-sinh mode
                RCLMIN .02
                RCL+ X      // X = 2*MIN(1, bma2)
                SKIP 008    // goto remaining tm calculation
                FS? .03     // neither SS nor ES with 0 limit?
                  SKIP 004  // no, must be SS or ES with 0 limit
                # 1/2       // start tm computation for ES with != 0 limit
                RCL/ .03    // RCP/ bpa2
                ABS         // X = ABS(1/(2*bpa2))
                SKIP 002    // goto remaining tm calculation
                RCL .04     // start tm computation for SS & ES0 cases
                SQRT        // X = SQRT(eps)
                RCL/ .04    // continue tm computation for all cases
                LN
                FC? .05     // not in TS mode?
                  RCL+ X    // 2*ln(...) for all cases except TS
                RCL+ X      // 4*ln(...) for all cases except TS (2*LN(...))
                PI
                /
                LN
                STO .07     // tm done
                // maximum t (tm) ready ********************************
                // level loop ******************************************
                # 002
                STO .08     // h = 2
        LBL 20
                # 000
                STO .09     // ssp = 0
                # 001
                STO .10     // j = 1
                # 002
                STO/ .08    // h /= 2
                RCL .08     // X = t
                // j loop ++++++++++++++++++++++++++++++++++++++++++++++
                // compute abscissas and weights  ----------------------
        LBL 30
                COSH        // cosh(t) (cosh is much faster than sinh/tanh)
                STO .11     // save for later
                x^2
                DEC X
                SQRT
                PI
                *
                # 002
                /           // pi/2*sqrt(cosh(t)^2 - 1) = pi/2*sinh(t)
                FS? .02     // ES mode?
                  SKIP 002  // yes, want EXP
                COSH        // no, want COSH
                SKIP 001
                EXP
                FILL
                FC? .05     // TS mode?
                  SKIP 002  // no, done with weight (for ES, SS cases)
                x^2         // end weight computation for TS mode
                1/x
                STO .12     // save weight
                DROP        // stack filled with exp/cosh of pi/2*sinh(t)
                FS? .02     // ES mode?
                  SKIP 003  // yes, skip this part
                x^2         // no ES mode, get r (the abscissa in the
                DEC X       // normalized domain)
                SQRT
                FS? .05     // TS mode?
                  RCL/ Y    // yes, adjust r
                FS? .04     // ES mode -infinity?
                  +/-       // yes, adjust r
                STO .13     // save normalized abscissa
                // done with abscissas and weights  --------------------
                // evaluate integrand ----------------------------------
                RCL* .02    // r*(b - a)/2
                RCL+ .03    // (b + a)/2 + r*(b - a)/2
                FILL
                XEQa        // f(bpa2 + bma2*r)
                SPEC?       // do not stop in error (flag D is set)
                  # 000
                RCL* .12    // fplus*w
                x<> .13     // p = fplus*w stored, r in X
                RCL .03     // (b + a)/2
                RCL .02     // (b - a)/2
                FS? .02     // ES mode?
                  SKIP 003
                RCL* Z      // no ES mode, "normal" abscissa
                -
                SKIP 002
                RCL/ Z      // ES mode, "inverse" abscissa
                +           // X = (b + a)/2 ± (b - a)/2*(r or 1/r)
                FILL
                XEQa        // f(bpa2 ± bma2*(r or 1/r))
                SPEC?       // do not stop in error (flag D is set)
                  # 000
                FS? .02     // ES mode?
                  SKIP 002
                RCL* .12    // no ES mode, normal weight, fminus*w
                SKIP 001
                RCL/ .12    // ES mode, inverse weight, fminus/w
                RCL+ .13    // p = fplus*w + fminus*(w or 1/w)
                RCL* .11    // X = ch*p
                STO+ .09    // ssp += p
                // done with integrand  --------------------------------
                // level 0 is special ----------------------------------
                FC? .06     // is level == 0?
                  GTO 40    // yes, go to next level
                INC .10     // no, sweep odd abscissas, j += 1
                ABS         // X = ABS(p)
                RCL .09     // ssp
                RCL* .04    // eps
                ABS         // X = ABS(ssp*eps), Y = ABS(p)
                x>=? Y      // ABS(p) <= ABS(eps*eps)?
                  GTO 50    // then done with level
                // and check user interrupt ............................
                KEY? Y      // any key?
                  SKIP 005  // no, continue
                # 035       // yes, check keycode for <-
                x!=? Z      // key is not [<-]?
                  SKIP 002  // it is not [<-], continue
                RCL .15     // it is, recall the last computation
                GTO 70      // and exit
                // if level > 0 done  ----------------------------------
        LBL 40
                INC .10     // j += 1
                RCL .10     // X = j
                RCL* .08    // X = t = j*h
                x<=? .07    // t <= tm?
                  GTO 30    // yes, continue j loop
                // done with j loop ++++++++++++++++++++++++++++++++++++
        LBL 50
                RCL .09     // no, update ss. X = ssp
                STO+ .14    // ss += ssp
                FS? .06     // is level > 0 ----------------------------
                  GTO 60    // yes, progress to next level
                # 001       // no, add 1st series term
                FS? .04     // left?
                  +/-
                RCL .03     // (b + a)/2
                FS? .02     // ES mode?
                  +
                FILL
                XEQa        // f(bpa2 (+/-1 in ES mode))
                SPEC?       // do not stop in error (flag D is set)
                  # 000
                STO+ .14    // ss += f(bpa2)
                // done with level 0  ----------------------------------
                // apply constant coeffs to sum ------------------------
        LBL 60
                RCL .14
                RCL* .02    // ss*bma2
                RCL* .08    // ss*bma2*h
                PI
                *           // ss*bma2*h*pi
                # 002
                /           // ss*bma2*h*pi/2
                FS? .01     // reverse?
                  +/-       // yes,so change sign
                // done with constant coeffs  --------------------------
                // show progress  --------------------------------------
                TOP?        // show value of approximation
                  MSG 25    // MSG_INTEGRATE
                // check convergence  ----------------------------------
                RCL .09
                RCL+ X      // 2*ssp
                RCL- .14    // 2*ssp - ss
                ABS         // ABS(2*ssp - ss)
                RCL .05
                RCL* .14    // thr*ss
                ABS         // ABS(thr*ss)
                x>? Y       // ABS(2*ssp - ss) < ABS(thr*ss)?
                  SKIP 004  // yes, computation done
                z<> .15     // no, save ss (with all coeffs applied)
                SF .06      // mark level 0 done,
                DSL .06     // update level &...
                  GTO 20    // loop.
                RCL Z       // recall result
                // done with level loop ********************************
                // compute error estimation and check for bad results **
                // here with result in X
        LBL 70
                RCL X       // stack: ss-ss-?-?
                ABS         // stack: |ss|-ss-?-?
                RCL Y       // stack: ss-|ss|-ss-?
                RCL- .15    // stack: (ss-ss1)-|ss|-ss-?
                ABS         // stack: err-|ss|-ss-?
                RCL X       // stack: err-err-|ss|-ss
                SDL 001     // stack: 10*err-err-|ss|-ss
                x<? Z       // 10*err < |ss|?
                  SKIP 004  // yes, assume result is OK
                x<> Z       // no, bad result. stack: |ss|-err-10*err-ss
                +           // stack: new_err-10*err-ss-ss
                # 000       // stack: 0-new_err-10*err-ss
                STO T       // stack: 0-new_err-10*err-0
                x<> T       // stack: 0-err-... or ss-err-...
                // done with bad results  ******************************
                //  exit  **********************************************
        LBL 80
                TOP?        // show new computed value (0)
                  MSG 25    // MSG_INTEGRATE
                cRCL .00    // stack: b-a-s-ss1
                // [->]A..D    // try to restore ABCD
                STO L       // b into L
                cx<> Z      // stack: s-ss1-b-a
                FC? .00     // restore user's D flag (DF)
                  CF D
              RTN           // bye
                // exit cause a == b  **********************************
        LBL 90
                # 000       // error & result both 0
                SKIP 002
                // exit cause a or b is NaN ****************************
        LBL 95
                # NaN       // error & result both NaN
                RCL X
                GTO 80      // exit
              END
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
(34S) Integration (using the double-exponential method) - emece67 - 03-24-2017 02:21 PM



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