The Museum of HP Calculators

HP Forum Archive 20

 IntegrationMessage #1 Posted by Marcus von Cube, Germany on 19 Jan 2012, 12:49 p.m. I've stumbled across a problem with the Romberg integrator in WP 34S. Just as a test if I had made a mess of the firmware again - quite a common situation in the last few weeks - I started with the integral from 0 to 1 of 1/x. This is infinite but the poor calculator was busily trying to approximate the result until the algorithm failed with a programming error (register limit exceeded). We tackled the problem by limiting the number of iterations and now the device gives up after a few minutes (a few seconds on the emulator). I then tried the same with my 15C LE in FIX 4 mode: After 15 minutes or more the display still showed "running" and I stopped it to save on the batteries. My question is, how do other incarnations of the integrate function fare on this task?

 Re: IntegrationMessage #2 Posted by Norman Dziedzic on 19 Jan 2012, 1:22 p.m.,in response to message #1 by Marcus von Cube, Germany 50g seems to hang although I only waited about a minute on 'Integral(0,0,INV(X),X)' From inside equation writer, entered the integral, highlighted the whole thing and got + Infinity very quickly.

 Re: IntegrationMessage #3 Posted by Marcus von Cube, Germany on 19 Jan 2012, 1:42 p.m.,in response to message #2 by Norman Dziedzic Quote: Integral(0,0,INV(X),X) I assume you mean Integral(0,1,INV(X),X). The immediate return with +infinity may be a result of the CAS knowing that the indefinite integral of 1/x is lnx which evaluates to infinity on the lower bound.

 Re: IntegrationMessage #4 Posted by Howard Owen on 19 Jan 2012, 8:38 p.m.,in response to message #3 by Marcus von Cube, Germany If you are gong to special case the integral for the purpose of limiting the iterations, why not do what the 50G does and shortcut the processing up front? Am I misunderstanding what's happening in your code?

 Re: IntegrationMessage #5 Posted by Paul Dale on 19 Jan 2012, 8:39 p.m.,in response to message #4 by Howard Owen We don't and won't do symbolic integration. It looks like the 50G knows that the integral of 1/x dx is Ln(x) and it is using that to short cut the evaluation. - Pauli

 Re: IntegrationMessage #6 Posted by Kiyoshi Akima on 19 Jan 2012, 1:57 p.m.,in response to message #2 by Norman Dziedzic My 50g hits the limit of 65535 function evaluations and returns a result of 22.8202972649 with IERR of -2.22365339103E-10 . I don't have my 15C+ with me, but I believe it will hit its limit at 16383 function evaluations.

 Re: IntegrationMessage #7 Posted by Marcus von Cube, Germany on 19 Jan 2012, 2:17 p.m.,in response to message #6 by Kiyoshi Akima From a recent mail from Pauli: Quote: The integrator does 8191 function evaluations maximum now.

 Re: IntegrationMessage #8 Posted by Paul Dale on 19 Jan 2012, 4:20 p.m.,in response to message #7 by Marcus von Cube, Germany This is more a bug on my part -- it does 8191 evaluations but only the first 4096 are actually used. The "we've been going too long and need to stop test" is done too late. We can set this iteration limit to any power of two we want and I will look into the earlier termination. - Pauli

 Re: IntegrationMessage #9 Posted by Marcus von Cube, Germany on 20 Jan 2012, 11:40 a.m.,in response to message #8 by Paul Dale I can't see the bug Pauli mentioned but I have increased the limit to 16383 again which is good for just below 5 minutes of execution time for 1/x (timed on my crystal-less 20b with TICKS). When the integrator is called from the keyboard it will show a progress message. In a prior version I had added an exit similar to SLV ("Solution not Found") and a skip next instruction return. It's back to how it was. Convergence can be easily tested with an x[approx]? Y test in a program or by manually comparing the values in x and y. Edited for typos. Edited: 20 Jan 2012, 12:00 p.m.

 Re: IntegrationMessage #10 Posted by peacecalc on 19 Jan 2012, 4:44 p.m.,in response to message #6 by Kiyoshi Akima Hello folks, my hp 50g shows an astonishing behavier: If I start in exact mode with: for example 4: 192E-14 3: 1 2: '1/X' 1: X Start with the integration symbol (RS Integral symbol), I get the question: approx. mode on? I answer: yes and I get immediatly the integration value: 2,6979E1 If I start within the approximation mode with the "same" (of course they are not interpretated as exact values) input numbers, I have to wait for the result some minutes. This looks like a good work around against a long integration duration: just start in exact mode. Sincerely peacecalc Edited: 19 Jan 2012, 4:45 p.m.

 Re: IntegrationMessage #11 Posted by peacecalc on 19 Jan 2012, 5:07 p.m.,in response to message #6 by Kiyoshi Akima Hello Kiyoshi Akima, where on the calc hp 50g you can see how many times the function is calculated for the integration? The "IERR" variable I've already found. Sincerely peacecalc

 Re: IntegrationMessage #12 Posted by Kiyoshi Akima on 19 Jan 2012, 5:17 p.m.,in response to message #11 by peacecalc Quote: where on the calc hp 50g you can see how many times the function is calculated for the integration? The "IERR" variable I've already found. Rather than integrating the program ```<< INV >> ``` I merely integrate the program ```<< 'C' INCR DROP INV >> ``` with the global variable 'C' set to some known initial value. The same thing can be done on RPN machines by including something along the lines of ISG 0 in the function, followed by a NOP if necessary. For this particular problem, another possible coding would be ```LBL A 1 STO+ 0 x<>y / RTN ```

 Re: IntegrationMessage #13 Posted by peacecalc on 20 Jan 2012, 8:29 a.m.,in response to message #12 by Kiyoshi Akima Hello Kiyoshi Akima, thanx for your fast reply. The little proggie only works in approxmation mode. I have a theory why the calculation works better in exact mode. The calc find the antiderivate, but cannot express it with the boundaries, so the calc forced the user to switch to approximation mode. That is shurly faster then the pure approximation. Sincerely peacecalc Edited: 20 Jan 2012, 8:29 a.m.

 Re: IntegrationMessage #14 Posted by Kiyoshi Akima on 20 Jan 2012, 12:54 p.m.,in response to message #6 by Kiyoshi Akima The 15C+ takes about an hour and a half to hit its limit of 65535 evaluations.

 Re: IntegrationMessage #15 Posted by Oliver Unter Ecker on 19 Jan 2012, 1:31 p.m.,in response to message #1 by Marcus von Cube, Germany ND1: << 0 1 '1/x' 'x' Integral >> => bad argument [no wait time] W|A: integrate 1/x from 0 to 1 => (integral does not converge) Edited: 19 Jan 2012, 1:33 p.m.

 Re: IntegrationMessage #16 Posted by Marcus von Cube, Germany on 19 Jan 2012, 1:43 p.m.,in response to message #15 by Oliver Unter Ecker Has ND1 a CAS of some sort?

 Re: IntegrationMessage #17 Posted by Oliver Unter Ecker on 19 Jan 2012, 2:51 p.m.,in response to message #16 by Marcus von Cube, Germany The currently released version has no CAS. ND1 (and ND0, the free version) can do numeric integration, deal with algebraic expressions (by themselves or in vectors/matrices), has PEVAL, can use OBJ-> to "unpeel" expressions, and can do some things like prime factorization, BigInts, and other stuff listed by the 50g AUR in CAS section. The last couple of months, I've been working on integrating with a commercial CAS, which will become available within 1-2 months. (Beta testers wanted.) Details are being finalized, but I'll be very excited when I can make the announcement. It's not just any CAS...

 Re: IntegrationMessage #18 Posted by Marcus von Cube, Germany on 19 Jan 2012, 3:00 p.m.,in response to message #17 by Oliver Unter Ecker How does ND1 then find out about the integral being infinite? Or is this just a failure because it attempts to evaluate the function at its endpoints?

 Re: IntegrationMessage #19 Posted by Oliver Unter Ecker on 19 Jan 2012, 4:51 p.m.,in response to message #18 by Marcus von Cube, Germany Sorry, I didn't pick up on your question. To avoid singularities, the endpoints are never evaluated. In this specific case, an overflow in the interior of the interval is causing NaN values, which leads to the reported error (which should be changed to something more appropriate).

 Re: IntegrationMessage #20 Posted by Gerson W. Barbosa on 19 Jan 2012, 2:26 p.m.,in response to message #1 by Marcus von Cube, Germany HP-42S: It was still 'Integrating' when I stopped it after thirty minutes. (ACC=0.0001) Free42 Decimal v1.4.62: 26.97737115553 (ACC<=0.0001). (This corresponds to a lower limit of about 1.92e-12).

 Re: IntegrationMessage #21 Posted by Valentin Albillo on 19 Jan 2012, 3:02 p.m.,in response to message #1 by Marcus von Cube, Germany Hi, Marcus: For the HP-71B, you get in a short while: ```> INTEGRAL(0,1,0,1/IVAR) 22.8202902142 > IBOUND -2.22365339102E-11 ``` where the negative IBOUND means the specified precision (1E-12) wasn't achieved. Conversely: ```> INTEGRAL(1,2,0,1/IVAR) .69314718056 > LN(2) .69314718056 >IBOUND 6.93146313016E-13 ``` in a few seconds (positive IBOUND, so 12 digits were indeed achieved). Best regards from V.

 Re: IntegrationMessage #22 Posted by Paul Dale on 19 Jan 2012, 4:23 p.m.,in response to message #21 by Valentin Albillo I think we can work out the maximum number of iterations each integrator uses based on the result of this integral. At least the results being mentioned look rather close to some of those I got with larger iteration limits. - Pauli

 Re: IntegrationMessage #23 Posted by Valentin Albillo on 19 Jan 2012, 6:26 p.m.,in response to message #22 by Paul Dale Quote: I think we can work out the maximum number of iterations each integrator uses based on the result of this integral. It's well known that the modified iterative Romberg algorithm as implemented in the HP-71B does a maximum of 65,535 function evaluations, which can be simply demonstrated like this: ```>LIST 10 DEF FNF(X) @ N=N+1 @ FNF=1/X @ END DEF 20 N=0 @ DISP INTEGRAL(0,1,0,FNF(IVAR)),N >RUN 22.8202902142 65535 ``` Best regards from V. ``` ```

 Re: IntegrationMessage #24 Posted by John B. Smitherman on 19 Jan 2012, 6:07 p.m.,in response to message #1 by Marcus von Cube, Germany Marcus, my TI83+SE gives the error message "TOL NOT MET" within 10 seconds. Regards, John

 Re: IntegrationMessage #25 Posted by From Hong Kong on 19 Jan 2012, 11:22 p.m.,in response to message #1 by Marcus von Cube, Germany TI-NSpire CX CAS gave infinity immediately.

 Re: IntegrationMessage #26 Posted by Marcus von Cube, Germany on 20 Jan 2012, 3:23 a.m.,in response to message #25 by From Hong Kong Can you try in approximate mode again?

 Re: IntegrationMessage #27 Posted by From Hong Kong on 20 Jan 2012, 4:43 a.m.,in response to message #26 by Marcus von Cube, Germany I'll try and get back to you soon. TI-85 and TI-86 gave "ERROR 33 TOL NOT MET" after about 40 seconds. FX-4500P returned "Ma ERROR" immediately. Edited: 20 Jan 2012, 5:00 a.m.

 Re: IntegrationMessage #28 Posted by From Hong Kong on 20 Jan 2012, 5:08 a.m.,in response to message #26 by Marcus von Cube, Germany It gave the same result immediately. Edited: 20 Jan 2012, 5:08 a.m.

Go back to the main exhibit hall