Post Reply 
Numerical integration vs. integrals that are zero
03-05-2019, 02:21 PM
Post: #1
Numerical integration vs. integrals that are zero
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.
Visit this user's website Find all posts by this user
Quote this message in a reply
03-05-2019, 03:51 PM
Post: #2
RE: Numerical integration vs. integrals that are zero
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.
Find all posts by this user
Quote this message in a reply
03-05-2019, 07:47 PM
Post: #3
RE: Numerical integration vs. integrals that are zero
(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
Find all posts by this user
Quote this message in a reply
03-05-2019, 08:36 PM (This post was last modified: 03-05-2019 09:04 PM by Thomas Okken.)
Post: #4
RE: Numerical integration vs. integrals that are zero
(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)
Visit this user's website Find all posts by this user
Quote this message in a reply
03-06-2019, 10:28 AM
Post: #5
RE: Numerical integration vs. integrals that are zero
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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
Find all posts by this user
Quote this message in a reply
03-06-2019, 11:08 AM
Post: #6
RE: Numerical integration vs. integrals that are zero
"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.)
Visit this user's website Find all posts by this user
Quote this message in a reply
03-06-2019, 12:55 PM
Post: #7
RE: Numerical integration vs. integrals that are zero
(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
...
Find all posts by this user
Quote this message in a reply
03-06-2019, 02:19 PM (This post was last modified: 03-06-2019 03:03 PM by Thomas Klemm.)
Post: #8
RE: Numerical integration vs. integrals that are zero
(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
Find all posts by this user
Quote this message in a reply
03-06-2019, 02:50 PM
Post: #9
RE: Numerical integration vs. integrals that are zero
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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
Find all posts by this user
Quote this message in a reply
03-06-2019, 05:23 PM
Post: #10
RE: Numerical integration vs. integrals that are zero
(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)?

— Ian Abbott
Find all posts by this user
Quote this message in a reply
03-06-2019, 05:27 PM
Post: #11
RE: Numerical integration vs. integrals that are zero
(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(fab):
    
c= (b-a)/4, (b+a)/2
    def t
(uk=3.*c):
        
eps 1.0 u*u      # return 0 if eps=0
        
return eps and k*eps f(c*u*(eps 2.) + d)
    return 


With nonlinear defined, integrate(f, a, b) = integrate(nonlinear(f, a, b), -1, 1)
Find all posts by this user
Quote this message in a reply
03-06-2019, 06:05 PM (This post was last modified: 03-06-2019 06:12 PM by ijabbott.)
Post: #12
RE: Numerical integration vs. integrals that are zero
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.)

— Ian Abbott
Find all posts by this user
Quote this message in a reply
03-06-2019, 06:54 PM (This post was last modified: 03-07-2019 12:49 AM by Albert Chan.)
Post: #13
RE: Numerical integration vs. integrals that are zero
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
Find all posts by this user
Quote this message in a reply
03-06-2019, 10:55 PM
Post: #14
RE: Numerical integration vs. integrals that are zero
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.
Visit this user's website Find all posts by this user
Quote this message in a reply
03-07-2019, 02:42 AM
Post: #15
RE: Numerical integration vs. integrals that are zero
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.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 2 Guest(s)