▼
Posts: 11
Threads: 3
Joined: May 2012
Hello,
Sorry if that has been adressed before.
I tried to integrate the 2nd Lengendre polynom over 1 ... 1, which should give a result close to zero. With "ALL 00" the integration seems to be trapped in an infinite loop. The problem doesn't occur with "FIX 11".
Here is the function, that I'm integrating:
0001 ****LBL'LE2'
0002 **LBL A
0003 2
0004 x[<>] Y
0005 P[subn]
0006 END
The integration is done the following way:
1
ENTER[^]
+/
[integral]'LE2'
I'm using Version 3458 with printer support.
BTW: The WP 34S is a really fine machine. Thanks to all who made it possible!
Regards
Bernd
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
It isn't in an infinite loop. The Romberg integration is having difficulties and it will call your user subroutine 16383 times before returning the answer 6.612851843479552E16. The exact integral here is (x^{3}  x) / 2 + C, when taken over the interval specified gives zero.
The integration routine has difficulties with integrals that evaluate to zero when in ALL mode. The problem being related to the approximately equal function which rounds to the current display mode  with zero, the rounding doesn't do much of anything except in FIX mode.
Interestingly, version 2 firmware, which uses a fixed point GaussKronrod quadrature, will get the correct result very quickly. This quadrature is known exact for polynomials of such low degree.
 Pauli
Edited: 10 Oct 2013, 5:38 a.m.
Posts: 3,229
Threads: 42
Joined: Jul 2006
A work around would be to break the interval of integration somewhere that doesn't produce a zero result. E.g. integrate 1 .. 0.5 and 0.5 .. 1 and add the results together. Unfortunately, 1 .. 0 and 0 .. 1 both run into the same problem.
 Pauli
▼
Posts: 11
Threads: 3
Joined: May 2012
Thanks for the answer.
I understand that that 16383 is not infinite. But it is quite close  at least with the real calc running on battries :)
At first I thought that this problem has to do with the ALL setting, but I noticed that this problem occurs also with SCI 10. So maybe for me the best solution is to set the accuracy not too high when calculating integrals.
Regards
Bernd
▼
Posts: 653
Threads: 26
Joined: Aug 2010
Assuming the current implementation, the problem will occur with any display setting except FIX, i.e. in ALL, SCI and ENG.
Try FIX 11 and the result (approx. 1,07 E16) appears after a few seconds.
Dieter
▼
Posts: 1,216
Threads: 75
Joined: Jun 2011
Quote:
Assuming the current implementation, the problem will occur with any display setting except FIX, i.e. in ALL, SCI and ENG.
Well, then the current implementation is not very clever  needing >16000 steps for a Romberg integration is definitely unacceptable.
Franz
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
Anyone is most welcome to submit a replacement integration routine. The source code is in trunk/xrom/integrate.wp34s and is no more than a standard keystroke program.
I've implemented two different integrators so far and don't really want to do a third.
 Pauli
▼
Posts: 653
Threads: 26
Joined: Aug 2010
I tried the originally posted function on my trusted 35s with P_{2}(x) = (3x^{2}  1) / 2. Even in ALL mode it came back within a few seconds with a result of 8,16 E13, i.e. approx. zero for a 12digit machine like this one. HP's integration algorithm uses a modified Romberg method, so there must be a way to overcome the limitations of the current 34s implementation.
The crucial point here is the "approx. equal" test in the current 34s Romberg code. In FIX n mode, this simply tests if the two last approximations agree if rounded to n decimals. This also works for values close to zero  for instance 3E14 and 6E15 are both rounded to zero, the two values agree, et voilà: the iteration terminates.
On the other hand, ALL, SCI and ENG will round the two last approximations to a certain number of significant digits. So if the iteration oscillates through several values very close to zero (3E16, 4E17, 1E16, ...) they will almost never agree in some or even all significant digits, and the iteration will not terminate (at least not before 16000+ function calls).
I examined the current 34s integration code, transferred it to the emulator (with a few adjustments, e.g. it uses the regular registers 010, and 11ff for the Romberg matrix), and finally I changed the exit condition in the following way:
After the call of the user function and the summation via STO+ r_Sk, also the absolute value of this increment is accumulated in a new register r_SkAbs, which of course has been set to zero at the beginning. This way, SkAbs times deltaU represents a kind of average function value over the interval. If now the final x[approx]y test fails, but the last approximation is very small compared to this average, the iteration can terminate:
x[approx]? Y
JMP int_done
RCL r_SkAbs
RCL[times] r_deltau
/
ABS
Num 1
ULP
x<? Y
JMP int_next_size
/* else exit */
/* Two matches in a row is goodness. */
int_done::
[cmplx]RCL cr_limits
STO L
[cmplx]x[<>] Z
TOP?
MSG MSG_INTEGRATE
RTN
You should check the way ULP works in XROM code  it should return 1E15 in SP and 1E33 in DP mode.
In ALL 03 mode, the integral of the originally posted function is approximated this way:
0,0634765625
2,0751953125 E4
9,15527343744 E6
8,021902 E17
... and after this the iteration terminates with a final (standard precision) result of 7,81753878289 E17.
On the other hand, an integral that really is that small is evaluated correctly without premature break. For instance, the integral of sin(x) * 10^{20} for x = 0..90 degrees (!) returns the following approximation sequence:
5,68490384437 E19
5,73024595021 E19
5,72957175023 E19
5,72957800022 E19
5,72957795128 E19
5,72957795131 E19
This matches the true result of 10 ^{20} * 180/Pi.
Might this be a feasible (and correct) solution?
Dieter
Addendum:
The x[approx] test compares two rounded values. It seems that, just like a manual ROUND command, this way at most 12 digits get compared. So the result will usually not carry 16 valid digits in SP mode, let alone 34 digits in DP mode: everything beyond 12 digits is not guaranteed. IMHO this is one more reason to get rid of the x[approx] test in the 34s integration routine.
Edited: 13 Oct 2013, 6:10 p.m.
▼
Posts: 4,587
Threads: 105
Joined: Jul 2005
Hallo Dieter,
Thanks for investigating. I leave the case to Pauli.
d:)
Posts: 3,229
Threads: 42
Joined: Jul 2006
I've add these changes but can't build new images at the moment. Assuming I got everything correct.
 Pauli
▼
Posts: 653
Threads: 26
Joined: Aug 2010
There still is the question of how the ULP function works in SP resp. DP mode for XROM code. Since the x[approx] test disregards any digits beyond the 12th anyway, I would suggest replacing the 1 ULP sequence by a simple 1 SDR 15, i.e. a constant 1E15.
A better solution would get rid of the x[approx] test completely.
Dieter
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
ULP should work the same in xrom code as elsewhere. The only thing xrom changes when activated is to force a four level stack. The xIN command does switch to double precision and an eight level stack, but the integration code doesn't call this function.
The integration code uses the approx test so that it has some kind of compatibility with HP's earlier integration code  the screen settings changed the result. We could lose this compatibility and do something different I guess. The approx test only fails when the integral is convergent to zero and the code you've already supplied seems to handle this case.
 Pauli
▼
Posts: 3,283
Threads: 104
Joined: Jul 2005
I've uploaded a new build (3460) to SF with Pauli's latest version of the integrator for testing.
▼
Posts: 653
Threads: 26
Joined: Aug 2010
Say hello to v. 3461  see below. ;)
Dieter
▼
Posts: 3,283
Threads: 104
Joined: Jul 2005
Quote:
Say hello to v. 3461
If Pauli changes the code and I'll do the build, the final revision will be 3462.
Posts: 653
Threads: 26
Joined: Aug 2010
Pauli,
Quote:
ULP should work the same in xrom code as elsewhere.
I thought XROM routines always run in DP mode, so that ULP(1) will always return 1E33. But anyway  here is a bugfix.
There still is an error that returns the wrong values on the stack if the integral is close to zero. Also, since the x[approx] test ignores everything beyond the first 12 digits anyway, I now prefer this version:
x[approx]? Y
JMP int_done
ENTER
RCL/ r_SkAbs
RCL/ r_deltau
ABS
SDL 15
x>1?
JMP int_next_size
DROP
/* else exit */
/* Two matches in a row is goodness. */
int_done::
[cmplx]RCL cr_limits
STO L
[cmplx]x[<>] Z
TOP?
MSG MSG_INTEGRATE
RTN
If you want to save 15 powers of 10 in the result, here is an alternative code:
x[approx]? Y
JMP int_done
ENTER
RCL/ r_SkAbs
RCL/ r_deltau
ABS
Num 1
SDR 15
x<? Y
JMP int_next_size
[cmplx]DROP
/* else exit */
/* Two matches in a row is goodness. */
int_done::
[cmplx]RCL cr_limits
STO L
[cmplx]x[<>] Z
TOP?
MSG MSG_INTEGRATE
RTN
Time for a new build. :)
Dieter
Edited: 14 Oct 2013, 2:45 p.m.
▼
Posts: 653
Threads: 26
Joined: Aug 2010
An old German proverb says: "Das einzig Beständige ist der Wandel". So here's a more conservative version that uses an additional (initially cleared) flag that is set for a final iteration. In other words: if the algorithm has detected that the integral is very close to zero, it runs one final iteration  just to be sure. This way the last two closetozeroresults are returned in X and Y.
FC?C flag_final /* "final" flag set */
x[approx]? Y /* OR x~y ? */
JMP int_done
RCL/ r_SkAbs
RCL/ r_deltau
ABS
Num 1
SDR 15
x>? Y /* integral close to zero? */
SF flag_final /* then do one final iteration */
JMP int_next_size
/* Two matches in a row is goodness. */
int_done::
[cmplx]RCL cr_limits
STO L
[cmplx]x[<>] Z
TOP?
MSG MSG_INTEGRATE
RTN
Dieter
Edited: 14 Oct 2013, 3:09 p.m.
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
I've implemented this version.
 Pauli
▼
Posts: 3,283
Threads: 104
Joined: Jul 2005
A new revision (3463) is built and online on SF.
Posts: 10
Threads: 3
Joined: Feb 2013
Just thought I would point out that, according to the printed version of the WP 34S Owner's Manual, integration returns "an upper limit of the uncertainty in Y". There's also a footnote that mentions the 34C and 15C Owner's Handbooks. The 34C and 15C do return a calculated uncertainty in Y, and their handbooks include a discussion of how this is calculated.
▼
Posts: 653
Threads: 26
Joined: Aug 2010
That's right  the 34C and 15C (and newer machines like the 35s, for that matter) return an estimate for the uncertainty in Y. However, the 34s does not. Neither before nor after the last modification of the integration routine. Actually the exit code remained unchanged. On exit, the stack contains these values:
T lower limit a
Z upper limit b
Y previous approximation
X result (final approximation)
L upper limit b
Yes, the 34s manual contains various references to other HP manuals, and there are cases where the 34s behaves differently compared to earlier HPs. This is one example, and it should be corrected in the manual. However, my PDF manual does not seem to contain the quoted "upper limit of the uncertainty in Y".
Maybe someone knows a good way to estimate the integral's uncertainty for the Romberg algorithm used here.
Dieter
▼
Posts: 3,229
Threads: 42
Joined: Jul 2006
STO Y would give an uncertainty estimate of sorts.
 Pauli
Posts: 3,229
Threads: 42
Joined: Jul 2006
Quote: I thought XROM routines always run in DP mode, so that ULP(1) will always return 1E33.
It is the xIN that starts most XROM commands that switches to DP mode. Before that it stays in whatever mode the device is currently in. Solve, integrate, product, sum and the derivatives don't call xIN anywhere so they stay single precision and don't gain the benefits of automatic argument processing and Last X handling.
 Pauli
Posts: 653
Threads: 26
Joined: Aug 2010
I wrote:
Quote:
... and after this the iteration terminates with a final (standard precision) result of 7,81753878289 E17.
This is caused by an error that is corrected in the bugfix I posted below. With this correction, the last approximation 8,021902 E17 is returned in X (and the second last 9,155 E6 in Y).
Dieter
▼
Posts: 239
Threads: 55
Joined: Sep 2006
Hi,
This is what I like of this project and why I continue to enjoy working (and playing, of course) with my trusted WP34s : the constant and really fast improvement of an already marvelous little machine.
Thank you Paul, Walter, Marcus and Dieter
Regards,
Miguel
▼
Posts: 653
Threads: 26
Joined: Aug 2010
Thank you very much, Miguel.
We should now discuss whether the new version includes a valid and realiable exit condition. In the updated algorithm this is done with an additional approximation for the integral of f(x):
Integral f(x) ~= 1,5 * r_ba4 * r_SkAbs * r_deltau
The basic idea now is simple: compare this integral of f(x) with the integral of f(x). If the latter is, say, 15 magnitudes smaller than the former, the algorithm assumes that the true value of the integral is zero, so it exits after one final iteration.
Since we compare just magnitudes, the factor 1,5 in the formula is omitted in the program. Also, it does not divide by r_ba4, which may cause problems if the interval is very small or very large (ba << 1 or >> 1). So I think this should be added:
...
FC?C flag_final /* "final" flag set */
x[approx]? Y /* OR x~y ? */
JMP int_done
RCL/ r_SkAbs
RCL/ r_deltau
RCL/ r_ba4
ABS
...
Finally, there still is one case where the program may exit with an error. This would happen if the sum SkAbs of all sampled (f(x) is zero, i.e. f(x) is zero for each and every sampled xvalue. This results in a division by zero. On the other hand, I think that in this case the two initial approximations of the integral would both be zero as well so that the algorithm exits with a result of zero. At least I tried this with f(x) = 0 and it worked. ;)
Any objections or suggestions?
Dieter
