HP 42S, DM 42 Integral

[lrdheat]

On my 42S, the integral of SIN X from 0 to 2pi produces 3.674E-12 (very close to the actual answer of 0, using an accuracy of 1E-5) nearly immediately.

My DM 42 takes a long enough time where I give up. Same 1E-5 accuracy, even tried different accuracies ranging from .1 to 1E-12 on the DM 42.

My code is lbl “fx”
MVAR X
RCL X
SIN X
END

Both the 42S and DM 42 very quickly come up with 2.000

Any idea as to why the DM 42 has trouble with the domain of (0,2pi) ?

[lrdheat]

Meant to say both quickly come up with 2.000 for the integral from 0 to pi.

[Thomas Okken]

The fact that the actual integral is zero is causing the problem. INTEG waits for the relative error to drop below ACC, but when the true result is zero, the approximated result keeps getting smaller at the same rate as the absolute error estimate, and so the relative error doesn't get better. After a certain number of iterations, it just gives up, but since each iteration takes twice as long as the last one, the total number of times the function has to be evaluated gets large rather quickly.

In Free42, you won't notice this because it usually runs on hardware that is so fast that even evaluating the function tens of thousands of times only takes a moment, but on the DM42 this does become noticeable. Clearly the HP-42S has a more aggressive termination condition for cases like this.

[Albert Chan]

(02-12-2021 10:58 PM)Thomas Okken Wrote:  The fact that the actual integral is zero is causing the problem. INTEG waits for the relative error to drop below ACC, but when the true result is zero, the approximated result keeps getting smaller at the same rate as the absolute error estimate, and so the relative error doesn't get better.

More likely, integral results are bouncing around 0, thus unable to detect convergence.
In other words, all it see is noise.

[Thomas Okken]

(02-12-2021 11:32 PM)Albert Chan Wrote:  

(02-12-2021 10:58 PM)Thomas Okken Wrote:  The fact that the actual integral is zero is causing the problem. INTEG waits for the relative error to drop below ACC, but when the true result is zero, the approximated result keeps getting smaller at the same rate as the absolute error estimate, and so the relative error doesn't get better.

More likely, integral results are bouncing around 0, thus unable to detect convergence.
In other words, all it see is noise.

You're right. Adding print statements in INTEG showing the result and the error estimate after each iteration produces this:

Code:

               INTEG "X"

-5.98921311713ᴇ-34   ***
5.98921311713ᴇ-34    ***

-1.31734759935ᴇ-33   ***
7.18426287633ᴇ-34    ***

-1.0279090934ᴇ-33    ***
 2.8943850595ᴇ-34    ***

3.55199210828ᴇ-34    ***
1.38310830422ᴇ-33    ***

2.36684546028ᴇ-34    ***
  1.185146648ᴇ-34    ***

-2.8689321452ᴇ-34    ***
5.23577760548ᴇ-34    ***

7.46276686603ᴇ-34    ***
1.03316990112ᴇ-33    ***

-2.3216434001ᴇ-34    ***
9.78441026613ᴇ-34    ***

-1.01058557743ᴇ-33   ***
7.78421237423ᴇ-34    ***

-6.08235288821ᴇ-34   ***
4.02350288611ᴇ-34    ***

-1.11859154415ᴇ-33   ***
 5.1035625533ᴇ-34    ***

-3.92290932651ᴇ-34   ***
  7.263006115ᴇ-34    ***

-6.63471685897ᴇ-35   ***
3.25943764062ᴇ-34    ***

4.65003711572ᴇ-34    ***
5.31350880161ᴇ-34    ***

5.01704822504ᴇ-34    ***
3.67011109323ᴇ-35    ***

-1.73454285248ᴇ-34   ***
6.75159107752ᴇ-34    ***

[Thomas Okken]

I am a bit surprised because I have observed the death spiral towards zero before, and I even thought it was with this specific test case. But apparently not!

[lrdheat]

Thanks for this! When I suspect this to be the case, I will add a 1 - to the code, and verify that it converges on -1.

[lrdheat]

Adding 1E-10 to the code made the DM 42 converge immediately on 2pi*1E-10 which is what I was expecting on this example...had no awareness of the potential of a zero solution spiral on the DM 42...I believe that I have encountered this on an earlier occasion, but did not attempt to see if this happened on the 42S. The designer of the 42S algorithm must have been aware of this possibility. It appears as if, unfortunately, this was not taken into account on the DM 42!

[lrdheat]

Wondering...would the DM 42 be running the same integration algorithm as the 42S? I’m curious as to if the zero spiral phenomenon is due to a different platform or algorithm. It appears as if when the emphasis of what had been a large HP calculator group was focused more on an engineering audience, a lot of thought went into how each function performed over a full range of scenarios that might be encountered. The DM 42 is very beautiful achievement, but perhaps did not have the luxury of having a large team to optimize all of the functions that are offered. Happy to have both.

[Thomas Okken]

The DM42 firmware is basically Free42. The INTEG algorithm in Free42 is similar to that used in the HP-42S, i.e. Romberg integration with non-uniform spacing of sample points, and very similar performance in terms of convergence and accuracy. The one significant difference is that the HP-42S INTEG algorithm is smarter about when to stop the iteration...

[J-F Garnier]

(02-14-2021 03:27 AM)Thomas Okken Wrote:  The INTEG algorithm in Free42 is similar to that used in the HP-42S, i.e. Romberg integration with non-uniform spacing of sample points, and very similar performance in terms of convergence and accuracy. The one significant difference is that the HP-42S INTEG algorithm is smarter about when to stop the iteration...

The Romberg algorithm, as implemented in HP calculators starting with the 34C has no problem with null integrals by design, there are no special cases.
In HP implementation, the termination condition is based on the sum of the absolute values of the samples. Is it the same in Free42?

We discussed the HP algorithm last year, see here for some indications on the terminating condition.

J-F

[Albert Chan]

(02-13-2021 01:38 AM)Thomas Okken Wrote:  I am a bit surprised because I have observed the death spiral towards zero before ...

It may be hard to see death spiral to 0, with INTEG u-transformation, and romberg extrapolation.
However, if we evaluate integral other ways, it is not hard to show one.

\(\displaystyle AGM(a,b) = {(a+b)(\pi/4) \over K\left(k = {a-b \over a+b}\right)} \)

AGM(a,b) = arithmetic-geometic mean of a, b (above is equation 8)
K = complete elliptic integral of first kind
k = elliptic modulus

(pi/2) / K(1) = AGM(2,0) = AGM(1,0) = AGM(1/2,0) = AGM(1/4,0) = ... = 0

Arithmetic mean forms a death-spiral toward zero.
Thus, it may be better to test geometric mean convergence, AGM(2,0) = AGM(1,0) = 0

[Thomas Okken]

(02-14-2021 08:57 AM)J-F Garnier Wrote:  We discussed the HP algorithm last year, see here for some indications on the terminating condition.

Thanks, I will look into that. It doesn't look like it would help with the death spiral scenario, but for the case at hand, where the approximation just bounces around, a better error estimate would definitely improve things, and allow the iteration to end much more quickly.

[Albert Chan]

A trick to detect bouncing around zero, or (rare) spiraling to zero, is flushing tiny result to 0

With sucessive integral estimates both zeroes, it now can detect convergence.

[Thomas Okken]

The problem with that is that you don't know what "tiny" is.

[Albert Chan]

(02-15-2021 02:06 PM)Thomas Okken Wrote:  The problem with that is that you don't know what "tiny" is.

May be we do, if we estimate size of M = ∫(|f(x)|, x= a .. b)

OP example, F = ∫(sin(x), x= 0 .. 2*pi)

M = ∫(|sin(x)|, x= 0 .. 2*pi) = 4. (rough estimate OK)

That meant rounding errors of M might be maximum of 4e-33 (Free-42 Decimal)
x = a + k*h might also generate some errors, so above estimate of 'tiny' is not too bad.

Any estimated |F| below 'tiny' might be considered 0.