[wp34s] New integration program

[emece67]

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.

[emece67]

Find attached my findings when comparing this new program against the built in integration program.


  Tanh-Sinh_comparison.pdf (Size: 21.23 KB / Downloads: 72)

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.

[Paul Dale]

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

[Massimo Gnerucci]

(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.

[emece67]

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.

[emece67]

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.

[emece67]

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.

[emece67]

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.

[nova]

(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

[emece67]

(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.

[Albert Chan]

(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 [-∞,+∞]

---

>>> quadde(f, [-40,40])
(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)
>>> quadde(fold, [0,40])
(mpf('0.81659478386385076'), mpf('2.0486468077507425e-10'), 161, 5, 0)

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

>>> quadde(fold, [0,inf])
(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

[emece67]

(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.