Numerical integration vs. integrals that are zero

[Thomas Okken]

I noticed that INTEG on the HP-42S terminates suspiciously quickly when evaluating the integral of SIN, in DEG mode, from 0 to 360. It returns zero, which is correct, but how?

What I find puzzling is that you specify the desired accuracy to the HP-42S integrator through the ACC variable, and ACC is a relative error. The thing with integrals that are exactly zero is that while successive approximations get closer and closer to zero, the absolute error shrinks at roughly the same rate as the estimated integral, and so the relative error doesn't improve. With ACC set to, say, 1e-6, you'd expect INTEG to keep running until it starts to underflow -- which would take forever since each successive iteration takes twice as many evaluations of the integrand than the previous one.

Free42 handles this situation by simply giving up after 20 iterations. That means up to 524,287 evaluations of the integrand, which isn't too bad given how fast PCs and smartphones are nowadays. The HP-42S doesn't seem to have such a hard limit, though: I tried integrating RAN in Emu42, and when I gave up, it had already clocked more than 9 million evaluations, which on a real 42S would presumably have taken hours or even days.

I assume what happens is that INTEG doesn't just check the relative error against ACC, but also checks whether the estimate of the integral just keeps getting smaller. I'm curious what the termination condition would be, though. I can't find any indication of this in the HP-42S Owner's Manual, nor in the Programming Examples and Techniques book.

Does anyone have any insights as to how the HP-42S integrator, or other versions of HP INTEG, deal with relative errors vs. integrals that are zero?

P.S. Because of symmetry, every iteration of the integral of SIN from 0 to 360 could be expected to return zero, but INTEG terminates quickly even when the situation isn't quite so tidy, like

Code:

00 { 33-Byte Prgm }
01▸LBL "S"
02 MVAR "X"
03 180
04 RCL "X"
05 X≤Y?
06 GTO 00
07 2
08 ×
09 180
10 -
11 SIN
12 2
13 ×
14 RTN
15▸LBL 00
16 SIN
17 END

from 0 to 270.

[ttw]

Relative error doesn't have much if any meaning when the exact answer is zero. There are several possibilities in the literature but all have problems. I once used (but I can't remember, I'll try to derive it again) a formula that interpolated between relative and absolute error; it computed relative error for large magnitude results and absolute error for small ones.

[Thomas Klemm]

(03-05-2019 02:21 PM)Thomas Okken Wrote:  Does anyone have any insights as to how the HP-42S integrator, or other versions of HP INTEG, deal with relative errors vs. integrals that are zero?

Maybe it is related to this:

Quote:The Romberg terminates when three passes evaluate to the same value under the current display setting.

Just out of curiosity: What happens with the integral of \(\sin(x)\) if you shift the interval from say 7 to 367?
I don't have my HP-42S here to test but the HP-48GX returns -2.61E-12 instead of 0.

Cheers
Thomas

[Thomas Okken]

(03-05-2019 07:47 PM)Thomas Klemm Wrote:  

(03-05-2019 02:21 PM)Thomas Okken Wrote:  Does anyone have any insights as to how the HP-42S integrator, or other versions of HP INTEG, deal with relative errors vs. integrals that are zero?

Maybe it is related to this:

Quote:The Romberg terminates when three passes evaluate to the same value under the current display setting.

The 42S ignores the display setting, but it could check for three passes evaluating to the same result. With limits 0 to 360, I assume every pass evaluates to exactly zero, and it does stop after 3 passes (7 evaluations of the integrand), with an absolute error of 2.19e-4 while ACC=1e-6.

(03-05-2019 07:47 PM)Thomas Klemm Wrote:  Just out of curiosity: What happens with the integral of \(\sin(x)\) if you shift the interval from say 7 to 367?
I don't have my HP-42S here to test but the HP-48GX returns -2.61E-12 instead of 0.

That depends on ACC. With ACC=1e-6, it returns 3.03e-7 with an absolute error of 2.29e-4, so a relative error of 757, after 6 steps, or 63 evaluations of the integrand.

With ACC=1e-3 => -7.56e-4 ± 2.30e-1 (4 steps)
With ACC=1e-6 => 3.03e-7 ± 2.29e-4 (6 steps)
With ACC=1e-9 => 3.42e-11 ± 2.29e-7 (7 steps)
With ACC=1e-12 => 1.19e-10 ± 2.29e-10 (8 steps)
With ACC=1e-15 => 1.19e-10 ± 2.29e-10 (8 steps)

[Werner]

Hi, Thomas.

With SIN(x) from 0 to 360 degrees, the sample points are

first iteration:
180
second iteration:
56.25
303.75
third iteration
113.90625
246.09375
15.46875
344.53125

The sines of which all, unsurprisingly (because of error-free range reduction), sum up to zero, and as already pointed out, three consecutive iterations that fall below the threshold cause the integration to end.

BTW your second example doesn't work for me? With Acc=1e-4, it iterates 8 times (255 samples) and returns 1.74395e-2 as a result. That's hardly suspiciously quick?

Cheers, Werner

[Thomas Okken]

"Suspiciously" quick was a poor choice of words, but what I meant is that it terminates without attaining the requested estimated relative error, and without hitting some arbitrary maximum number of iterations. For integrals like these, Free42 keeps going until it reaches 20 iterations, and I'm looking for a strategy to make it terminate more quickly, but without changing the UI. (Mathematica lets you choose between specifying a relative or absolute error, which would take care of this problem, but the 42S has only the single parameter ACC.)

[Albert Chan]

(03-06-2019 10:28 AM)Werner Wrote:  With SIN(x) from 0 to 360 degrees, the sample points are

first iteration:
180
second iteration:
56.25
303.75
third iteration
113.90625
246.09375
15.46875
344.53125

The sines of which all, unsurprisingly (because of error-free range reduction), sum up to zero ...

What is the formula for fourth, fifth ... iterations ?

Why uneven intervals ? Is it for faster convergence ?

You can still have error-free range reduction with even intervals.
180
90 270
45 315 135 225
...

[Thomas Klemm]

(03-06-2019 12:55 PM)Albert Chan Wrote:  What is the formula for fourth, fifth ... iterations ?

Cf. About the Algorithm

4th iteration
  4.04296875
 33.22265625
 83.49609375
146.42578125
213.57421875
276.50390625
326.77734375
355.95703125

5th iteration
  1.03271484375
  8.89892578125
 23.62060546875
 44.14306640625
 69.41162109375
 98.37158203125
129.96826171875
163.14697265625
196.85302734375
230.03173828125
261.62841796875
290.58837890625
315.85693359375
336.37939453125
351.10107421875
358.96728515625

[Werner]

Albert, the sample points are spaced non-uniformly, by substituting a polynomial of 3rd degree for the variable of substitution - I have no further info on this.
In his paper 'Handheld Calculator evaluates Integral', W. Kahan uses

u = 1.5*v - 0.5*v^3

as substitution for the integral

int(-1,1,f(u),u)

Cheers, Werner

[ijabbott]

(03-06-2019 02:19 PM)Thomas Klemm Wrote:  

(03-06-2019 12:55 PM)Albert Chan Wrote:  What is the formula for fourth, fifth ... iterations ?

Cf. About the Algorithm

4th iteration
  4.04296875
 33.22265625
 83.49609375
146.42578125
213.57421875
276.50390625
326.77734375
355.95703125

5th iteration
  1.03271484375
  8.89892578125
 23.62060546875
 44.14306640625
 69.41162109375
 98.37158203125
129.96826171875
163.14697265625
196.85302734375
230.03173828125
261.62841796875
290.58837890625
315.85693359375
336.37939453125
351.10107421875
358.96728515625

So if I understand correctly, \(v\) is chosen uniformly from the interval \([-1,1]\) (using \(\{0\}\) for the first iteration, \(\{-\frac{1}{2}, \frac{1}{2}\}\) for the second iteration, \(\{-\frac{3}{4}, -\frac{1}{4}, \frac{1}{4}, \frac{3}{4}\}\) for the third iteration, etc.), converted to \(u\) via \(u=\frac{3}{2}v-\frac{1}{2}v^3\) and then converted to \(x\) via \(x=\frac{b-a}{2}u+\frac{b+a}{2}\) (where \(a\) is the lower limit and \(b\) is the upper limit)?

[Albert Chan]

(03-06-2019 02:19 PM)Thomas Klemm Wrote:  Cf. About the Algorithm

Noticed a typo, transformed integral involve 2 substitutions, not 1:

x = (b-a)/2 * v + (b+a)/2, to transform integration range from [a,b]  to [-1,1]
v = 3u/2 - u³/2, nonlinear transformed, integration range from [-1,1] to (-1, 1)

PHP Code:

def nonlinear(f, a, b):
    c, d = (b-a)/4, (b+a)/2
    def t(u, k=3.*c):
        eps = 1.0 - u*u      # return 0 if eps=0
        return eps and k*eps * f(c*u*(eps + 2.) + d)
    return t 


With nonlinear defined, integrate(f, a, b) = integrate(nonlinear(f, a, b), -1, 1)

[ijabbott]

Trying out Google's function graphing capabilities, here is the non-linear transformation:

https://www.google.com/search?q=1.5x-0.5...om+-1+to+1

(EDITED to remove spurious Google tracking stuff.)

[Albert Chan]

Hi, ijabbott

Transformation is actually a bit more involved, not just x, but dx too

x = (a+b)/2 + (b-a)/2 * (3u/2 - u³/2)

dx = 3/4 * (b-a) * (1 - u²) du

With a=0, b=2*Pi, after simplication, we got:

integrate(sin(x), x = 0 to 2*Pi) = integrate(3/2 pi (1-u^2) * sin(pi/2 u(u^2-3)), u = -1 to 1) = 0

integrate(sin(x), x = 0° to 360°) = 180/Pi * 0 = 0

[Thomas Okken]

Trying to get a handle on the behavior of INTEG on the HP-42S, I had it integrate SIN (in DEG mode) from 7 to 367, with ACC = 1, 0.1, 0.01, ..., 1e-15, and charted the results, the estimated absolute error, and the number of points that were sampled:

Code:

    ACC              INTEG                  EPS           PTS

1               -6.63652681894E-1   225.790474144           7
0.1             -6.63652681894E-1    22.5790474144          7
0.01            -7.5600450972E-4      2.2999768795         15
0.001           -7.5600450972E-4      0.22999768795        15
0.0001          -7.5600450972E-4      2.2999768795E-2      15
0.00001          1.872373932E-5       2.29765529718E-3     31
0.000001         3.02829732E-7        2.29202676626E-4     63
0.0000001        1.0434348E-8         2.29202676626E-5     63
0.00000001       3.42144E-11          2.29214420636E-6    127
0.000000001      3.42144E-11          2.29214420636E-7    127
0.0000000001    -6.57E-12             2.29214420636E-8    127
0.00000000001    1.191528E-10         2.29189856856E-9    255
1.E-12           1.194876E-10         2.29189856856E-10   255
1.E-13           1.194876E-10         2.29189856856E-10   255
1.E-14           1.194876E-10         2.29189856856E-10   255
1.E-15           1.194876E-10         2.29189856856E-10   255

I'm not showing the actual sample points; they are distributed as described before by Werner and Thomas K.

There is something odd here: between 1e-6 and 1e-7, the number of points sampled remains the same, 63, but the evaluated integral is different. How can that be? And this happens again between 1e-9 and 1e-10: the same 127 points are sampled, but the results are different. And EPS is clearly also not just a function of the points being sampled.

[Albert Chan]

Excel simulation, sin(x), x = 7 to 367 degree, u transformed integral:

Code:

n    Simpson  Romberg extrapolation
1 -43.872964
3   7.312644  10.725017
7   0.001692  -0.485705  -0.663653
15 -0.000424  -0.000565   0.007136  0.009766
31 -0.000007   0.000021   0.000030  0.000002  -0.000007

It is likely area algorithm is not Romberg extrapolation.