[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
 emece67 Senior Member Posts: 379 Joined: Feb 2015
[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.
 « Next Oldest | Next Newest »

 Messages In This Thread [wp34s] New integration program - emece67 - 03-08-2017 05:46 PM RE: [wp34s] New integration program - emece67 - 03-10-2017, 01:35 PM RE: [wp34s] New integration program - Paul Dale - 03-11-2017, 12:22 PM RE: [wp34s] New integration program - Massimo Gnerucci - 03-11-2017, 12:45 PM RE: [wp34s] New integration program - emece67 - 03-11-2017, 03:10 PM RE: [wp34s] New integration program - emece67 - 03-12-2017, 11:07 PM RE: [wp34s] New integration program - emece67 - 03-17-2017, 06:16 PM RE: [wp34s] New integration program - nova - 08-11-2019, 05:33 AM RE: [wp34s] New integration program - emece67 - 04-03-2021, 01:40 PM RE: [wp34s] New integration program - emece67 - 03-24-2017, 05:12 PM RE: [wp34s] New integration program - Albert Chan - 04-04-2021, 03:33 PM RE: [wp34s] New integration program - emece67 - 04-04-2021, 05:09 PM

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