WP34s integration trapped in infinite loop « Next Oldest | Next Newest »

 ▼ Bernd Grubert Unregistered Posts: 11 Threads: 3 Joined: May 2012 10-10-2013, 04:32 AM 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[sub-n] 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 ▼ Paul Dale Unregistered Posts: 3,229 Threads: 42 Joined: Jul 2006 10-10-2013, 05:37 AM 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.612851843479552E-16. The exact integral here is (x3 - 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 Gauss-Kronrod 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. Paul Dale Unregistered Posts: 3,229 Threads: 42 Joined: Jul 2006 10-10-2013, 05:40 AM 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 ▼ Bernd Grubert Unregistered Posts: 11 Threads: 3 Joined: May 2012 10-10-2013, 04:55 PM 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 ▼ Dieter Unregistered Posts: 653 Threads: 26 Joined: Aug 2010 10-10-2013, 05:19 PM 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 E-16) appears after a few seconds. Dieter ▼ fhub Unregistered Posts: 1,216 Threads: 75 Joined: Jun 2011 10-10-2013, 05:38 PM 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 ▼ Paul Dale Unregistered Posts: 3,229 Threads: 42 Joined: Jul 2006 10-10-2013, 06:18 PM 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 ▼ Dieter Unregistered Posts: 653 Threads: 26 Joined: Aug 2010 10-13-2013, 12:10 PM I tried the originally posted function on my trusted 35s with P2(x) = (3x2 - 1) / 2. Even in ALL mode it came back within a few seconds with a result of 8,16 E-13, i.e. approx. zero for a 12-digit 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 3E-14 and 6E-15 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 (3E-16, 4E-17, 1E-16, ...) 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 0-10, 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] Z TOP? MSG MSG_INTEGRATE RTN You should check the way ULP works in XROM code - it should return 1E-15 in SP and 1E-33 in DP mode. In ALL 03 mode, the integral of the originally posted function is approximated this way: 0,0634765625 -2,0751953125 E-4 -9,15527343744 E-6 -8,021902 E-17 ... and after this the iteration terminates with a final (standard precision) result of 7,81753878289 E-17. 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 E-19 5,73024595021 E-19 5,72957175023 E-19 5,72957800022 E-19 5,72957795128 E-19 5,72957795131 E-19 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. ▼ Walter B Unregistered Posts: 4,587 Threads: 105 Joined: Jul 2005 10-13-2013, 07:13 PM Hallo Dieter, Thanks for investigating. I leave the case to Pauli. d:-) Paul Dale Unregistered Posts: 3,229 Threads: 42 Joined: Jul 2006 10-13-2013, 08:01 PM I've add these changes but can't build new images at the moment. Assuming I got everything correct. - Pauli ▼ Dieter Unregistered Posts: 653 Threads: 26 Joined: Aug 2010 10-14-2013, 08:13 AM 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 1E-15. A better solution would get rid of the x[approx] test completely. Dieter ▼ Paul Dale Unregistered Posts: 3,229 Threads: 42 Joined: Jul 2006 10-14-2013, 08:46 AM 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 ▼ Marcus von Cube, Germany Unregistered Posts: 3,283 Threads: 104 Joined: Jul 2005 10-14-2013, 11:40 AM I've uploaded a new build (3460) to SF with Pauli's latest version of the integrator for testing. ▼ Dieter Unregistered Posts: 653 Threads: 26 Joined: Aug 2010 10-14-2013, 02:45 PM Say hello to v. 3461 - see below. ;-) Dieter ▼ Marcus von Cube, Germany Unregistered Posts: 3,283 Threads: 104 Joined: Jul 2005 10-14-2013, 03:32 PM Quote: Say hello to v. 3461 If Pauli changes the code and I'll do the build, the final revision will be 3462. Dieter Unregistered Posts: 653 Threads: 26 Joined: Aug 2010 10-14-2013, 02:43 PM 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 1E-33. 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] Z TOP? MSG MSG_INTEGRATE RTN Time for a new build. :-) Dieter Edited: 14 Oct 2013, 2:45 p.m. ▼ Dieter Unregistered Posts: 653 Threads: 26 Joined: Aug 2010 10-14-2013, 03:03 PM 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 close-to-zero-results 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. ▼ Paul Dale Unregistered Posts: 3,229 Threads: 42 Joined: Jul 2006 10-14-2013, 05:35 PM I've implemented this version. - Pauli ▼ Marcus von Cube, Germany Unregistered Posts: 3,283 Threads: 104 Joined: Jul 2005 10-15-2013, 01:28 PM A new revision (3463) is built and online on SF. David Maier Unregistered Posts: 10 Threads: 3 Joined: Feb 2013 10-18-2013, 10:47 AM 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. ▼ Dieter Unregistered Posts: 653 Threads: 26 Joined: Aug 2010 10-18-2013, 03:05 PM 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 ▼ Paul Dale Unregistered Posts: 3,229 Threads: 42 Joined: Jul 2006 10-18-2013, 08:28 PM STO- Y would give an uncertainty estimate of sorts. - Pauli Paul Dale Unregistered Posts: 3,229 Threads: 42 Joined: Jul 2006 10-14-2013, 05:27 PM Quote:I thought XROM routines always run in DP mode, so that ULP(1) will always return 1E-33. 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 Dieter Unregistered Posts: 653 Threads: 26 Joined: Aug 2010 10-14-2013, 02:50 PM I wrote: Quote: ... and after this the iteration terminates with a final (standard precision) result of 7,81753878289 E-17. This is caused by an error that is corrected in the bugfix I posted below. With this correction, the last approximation -8,021902 E-17 is returned in X (and the second last -9,155 E-6 in Y). Dieter ▼ Miguel Toro Unregistered Posts: 239 Threads: 55 Joined: Sep 2006 10-16-2013, 09:31 AM 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 ▼ Dieter Unregistered Posts: 653 Threads: 26 Joined: Aug 2010 10-17-2013, 08:50 AM 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 (|b-a| << 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 x-value. 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

 Possibly Related Threads… Thread Author Replies Views Last Post Integration question and "RPN" mode comment Craig Thomas 16 5,843 12-05-2013, 02:32 AM Last Post: Nick_S HP Prime integration Richard Berler 1 1,205 10-06-2013, 10:52 PM Last Post: Helge Gabert New community-maintained version of "Calculators benchmark: add loop" Pier Aiello 20 6,214 09-12-2013, 02:42 AM Last Post: Pier Aiello integration on 39gII emulator Wes Loewer 29 7,089 06-07-2013, 05:58 PM Last Post: Chris Smith WP-34S Integration Richard Berler 15 3,757 03-08-2013, 02:29 AM Last Post: Walter B HP 34S integration Richard Berler 16 4,183 02-18-2013, 04:42 PM Last Post: Marcus von Cube, Germany Calculator Speed Benchmark (Add Loop) Thomas Chrapkiewicz 2 1,447 01-20-2013, 11:24 AM Last Post: Thomas Chrapkiewicz [WP34S] WP34S firmware on the AT91SAM7L-STK dev kit? jerome ibanes 1 1,215 10-04-2012, 04:59 PM Last Post: Paul Dale [WP34S] Speeding up the Romberg Integration Les Wright 14 4,066 05-31-2012, 03:39 PM Last Post: Marcus von Cube, Germany HP-IL: Using 82162A Printer and PIL-Box in Same Loop Les Wright 12 3,695 05-17-2012, 11:58 PM Last Post: Les Wright

Forum Jump: