The Museum of HP Calculators

HP Forum Archive 20

[ Return to Index | Top of Index ]

Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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: Integration
Message #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.


[ Return to Index | Top of Index ]

Go back to the main exhibit hall