Post Reply 
[wp34s] New integration program
03-08-2017, 05:46 PM (This post was last modified: 03-13-2017 05:59 AM by emece67.)
Post: #1
[wp34s] New integration program
Hi all,

In another thread I presented a description of an algorithm (not mine) to approximate definite integrals different to the Romberg method classically used by hp machines: the tanh-sinh method.

After some time I was able to write a program for the wp34s machine that uses this method to compute definite integrals. Find it below.

In my tests this method is much faster than the Romberg one & seems to manage appropriately a wider class of integrands (infinite integrands at one or both interval edges, infinite derivatives at one or both interval ends).

I really think that this method may replace the Romberg one even in small handheld calculators, but it is the first implementation of such method I know in a calculator & much more tests are needed to check its behaviour, so beta testers are welcome.

Regards.

Code:

// Tanh-Sinh Integration for the wp34s calculator
//
// v1.0l-372 (20170313) by M. César Rodríguez GGPL
//
// USAGE:
//  Inputs:
//    - integration limits on X (upper limit) and Y (lower limit)
//    - name of integrand program in the alpha register
//  Outputs:
//    - approximate value of integral on X
//    - updates L
//    - corrupts stack (works in both SSIZE8 and SSIZE4)
//  Remarks:
//    - 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 after & the last approximation
//      will be returned
//
// TODO: document
//
// 203 steps (118)
// 15 local registers (23)

              LBL'TSI'      // Tanh-Sinh Integration
                INTM?       // check for invalid mode
                  ERR 13    // ERR_BAD_MODE
                LocR 15     // Local registers and flags:
                            //  .00       X backup (to return in L)
                            //  .01       bpa2  ((b + 2)/2)
                            //  .02       bma2  ((b - a)/2)
                            //  .03       thr   (convergence threshold)
                            //  .04       level (counter)
                            //  .05       jmax  (maximum index on each level)
                            //  .06       Nmax  (maximum index on all levels)
                            //  .07       step  (j step size for each level)
                            //  .08       hmin  (minimum space between
                            //                  abscissas in the t domain)
                            //  .09       ss    (value of integral)
                            //  .10       ssp   (partial sum at each level)
                            //  .11       j     (counter)
                            //  .12       w     (weight)
                            //  .13       stores f(bpa2 + bma2*r) and
                            //            (bpa2 + bma2*r) alternatively
                            //  .14       fname (name of integrand program)
                            //  F.00      backup of D flag


                STO .00     // save X (for L at exit)
                SPEC?       // check for invalid limits
                  JMP TSI_bad_rng       // SKIP 179
                x<>y
                SPEC?
                  JMP TSI_bad_rng       // SKIP 176
                x=? Y       // check for equal limits
                  JMP TSI_0             // SKIP 171
                FS?S D      // backup & set flag D, specials are not errors here
                  SF .00
                aSTO .14    // aSTO fname. Name of function to integrate
                cENTER      // let's compute bpa2 & bma2
                +
                # 002
                /
                STO .01     // bpa2 = (b + a) / 2
                DROP
                -
                # 002
                /
                STO .02     // bma2 = (b - a) / 2. Done
                # 34        // compute precision parameters
                +/-
                DBL?        // in DBL mode use all digits
                  JMP TSI_mps           // 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 -the number of significant digits + 1
                # 002
                -           // -the number of significant digits + 3
TSI_mps::       ENTER
                10^x
                STO .06     // epsilon
                SQRT
                SDL 003
                STO .03     // convergence threshold = 1000 * sqrt(epsilon)
                DROP
                ABS
                LB
                ROUNDI
                # 002
                +
                # 006
                MIN         // maxlevel must be < 6 cause of DSE counter step
                STO .04     // maxlevel = round(log2(digits)) + 2
                2^x
                STO .07     // step = 2^maxlevel
                1/x
                STO .08     // hmin = 2^-maxlevel
                # 001       // let's compute the maximum t value that yields
                RCL .02     // {r, x} values strictly below {1, b} when
                ABS         // working at given precission,
                MIN         // t ~ asinh(ln(2*min(1, abs(bma2))/eps)/pi)
                RCL+ X
                RCL/ .06    // RCL/ epsilon
                LN
                PI
                /
                ASINH       // done with t
                RCL* .07    // RCL* step
                FLOOR
                STO .06     // Nmax (tmax = Nmax*step)
                RCL .07     // RCL
                IDIV
                RCL* L
                STO .05     // jmax = step * (Nmax // step)
                # 000
                STO .09     // ss = 0
                RCL .04     // RCL maxlevel
                SDR 003
                STO .04     // level loop, use ISG (0.maxlevel)
TSI_lvl_loop::  # 000
                STO .10     // ssp = 0
                RCL .07     // step
                SDR 005     // step/100000
                RCL+ .05    // jmax + step/100000
                STO .11     // j counter, use DSE (jmax.000step)
TSI_j_loop::    RCL .11     // RCL j
                IP          // discard cccff DSE fields
                RCL* .08    // t = j*hmin
                COSH        // cosh(t)
                ENTER
                x^2
                DEC X
                SQRT        // sqrt(cosh^2(t)) - 1 = sinh(t)
                PI
                *
                # 002
                /           // pi/2*sinh(t)
                cENTER
                COSH        // cosh(pi/2*sinh(t))
                x^2         // cosh^2(pi/2*sinh(t))
                /           // cosh(t)/cosh^2(pi/2*sinh(t)) = w
                STO .12     // w
                RCL L       // cosh^2(pi/2*sinh(t))
                1/x         // 1/cosh^2(pi/2*sinh(t))
                # 001
                x<>y
                -           // 1 - 1/cosh^2(pi/2*sinh(t))
                SQRT        // sqrt(...) = tanh(pi/2*sinh(t)) = r
                RCL* .02    // r*(b - a)/2
                STO .13     // save r*(b - a)/2
                RCL+ .01    // (b + a)/2 + r*(b - a)/2 = x
                FILL
                aXEQ .14    // f(bpa2 + bma2*r)
                SPEC?       // do not stop in error (flag D is set)
                  # 000
                x<> .13     // save f(xp), restore r*(b - a)/2
                RCL .01     // bpa2
                x<> Y
                -           // (b + a)/2 - r*(b - a)/2 = x
                FILL
                aXEQ .14    // f(bpa2 - bma2*r)
                SPEC?
                  # 000
                RCL+ .13    // f(xp) + f(xm)
                RCL* .12    // w*(f(xp) + f(xm))
                STO+ .10    // ssp += w*(f(fp) + f(fm))
                // check user interrupt
                KEY? Y      // any key? No, continue
                  JMP TSI_nok           // SKIP 005
                # 035       // keycode for <-
                x!=? Z      // key is not [<-]?
                  SKIP 002  // it is not [<-], continue
                RCL .09     // it is [<-], recall ss and exit
                JMP TSI_converged       // SKIP 045
TSI_nok::       DSE .11     // j loop
                JMP TSI_j_loop          // BACK 049
                RCL .10     // ssp
                STO+ .09    // ss += ssp
                RCL .04     // level
                # 001
                x>? Y       // level == 0?
                  JMP TSI_l0            // SKIP 003
                # 002       // no
                STO/ .07    // step /= 2, levels 0 & 1 use the same step
                JMP TSI_chk_conv        // SKIP 006
TSI_l0::        RCL .01     // yes, bpa2
                FILL
                aXEQ .14    // f(bpa2)
                SPEC?
                  # 000
                STO+ .09    // ss += f(bpa2) (only at level 0)
TSI_chk_conv::  GSB TSI_factors         // BSRF 044
                TOP?        // show value of approximation
                  MSG 25    // MSG_INTEGRATE
                RCL .10     // ssp
                x=+0?
                  JMP TSI_converged     // SKIP 022
                x=-0?
                  JMP TSI_converged     // SKIP 020
                # 002
                *           // 2*ssp
                RCL/ .09    // 2*ssp/ss
                DEC X       // 2*ssp/ss - 1
                ABS         // abs(2*ssp/ss - 1)
                RCL .03     // threshold
                x>? Y       // converged? (abs(2*ssp/ss - 1) < threshold)
                  JMP TSI_converged     // SKIP 012
                RCL .06     // no, compute next loop. Nmax
                RCL .07     // step
                # 002
                /
                STO .05     // step/2
                -           // Nmax - step/2
                RCL .07     // step
                IDIV        // (Nmax - step/2) // step
                RCL* .07    // step * ((Nmax - step/2) // step)
                STO+ .05    // jmax = step/2 + step * ((Nmax - step/2) // step)
                ISG .04     // level loop
                  JMP TSI_lvl_loop      // BACK 098
TSI_converged:: GSB TSI_factors         // BSRF 016
                FC? .00     // restore user's D flag
                  CF D
                RCL .00     // RCL initial X
                STO L
                x<>y
              RTN
TSI_0::         DROP        // here if limits are equal
                # 000       // result is 0
                SKIP 002
TSI_bad_rng::   DROP        // here if limit is special
                # NaN       // result is NaN
                y<> .00     // restore argument
                y<> L       // save in L
                x<> Y       // restore stack
                DROP
              RTN           // all done
TSI_factors::     RCL .09   // multiply sum by constant coeffs to get approx.
                  RCL* .02  // ss*bma2
                  RCL* .08  // ss*bma2*hmin
                  RCL* .07  // ss*bma2*hmin*step
                  PI
                  *         // ss*bma2*hmin*step*pi
                  # 002
                  /         // ss*bma2*hmin*step*pi/2
                RTN
            END

NOTE: updated to v1.0l (some tweaks in order to speed-up execution, get more precision in some difficult cases and a more responsive user interrupt)

Note 2: updated to work in DBLON mode.

César - Information must flow.
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
[wp34s] New integration program - emece67 - 03-08-2017 05:46 PM
RE: [wp34s] New integration program - nova - 08-11-2019, 05:33 AM



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