WP34s integration trapped in infinite loop



#27

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


#28

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.

#29

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


#30

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


#31

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


#32

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


#33

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


#34

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<? 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 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.


#35

Hallo Dieter,

Thanks for investigating. I leave the case to Pauli.

d:-)

#36

I've add these changes but can't build new images at the moment. Assuming I got everything correct.


- Pauli


#37

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


#38

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


#39

I've uploaded a new build (3460) to SF with Pauli's latest version of the integrator for testing.


#40

Say hello to v. 3461 - see below. ;-)

Dieter


#41

Quote:
Say hello to v. 3461

If Pauli changes the code and I'll do the build, the final revision will be 3462.
#42

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<? 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.


#43

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.


#44

I've implemented this version.


- Pauli


#45

A new revision (3463) is built and online on SF.

#46

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.


#47

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


#48

STO- Y would give an uncertainty estimate of sorts.


- Pauli

#49

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

#50

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


#51

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


#52

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,947 12-05-2013, 02:32 AM
Last Post: Nick_S
  HP Prime integration Richard Berler 1 1,231 10-06-2013, 10:52 PM
Last Post: Helge Gabert
  New community-maintained version of "Calculators benchmark: add loop" Pier Aiello 20 6,410 09-12-2013, 02:42 AM
Last Post: Pier Aiello
  integration on 39gII emulator Wes Loewer 29 7,362 06-07-2013, 05:58 PM
Last Post: Chris Smith
  WP-34S Integration Richard Berler 15 3,887 03-08-2013, 02:29 AM
Last Post: Walter B
  HP 34S integration Richard Berler 16 4,315 02-18-2013, 04:42 PM
Last Post: Marcus von Cube, Germany
  Calculator Speed Benchmark (Add Loop) Thomas Chrapkiewicz 2 1,515 01-20-2013, 11:24 AM
Last Post: Thomas Chrapkiewicz
  [WP34S] WP34S firmware on the AT91SAM7L-STK dev kit? jerome ibanes 1 1,243 10-04-2012, 04:59 PM
Last Post: Paul Dale
  [WP34S] Speeding up the Romberg Integration Les Wright 14 4,247 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,809 05-17-2012, 11:58 PM
Last Post: Les Wright

Forum Jump: