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.
Search
Member List
Calendar
Help
/
[wp34s] New integration program
