(34S) Integration (using the double-exponential method)
This is my implementation of the double exponential integration method for the wp34s.
As stated in other threads, this method copes with discontinuities of the integrand or its derivatives at the integration interval ends, accepts infinities as interval ends and is, usually, much faster than the Romberg method.
There are here 3 versions of this program (the 3 versions are not identical in behaviour):
- The 1st one (v1.2r) is suited to be entered by hand (uses local labels for long jumps)
- The 2nd one (v1.0l) must be assembled (contains JMPs and GSBs) and is suited to be stored in RAM or LIB
- The 3rd one (v1.3x) is aimed to replace the built in integrator, must be build with the whole firmware
This is the version suited to be keyed in by hand.
Code:
// Double Exponential Integration for the wp34s calculator
//
// v1.2r-393 (20170327) by M. César Rodríguez GGPL
//
// USAGE:
// Inputs:
// - integration limits in X (upper limit) and Y (lower limit)
// - name of integrand program in the alpha register
// Outputs:
// - approximate value of integral in X
// - estimated error in Y
// - upper integration limit in Z
// - lower integration limit in T
// - updates L (with upper limit)
// - may corrupt ABCD
//
// Remarks:
// - accepts +infinite and -infinite as integration limits. If none
// of the limits is infinite, the method applied is the tanh-sinh
// one. If both are infinite then the sinh-sinh method is selected.
// Otherwise the exp-sinh method is used
// - 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 & the last approximation
// will be returned. The estimated error will be 0 in this case
// - if the computed error estimation is not much smaller than the computed
// result, it is assumed that all digits of the result are corrupted by
// roundoff. In such cases, the reported result is 0 and the reported
// error estimation equals the sum of the absolute values of the computed
// result and error. This usually happens when the integral evaluates to 0
// - all internal computations are carried out with flag D set (no
// matter the value the user set to it), thus many errors when
// evaluating the integrand (say 0/0 or 1/0, overflows or ...) will
// not raise an error, but return infinite or NaN instead. Such
// conditions will be detected by this program that will continue
// the computation simply ignoring the offending point
// - keep in mind that, when both integration limits are infinite (or
// one is infinite and the other 0) the integrand may be evaluated
// at really bigs/small points (up to 1e±201 in DBLOFF mode and
// 1e±3088 in DBLON mode). Ensure that your integrand program
// behaves as expected for such arguments (remember that this program
// works with flag D set, as stated previously)
// - the double exponential method relies on the function to be
// integrated being analytic over the integration interval, except,
// perhaps, at the interval ends. Thus, if it is known that the
// integrand, _or any of its derivatives_, has a discontinuity at some
// point inside the integration interval, it is advised to break the
// interval at such point and compute two separate integrals, one
// at each side of the problematic point. Better results will be
// get this way. Thus, beware of abs, fp, ip and such non analytic
// functions inside the integrand program. Other methods (Romberg)
// may behave better in this respect, but, on the other side, the
// double exponential method manages discontinuities (of the
// integrand or its derivatives) at the interval edges better and
// is usually way faster that the Romberg method
// - as many other quadrature methods (Romberg included), the
// double exponential algorithm does not manage well highly
// oscillatory integrands. Here, highly oscillatory means that the
// integrand changes sign many (or infinite) times along the
// integration interval
//
// Method used:
// The algorith used here is adapted from that described in:
// H. Takahasi and M. Mori. "Double Exponential Formulas for Numerical
// Integration." Publications of RIMS, Kyoto University 9 (1974),
// 721–741.
// Another useful reference may be:
// David H. Bailey, Xiaoye S. Li and Karthik Jeyabalan, "A comparison
// of three high-precision quadrature schemes," Experimental
// Mathematics, vol. 14 (2005), no. 3, pg 317-329.
//
// 275 steps
// 16 local registers
LBL'DEI' // Double-Exponential Integration
INTM? // check for invalid mode
ERR 13 // ERR_BAD_MODE
LocR 16 // Local registers and flags:
// local registers
// .00 XB backup of X (b)
// .01 YB backup of Y (a)
// .02 bma2 (b - a)/2, a & b are the integration limits
// .03 bpa2 (b + 2)/2
// .04 eps epsilon
// .05 thr convergence threshold
// .06 lvl level counter
// .07 tm maximum abscissa in the t domain
// .08 h space between abscissas in the t domain
// .09 ssp partial sum at each level
// .10 j abscissa counter
// .11 ch cosh(t)
// .12 w weight
// .13 rp abscissa (r) or fplus + fminus (p)
// .14 ss value of integral
// .15 ss1 previous ss
// local flags
// F.00 DF backup of D flag
// F.01 rev b is < a
// F.02 ES expsinh mode
// F.03 bpa2z bpa2 is zero in sinhsinh or expsinh modes
// F.04 left expsinh case with -infinite limit
// F.05 TS tanhsinh mode
// F.06 lg0 true if level > 0
// save stack **************************************************
FS?S D // backup & set flag D, specials are not errors here
SF .00
cSTO .00 // save XY (ba)
// A..D[->] // save ABCD
// check arguments ************************************
NaN? // check for invalid limits
GTO 95 // b is NaN, exit
x<>y
NaN?
GTO 95 // a is NaN, exit
x=? Y // check for equal limits
GTO 90 // a == b, return 0
// set mode flags, bma2 and bpa2 **********************
x>? Y // a > b?
SF .01 // yes, flag it
INF? // a == ±inf?
SF .02 // yes, flag it (temporarily in ES)
x<>y
INF? // b == ±inf?
SF .03 // yes, flag it (temporarily in bpa2z)
// set the ES, bpa2z, left & TS flags according to the
// a & b values ++++++++++++++++++++++++++++++++++++++++
FS? .02 // was a == ±inf?
FC? .03 // was b != ±inf?
SKIP 004 // skip if a != ±inf or b != ±inf
# 000 // here in the sinhsinh case --------------
# 002 // 2*bpa2 = 0, 2*bma2 = 2
CF .02 // as this flag was dirty
GTO 10 // done with sinhsinh
FC? .02 // was a != ±inf?
FS? .03 // was b == ±inf?
SKIP 007 // skip if a == ±inf or b == ±inf
cENTER // here in the tanhsinh case --------------
+ // 2*bpa2 ready
x<> Z
-
ABS // 2*bma2 ready
SF .05 // flag the tanhsinh case
GTO 10 // done with tanhsinh
INF? // here in the expsinh case ----------------
x<>y // now X is the finite limit
x>? Y // finite limit > infinite one?
SF .04 // yes, left case
RCL+ X // 2*bpa2 ready
CF .03 // as this flag was dirty
x=0? // finite limit == 0?
SF .03 // yes, flag it
SF .02 // flag the expsinh case
# 002 // 2*bma2 ready
LBL 10
# 002 // here with 2*bpa & 2*bma2 on stack
/ // and mode flags ready
STO .02 // save bma2
x<>y
RCL/ L
STO .03 // save bpa2
// done with mode flags, bpa2 and bma2 ****************
// compute precision parameters ************************
# 034 // 34 digits in DBLON mode
+/-
DBL? // in DBL mode use all digits
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 -(number of displayed digits + 1)
# 002
- // -(the number of displayed digits + 3)
ENTER // here X = -(number of displayed + 3) (or -34)
10^x
STO .04 // save epsilon
SQRT // compute the convergence threshold, .05,
SDL 001 // = 10*sqrt(epsilon)
STO .05 // save the convergence threshold
DROP
ABS // compute the maximum level
LB
ROUNDI
# 002
+ // maxlevel = round(log2(digits)) + 2
STO .06 // save level
FC? .03 // the sinhsinh case and the expsinh case
SKIP 004 // when the finite limit == 0 can use a
# 001 // smaller epsilon
# 000
NEIGHB
STO .04 // save epsilon for such cases
// precision parameters (epsilon, thr, level) ready ****
// compute maximum t value (tm) ************************
// tm is computed in such a way that all abscissas are
// representable at the given precision and that the
// integrand is never evaluated at the limits
FC? .05 // tanhsinh mode?
SKIP 004 // no, check other modes
# 001 // start tm computation for tanh-sinh mode
RCLMIN .02
RCL+ X // X = 2*MIN(1, bma2)
SKIP 008 // goto remaining tm calculation
FS? .03 // neither SS nor ES with 0 limit?
SKIP 004 // no, must be SS or ES with 0 limit
# 1/2 // start tm computation for ES with != 0 limit
RCL/ .03 // RCP/ bpa2
ABS // X = ABS(1/(2*bpa2))
SKIP 002 // goto remaining tm calculation
RCL .04 // start tm computation for SS & ES0 cases
SQRT // X = SQRT(eps)
RCL/ .04 // continue tm computation for all cases
LN
FC? .05 // not in TS mode?
RCL+ X // 2*ln(...) for all cases except TS
RCL+ X // 4*ln(...) for all cases except TS (2*LN(...))
PI
/
LN
STO .07 // tm done
// maximum t (tm) ready ********************************
// level loop ******************************************
# 002
STO .08 // h = 2
LBL 20
# 000
STO .09 // ssp = 0
# 001
STO .10 // j = 1
# 002
STO/ .08 // h /= 2
RCL .08 // X = t
// j loop ++++++++++++++++++++++++++++++++++++++++++++++
// compute abscissas and weights ----------------------
LBL 30
COSH // cosh(t) (cosh is much faster than sinh/tanh)
STO .11 // save for later
x^2
DEC X
SQRT
PI
*
# 002
/ // pi/2*sqrt(cosh(t)^2 - 1) = pi/2*sinh(t)
FS? .02 // ES mode?
SKIP 002 // yes, want EXP
COSH // no, want COSH
SKIP 001
EXP
FILL
FC? .05 // TS mode?
SKIP 002 // no, done with weight (for ES, SS cases)
x^2 // end weight computation for TS mode
1/x
STO .12 // save weight
DROP // stack filled with exp/cosh of pi/2*sinh(t)
FS? .02 // ES mode?
SKIP 003 // yes, skip this part
x^2 // no ES mode, get r (the abscissa in the
DEC X // normalized domain)
SQRT
FS? .05 // TS mode?
RCL/ Y // yes, adjust r
FS? .04 // ES mode -infinity?
+/- // yes, adjust r
STO .13 // save normalized abscissa
// done with abscissas and weights --------------------
// evaluate integrand ----------------------------------
RCL* .02 // r*(b - a)/2
RCL+ .03 // (b + a)/2 + r*(b - a)/2
FILL
XEQa // f(bpa2 + bma2*r)
SPEC? // do not stop in error (flag D is set)
# 000
RCL* .12 // fplus*w
x<> .13 // p = fplus*w stored, r in X
RCL .03 // (b + a)/2
RCL .02 // (b - a)/2
FS? .02 // ES mode?
SKIP 003
RCL* Z // no ES mode, "normal" abscissa
-
SKIP 002
RCL/ Z // ES mode, "inverse" abscissa
+ // X = (b + a)/2 ± (b - a)/2*(r or 1/r)
FILL
XEQa // f(bpa2 ± bma2*(r or 1/r))
SPEC? // do not stop in error (flag D is set)
# 000
FS? .02 // ES mode?
SKIP 002
RCL* .12 // no ES mode, normal weight, fminus*w
SKIP 001
RCL/ .12 // ES mode, inverse weight, fminus/w
RCL+ .13 // p = fplus*w + fminus*(w or 1/w)
RCL* .11 // X = ch*p
STO+ .09 // ssp += p
// done with integrand --------------------------------
// level 0 is special ----------------------------------
FC? .06 // is level == 0?
GTO 40 // yes, go to next level
INC .10 // no, sweep odd abscissas, j += 1
ABS // X = ABS(p)
RCL .09 // ssp
RCL* .04 // eps
ABS // X = ABS(ssp*eps), Y = ABS(p)
x>=? Y // ABS(p) <= ABS(eps*eps)?
GTO 50 // then done with level
// and check user interrupt ............................
KEY? Y // any key?
SKIP 005 // no, continue
# 035 // yes, check keycode for <-
x!=? Z // key is not [<-]?
SKIP 002 // it is not [<-], continue
RCL .15 // it is, recall the last computation
GTO 70 // and exit
// if level > 0 done ----------------------------------
LBL 40
INC .10 // j += 1
RCL .10 // X = j
RCL* .08 // X = t = j*h
x<=? .07 // t <= tm?
GTO 30 // yes, continue j loop
// done with j loop ++++++++++++++++++++++++++++++++++++
LBL 50
RCL .09 // no, update ss. X = ssp
STO+ .14 // ss += ssp
FS? .06 // is level > 0 ----------------------------
GTO 60 // yes, progress to next level
# 001 // no, add 1st series term
FS? .04 // left?
+/-
RCL .03 // (b + a)/2
FS? .02 // ES mode?
+
FILL
XEQa // f(bpa2 (+/-1 in ES mode))
SPEC? // do not stop in error (flag D is set)
# 000
STO+ .14 // ss += f(bpa2)
// done with level 0 ----------------------------------
// apply constant coeffs to sum ------------------------
LBL 60
RCL .14
RCL* .02 // ss*bma2
RCL* .08 // ss*bma2*h
PI
* // ss*bma2*h*pi
# 002
/ // ss*bma2*h*pi/2
FS? .01 // reverse?
+/- // yes,so change sign
// done with constant coeffs --------------------------
// show progress --------------------------------------
TOP? // show value of approximation
MSG 25 // MSG_INTEGRATE
// check convergence ----------------------------------
RCL .09
RCL+ X // 2*ssp
RCL- .14 // 2*ssp - ss
ABS // ABS(2*ssp - ss)
RCL .05
RCL* .14 // thr*ss
ABS // ABS(thr*ss)
x>? Y // ABS(2*ssp - ss) < ABS(thr*ss)?
SKIP 004 // yes, computation done
z<> .15 // no, save ss (with all coeffs applied)
SF .06 // mark level 0 done,
DSL .06 // update level &...
GTO 20 // loop.
RCL Z // recall result
// done with level loop ********************************
// compute error estimation and check for bad results **
// here with result in X
LBL 70
RCL X // stack: ss-ss-?-?
ABS // stack: |ss|-ss-?-?
RCL Y // stack: ss-|ss|-ss-?
RCL- .15 // stack: (ss-ss1)-|ss|-ss-?
ABS // stack: err-|ss|-ss-?
RCL X // stack: err-err-|ss|-ss
SDL 001 // stack: 10*err-err-|ss|-ss
x<? Z // 10*err < |ss|?
SKIP 004 // yes, assume result is OK
x<> Z // no, bad result. stack: |ss|-err-10*err-ss
+ // stack: new_err-10*err-ss-ss
# 000 // stack: 0-new_err-10*err-ss
STO T // stack: 0-new_err-10*err-0
x<> T // stack: 0-err-... or ss-err-...
// done with bad results ******************************
// exit **********************************************
LBL 80
TOP? // show new computed value (0)
MSG 25 // MSG_INTEGRATE
cRCL .00 // stack: b-a-s-ss1
// [->]A..D // try to restore ABCD
STO L // b into L
cx<> Z // stack: s-ss1-b-a
FC? .00 // restore user's D flag (DF)
CF D
RTN // bye
// exit cause a == b **********************************
LBL 90
# 000 // error & result both 0
SKIP 002
// exit cause a or b is NaN ****************************
LBL 95
# NaN // error & result both NaN
RCL X
GTO 80 // exit
END
|
Search
Member List
Calendar
Help
/
(34S) Integration (using the double-exponential method)
