[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: 378 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.
03-10-2017, 01:35 PM (This post was last modified: 03-12-2017 10:43 PM by emece67.)
Post: #2 emece67 Senior Member Posts: 378 Joined: Feb 2015
RE: [wp34s] New integration program
Find attached my findings when comparing this new program against the built in integration program. Tanh-Sinh_comparison.pdf (Size: 21.23 KB / Downloads: 67)

Finally I used a battery of 40 test integrals of various kinds.

On this table:
• "Better" means that the new program is more precise and faster
• "More precise" means that it is more precise but slower
• "Faster" means that it is equally precise but faster
• "Equal" means that it is both equally precise and fast
• "Slower" means that it is equally precise but slower
• "Less precise" means that it is less precise but faster
• "Worse" means that it is less precise and slower

When counting the correct digits I was only concerned about the 12 visible digits (if the expected result is 0, the table shows the absolute value of the exponent of the returned result). Time was measured with a handheld chrono, so don't expect high accuracy on such. Times marked as ">xxx" means that I stopped measuring the time when xxx seconds get reached. "est. time" means that I estimated the time required for the built-in integrator by measuring a few iterations on the real machine and using the simulator to get the number of total iterations needed. The machine was always in SLOW and DBLOFF modes.

Regards.

NOTE: updated to show results of last v1.0l version. Some typos fixed.
03-11-2017, 12:22 PM
Post: #3 Paul Dale Senior Member Posts: 1,750 Joined: Dec 2013
RE: [wp34s] New integration program
Any chance of adding the WP43 version 2.2 firmware to the results? The integration there was a non-adaptive Gauss-Kronrod 10/21 point quadrature -- it is exact for polynomials up to degree 19, although subject to numeric stability, but not so great for other functions.

Pauli
03-11-2017, 12:45 PM
Post: #4 Massimo Gnerucci Senior Member Posts: 2,446 Joined: Dec 2013
RE: [wp34s] New integration program
(03-11-2017 12:22 PM)Paul Dale Wrote:  Any chance of adding the WP43 version 2.2 firmware to the results?

So WP43 firmware is already at version 2.2?!?

Great news!

Just kidding, SCNR.

Greetings,
Massimo

-+×÷ ↔ left is right and right is wrong
03-11-2017, 03:10 PM
Post: #5 emece67 Senior Member Posts: 378 Joined: Feb 2015
RE: [wp34s] New integration program
V2.2 used the 10/21 Gauss-Kronrod method. Current v3.3 uses the Romberg Method. The table I posted compares my implementation of the Tanh-Sinh method against the Romberg method in v3.3.
03-12-2017, 11:07 PM (This post was last modified: 03-13-2017 06:03 AM by emece67.)
Post: #6 emece67 Senior Member Posts: 378 Joined: Feb 2015
RE: [wp34s] New integration program
I've just updated the program listing and the comparison table cause I've added some tweaks to the program.

Overall I think that it, definitely, wins against the built-in, Romberg-based integration program in the v3.3 wp34.

The main source of trouble for this algorithm seems to be such functions with infinite derivatives inside the integration interval (even if the function does not go to infinity).

There are other cases where the built-in integrator is better, but even in such cases this tanh-sinh method gives acceptable results (F22) or, simply, both methods are failing miserably (F30, F32).

In all these problem cases, switching to DBLON mode does not help.

On the other side, the huge speed-up in almost any situation of this tanh-sinh method and its better results when the integrand goes to infinity at the integration interval edges makes it, IMHO, a safe replacement for the built-in integrator in the wp34 calculator and, perhaps, in other machines with Romberg based integrators.

Even more, with minor changes, this method can cope with integrals where one of the limits is infinite (exp-sinh method) or both are (sinh-sinh method).

Regards.
03-17-2017, 06:16 PM
Post: #7 emece67 Senior Member Posts: 378 Joined: Feb 2015
RE: [wp34s] New integration program
I've continued looking for ways to improve this program. Although I've not yet implemented it for the wp34s, I do now have an emulation, in Python, of a new version of the program that also manages not only finite integration intervals, but both infinite and semi-infinite ones.

I'll post the wp34s version in a few days.

Regards.
03-24-2017, 05:12 PM (This post was last modified: 03-27-2017 05:52 PM by emece67.)
Post: #8 emece67 Senior Member Posts: 378 Joined: Feb 2015
RE: [wp34s] New integration program
I've just posted an implementation of the double-exponential method (that is: tanh-sinh, exp-sinh or sinh-sinh, depending on the integration interval ends being finite or not) in the General Software library.

This new program has grown in size (about 50 steps for a total of ~250) but now it copes with infinite or semi-infinite integration intervals. It can now manage things like:
$\int_{-1}^\infty{\ln(x+1)\over x^2+1}dx={3\pi\ln(2)\over8}\approx0.8165947838638507989377583368391052\ldots$
that has a discontinuity at one interval end (-1), being the other end infinite, getting as result:
• 0.816594783863850 in DBLOFF mode in 14 seconds
• 0.8165947838638507989377583368390855 in DBLON mode in 35 seconds

Although there are, of course, integrands that may drive this method in trouble, I'm amazed about its performance and capabilities.

Regards.
08-11-2019, 05:33 AM
Post: #9
 nova Junior Member Posts: 9 Joined: Aug 2019
RE: [wp34s] New integration program
(03-17-2017 06:16 PM)emece67 Wrote:  I've continued looking for ways to improve this program. Although I've not yet implemented it for the wp34s, I do now have an emulation, in Python, of a new version of the program that also manages not only finite integration intervals, but both infinite and semi-infinite ones.

I'll post the wp34s version in a few days.

Regards.

Yes, it's been a while since the above message was placed!
The post discusses a new Python version which manages semi-infinite and infinite intervals as well as finite intervals.
Is the Python code listing available?
Thanks
nova
04-03-2021, 01:40 PM
Post: #10 emece67 Senior Member Posts: 378 Joined: Feb 2015
RE: [wp34s] New integration program
(08-11-2019 05:33 AM)nova Wrote:  Is the Python code listing available?

Hi,

As I have stated in another thread, my original Python code for the double-exponential quadrature is now on GitHub.

It has suffered some refactoring from its original shape (that I used to derive the wp34s code).

Regards.
04-04-2021, 03:33 PM (This post was last modified: 04-10-2021 10:45 PM by Albert Chan.)
Post: #11
 Albert Chan Senior Member Posts: 1,843 Joined: Jul 2018
RE: [wp34s] New integration program
(03-24-2017 05:12 PM)emece67 Wrote:  $\int_{-1}^\infty{\ln(x+1)\over x^2+1}dx={3\pi\ln(2)\over8}\approx0.8165947838638507989377583368391052\ldots$

With this integral transformed, I noted some weekness of tanh-sinh.
(Or, it might be just a fluke with this example ...)

Rewrite above integral, then let x=e^y:

$$\displaystyle \int _0 ^∞ {\ln(x) \over 1+(x-1)^2} dx = \int _{-∞} ^∞ {y\;e^y \over 1+(e^y-1)^2} dy$$

>>> from mpmath import *
>>> import double_exponential # emece67 source from here
>>> quadde = lambda f,(a,b): double_exponential.double_exponential(f,a,b)

>>> f = lambda y: y*exp(y) / (1 + expm1(y)**2)
>>> quadde(f, [-inf,inf]) # method sinh-sinh
(mpf('0.81659478386382767'), mpf('7.7731406533665393e-8'), 73, 4, 2)

To use tanh-sinh, just convert back intervals: [-∞,+∞] -> [0,1]

>>> def f01(t): y=t*t-t; return f((2*t-1)/y) * (2*y+1)/(y*y)
...
>>> quadde(f01, [0,1]) # method tanh-sinh
(mpf('0.81659478390520546'), mpf('1.0424940875997102e-5'), 77, 5, 0)

It seems sinh-sinh does better than tanh-sinh, for interval [-∞,+∞]

---

(0, mpf('0.83732598514658618'), 87, 5, 0)

A bit unexpected, getting area of 0.0, with huge error estimate.

plot(f, [-40,40]), showed it has a shape of a heartbeat (inside ±5)
We should instead move dominant part to the edge.

>>> fold = lambda x: f(x) + f(-x)
(mpf('0.81659478386385076'), mpf('2.0486468077507425e-10'), 161, 5, 0)

Result is similar to exp-sinh, for [0,+∞]

(mpf('0.81659478386385087'), mpf('4.1397899552819695e-9'), 131, 4, 1)

Edit: I misinterpreted last argument, should be 0=tanh-sinh, 1=exp-sinh, 2=sinh-sinh
04-04-2021, 05:09 PM (This post was last modified: 04-04-2021 05:16 PM by emece67.)
Post: #12 emece67 Senior Member Posts: 378 Joined: Feb 2015
RE: [wp34s] New integration program
(04-04-2021 03:33 PM)Albert Chan Wrote:  To use tanh-sinh, just convert back intervals: [-∞,+∞] -> [0,1]

>>> def f01(t): y=t*t-t; return f((2*t-1)/y) * (2*y+1)/(y*y)
...
>>> quadde(f01, [0,1]) # method tanh-sinh
(mpf('0.81659478390520546'), mpf('1.0424940875997102e-5'), 77, 5, 0)

It seems exp-sinh does better than tanh-sinh, for interval [-∞,+∞]

I suppose it depends on the particular change of variable you use to get from [-∞, +∞] to [0, 1]. Maybe you can find cases when the tanh-sinh method behaves better in [0, 1] than transforming the problem to the exp-sinh case.

(04-04-2021 03:33 PM)Albert Chan Wrote:  >>> quadde(f, [-40,40])
(0, mpf('0.83732598514658618'), 87, 5, 0)

A bit unexpected, getting area of 0.0, with huge error estimate.

This implementation of the double-exponential algorithm is aimed at a handheld calculator, so it uses simple heuristics in some places. Here you are fooling one of such heuristics (related to the fact that the function vanishes rapidly away from 0). The first iterations of the algorithm "doesn't see" the interesting part of the function, and when later iterations find something, the algorithm thinks it is round-off noise. In some way this function is "bell-shaped", and this algorithm concentrates the sampling points at the extremes of the integration interval, so caution is advised for cases such this one. Integrating exp(-x^2) in [-10, 10] or [-inf, +inf] gives good results, but in [-20, 20] triggers the heuristic.

You have found a way to circumvent this, though.

Regards.
 « Next Oldest | Next Newest »

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