(34S) Integration (using the double-exponential method)

[emece67]

(03-26-2017 12:25 AM)Paul Dale Wrote:  The built in integrator doesn't corrupt the stack. It leaves the integral in X, the error estimate in Y and the original limits of integration in Z & T and also preserves the original X value in L. Refer to the 34C manual page 211. I'm don't remember what it does to A, B, C & D with stack size 8 -- XROM switches to four level stack mode on entry so they won't be directly altered, although the fills before calling user code might.

Surely that's the behavior of the 34C, but the 34s:

1. Returns the computed result in X
2. Returns in Y what seems to be some kind of a previous approximation to the result (it is not much smaller than the result, as expected for an error estimation, but approximately equal to it)
3. Returns the integration limits in Z & T
4. Returns the upper integration limit in L
5. Corrupts ABCD with some unknown data (no matter if the user is in SSIZE4 or SSIZE8, ABCD are always corrupted)
   

Thus, I've modified my program to behave more or less in the same way, so it now (v1.1r and v1.1x):

1. Returns the computed result in X (as before)
2. Returns in Y the previous approximation to the result
3. Returns the integration limits in Z & T
4. Returns the upper integration limit in L (as before)
5. Keeps ABCD untouched (no matter if the user is in SSIZE4 or SSIZE8)
   

(03-26-2017 11:36 AM)Dieter Wrote:  There is a logic that calculates an epsilon value depending on the display setting. This seems to work fine in FIX mode. But how does this work in SCI/ENG or ALL mode? Note that (unlike RDP) ROUND rounds 1/9 to at most 11 places. Or a certain number of significant digits in SCI/ENG mode.

Later, epsilon (eps) is mainly used to compute some other values (thr & tm). Such values depend on the number of significant digits displayed thus, epsilon is not an absolute threshold for anything, but a number carrying information about the number of significant digits displayed. This way, eps is more meaningful here in ALL/SCI/ENG modes than in FIX mode.

There was also an absolute comparison with eps that allows the program not to iterate (mainly when both integration ends are infinite) through a big number of abscissas whose weights are so small than they do not contribute to the final result. This was only an heuristic to speed-up computation, not the main convergence test. In any case, thanks to your comment, I've found that such test interferes with cases where the integrand is always very small, so I have changed such test and now all tests are relative.

(03-26-2017 11:36 AM)Dieter Wrote:  At LBL 40 the program checks if the integration variable is within the integration limits. This is done by calculating the current x and comparing it to tm in R.06. In such cases I always wonder if roundoff errors might not lead to an error here. Suppose tm is 2 and h*j rounds to 2,000000000000001 so that the test returns false. Can this possibly happen? A simple counter would be a safe solution.

j is a small integer only updated with INC. h is an integer power of 2 from 2 to 1/256 updated by continued division by 2, all these powers have an exact decimal representation in the wp34s. Thus, h*j is always exactly representable in the machine, so no concerns about roundoff here.

(03-26-2017 11:36 AM)Dieter Wrote:  Very small values of the integral, e.g. 1E–16, may be calculated because the actual value is zero, but just as well the correct value of the true integral. Is the program able to distinguish these two cases?

Being now all tests relative, integrands that are always small do not pose any special problem. But the program does not distinguish if the true integral is zero and the computed, small, result is due to round-off, or the computed, small, result is a nice approximation of a small but not zero integral.

(03-26-2017 11:36 AM)Dieter Wrote:  BTW, Pi/2 is directly available as CNST 86. ;-)

In my local builds that's not true, cause I have a modified CONST menu, so I try to avoid such constructs.