+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: Not HP Calculators (/forum-7.html)
+--- Forum: Not quite HP Calculators - but related (/forum-8.html)
+--- Thread: [wp34s] New integration program (/thread-7905.html)
[wp34s] New integration program - emece67 - 03-08-2017 05:46 PM
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
ENDNOTE: 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.
RE: [wp34s] New integration program - emece67 - 03-10-2017 01:35 PM
Find attached my findings when comparing this new program against the built in integration program.
[attachment=4577]
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.
RE: [wp34s] New integration program - Paul Dale - 03-11-2017 12:22 PM
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
RE: [wp34s] New integration program - Massimo Gnerucci - 03-11-2017 12:45 PM
(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.
RE: [wp34s] New integration program - emece67 - 03-11-2017 03:10 PM
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.
RE: [wp34s] New integration program - emece67 - 03-12-2017 11:07 PM
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.
RE: [wp34s] New integration program - emece67 - 03-17-2017 06:16 PM
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.
RE: [wp34s] New integration program - emece67 - 03-24-2017 05:12 PM
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.
RE: [wp34s] New integration program - nova - 08-11-2019 05:33 AM
(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
RE: [wp34s] New integration program - emece67 - 04-03-2021 01:40 PM
(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.
RE: [wp34s] New integration program - Albert Chan - 04-04-2021 03:33 PM
(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
RE: [wp34s] New integration program - emece67 - 04-04-2021 05:09 PM
(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.