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

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 (x³ - 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.

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

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

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

Quote:

---

Assuming the current implementation, the problem will occur with any display setting except FIX, i.e. in ALL, SCI and ENG.

---

Franz

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

I tried the originally posted function on my trusted 35s with P₂(x) = (3x² - 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

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

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

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.

Hallo Dieter,

Thanks for investigating. I leave the case to Pauli.

d:-)

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

- Pauli

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

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

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

Pauli,

Quote:

---

ULP should work the same in xrom code as elsewhere.

---

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

             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

Dieter

Edited: 14 Oct 2013, 2:45 p.m.

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

Dieter

I wrote:

Quote:

---

... and after this the iteration terminates with a final (standard precision) result of 7,81753878289 E-17.

---

Dieter

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.

Quote:

---

Say hello to v. 3461

---

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

I've implemented this version.

- Pauli

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

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

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

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

   ...

Any objections or suggestions?

Dieter

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.

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

Maybe someone knows a good way to estimate the integral's uncertainty for the Romberg algorithm used here.

Dieter

STO- Y would give an uncertainty estimate of sorts.

- Pauli