Threaded Mode | Linear Mode

Albert Chan · 02-02-2020, 03:31 PM

Tried on Emu71b, energy integral from thread Sharp EL-506 vs Sharo EL-516T

From the user's manual, about the algoritm, it uses Romberg's method, after a non-linear transform.

Code:

10 A=6371000 @ B=9.46073047258E+15 @ C=3.98589196E+17 @ A1=A-1

20 DEF FNA(X) @ N=N+1 @ FNA=C/EXP(X) @ END DEF

30 DEF FNB(X) @ N=N+1 @ X=EXP(X) @ FNB=C*X/(X+A1)^2 @ END DEF

40 INPUT "P=?";P

50 T=TIME @ N=0 @ I=INTEGRAL(LN(A),LN(B),P,FNA(IVAR))

60 DISP "IA(";N;") =";IROUND(I),IBOUND,TIME-T;"SEC"

70 T=TIME @ N=0 @ I=INTEGRAL(0,LN(B-A1),P,FNB(IVAR))

80 DISP "IB(";N;") =";IROUND(I),IBOUND,TIME-T;"SEC"

>RUN
P=?1E-5
IA( 127 ) = 62563050744 625265.300862 .11 SEC
IB( 255 ) = 62563050665 625640.053047 .3 SEC
>RUN
P=?1E-6
IA( 127 ) = 62563050647 62526.5300862 .11 SEC
IB( 255 ) = 62563050665 62564.0053047 .3 SEC
>RUN
P=?1E-7
IA( 255 ) = 62563050656 6255.39238586 .22 SEC
IB( 255 ) = 62563050665 6256.40053047 .31 SEC

Why does IA(127) have different area, for P=1E-5 vs P=1E-6 ?
Both are using exactly the same 127 sample points.

Does IBOUND (second column) has something to do with different IA(127) results ?
How is IBOUND intepreted ?
At what IBOUND numbers should an integral be splitted ?

Albert Chan · (This post was last modified: 02-03-2020 02:29 PM by Albert Chan.)

(02-02-2020 03:31 PM)Albert Chan Wrote: Why does IA(127) have different area, for P=1E-5 vs P=1E-6 ?
Both are using exactly the same 127 sample points.

Answering myself, first column = trapezoidal area of u-transformed FNA(X) Big Grin

Code:

Intervals       Richardson Extrapolations

1             0

2      25715027    34286703                

4   13721555536 18286835705 19503672306         

8   50890703124 63280418987 66279991205 67022472458     

16  59914609212 62922577908 62898721836 62845050894 62828668848 

32  61913859850 62580276729 62557456650 62552039743 62550890679 ...

64  62401515156 62564066925 62562986272 62563074044 62563117315 ...

128 62522713769 62563113307 62563049732 62563050740 62563050648 ...

Note: u-transformed integrand returned 0 for integrand limit of -1 and 1, thus 1 trapezoid have zero area.

Extending rows of 128 intervals, extrapolate to its limits, we have:

62522713769 Δ
62563113307 40399538
62563049732 -63575
62563050740 1007
62563050648 -91
62563050583 -65 <-\
62563050564 -19 <-- extrapolations does not improve estimate, actually made it worse
62563050559 -5 <-/

The numbers are done in Excel, which use IEEE-754 binary64.
Except for slight rounding errors, the bold numbers matched IA(127) for P=1E-5 and 1E-6.

J-F Garnier · 02-03-2020, 03:00 PM

(02-03-2020 01:58 PM)Albert Chan Wrote: Answering myself, first column = trapezoidal area of u-transformed FNA(X)

Code:

Intervals Richardson Extrapolations 1 0 2 25715027 34286703 4 13721555536 18286835705 19503672306 8 50890703124 63280418987 66279991205 67022472458 16 59914609212 62922577908 62898721836 62845050894 62828668848 32 61913859850 62580276729 62557456650 62552039743 62550890679 ... 64 62401515156 62564066925 62562986272 62563074044 62563117315 ... 128 62522713769 62563113307 62563049732 62563050740 62563050648 ...

The numbers are done in Excel, which use IEEE-754 binary64.
Except for slight rounding errors, the bold numbers matched IA(127) for P=1E-5 and 1E-6.

I'm interested in documenting the Math ROM algorithms, can you please post the Excel formulae you used to reconstruct the INTEGRAL results?

From the HP Journal article, July 1984 p31 on the Math ROM "A Closer Look at INTEGRAL":
"INTEGRAL proceeds by computing a sequence of weighted sums [...]. These sums are accumulated in an extended-precision, eight-level Romberg scheme.[...] These iterations produce a sequence of approximate integrals I(k,j), where k=0,1,... and j=0,1,...,min(k,7). Here, I(k,0) is the direct weighted sum of 2^k samples integration points, and I(k,j) is the Romberg extrapolation given by I(k,j) = (4^j.I(k,j-1) - I(k-1,j-1)) / (4^j-1)."

Is the 8-level Romberg scheme close to the Richardson Extrapolations you used?

About IBOUND:
IBOUND is an estimation of the absolute error of the computed integral. Some details on how it is computed are also mentioned in the article cited above. IBOUND is usually close to the user-provided target precision P times the integral estimation.

J-F

Albert Chan · (This post was last modified: 04-14-2021 12:52 AM by Albert Chan.)

(02-03-2020 03:00 PM)J-F Garnier Wrote: Is the 8-level Romberg scheme close to the Richardson Extrapolations you used?

Romberg's method ≡ Richardson Extrapolations applied to Trapezoids area.

I did the u-transformed trapezoids area in Lua, then copy/paste to Excel

Code:

function u(f, a, b)       -- u-transformed f

    local c, d = (b-a)/4, (b+a)/2

    local k = 3*c

    return function(u)      -- u = (-1, 1)

        local w = (1+u)*(1-u)

        return w==0 and 0 or k*w*f(c*u*(w+2)+d)

    end

end

function trapezoid(f, a, b, t0)

    local n = 2

    local h = (b-a)*0.5

    t0 = t0 or (f(a) + f(b)) * h

    local function iter()

        while true do

            local s = 0

            for i = 1, n, 2 do s = s + f(a + i*h) end

            local t1 = t0*0.5 + s*h

            coroutine.yield(t1)

            t0, n, h = t1, n*2, h*0.5

        end

    end

    return coroutine.wrap(iter)

end

lua> a, b, c = 6.371e6, 9.4607304725808e15, 3.98589196e17
lua> exp, log = math.exp, math.log
lua> function fna(x) return c/exp(x) end
lua> t = trapezoid(u(fna, log(a), log(b)), -1, 1)
lua> for i=1,7 do print(2^i, t()) end
2 25715027.086919405
4 13721555535.80971
8 50890703123.90066
16 59914609212.0083
32 61913859849.80031
64 62401515156.47649
128 62522713769.32443

Only 2 formulas used in Excel, one for the weight, second for the extrapolation. Others are just copies.
With this setup, Simpson's rule have weight of 3, Boole's rule have weight of 15, ...

Code:

   A   B                   C

                         =4*B1+3

 1   0       

 2   25715027.0869194    =B3+(B3-B2)/C$1

 4   13721555535.8097    =B4+(B4-B3)/C$1

 8   50890703123.9006    =B5+(B5-B4)/C$1

 16  59914609212.0083    =B6+(B6-B5)/C$1

 32  61913859849.8003    =B7+(B7-B6)/C$1

 64  62401515156.4764    =B8+(B8-B7)/C$1

 128 62522713769.3244    =B9+(B9-B8)/C$1

Albert Chan · (This post was last modified: 02-03-2020 11:40 PM by Albert Chan.)

(02-02-2020 03:31 PM)Albert Chan Wrote: From the user's manual, about the algoritm, it uses Romberg's method, after a non-linear transform.

There was an error in the manual. x = ½(3u-u³) substitution *required* integrand limit -1 to 1, not a to b

\(\large \int _{-1} ^1 f(x) dx = {3 \over 2}\int _{-1} ^ 1 (1-u^2) f \left(
{ u (3-u^2) \over 2} \right) du \)

For a general case, with integrand limit a to b, we have:

\(\large \int _a ^b f(x) dx = {3(b-a) \over 4}\int _{-1} ^ 1 (1-u^2) f \left(
{ 2(b+a) + (b-a) u (3-u^2) \over 4} \right) du \)

From the manual:

Quote:INTEGRAL uses extended precision. Internally, sums are accumulated in 15-digit numbers ...
up to 65,535 points can be sampled on each subinterval, thus computing the integral to greater precision.

Trivia:
For 2 intervals, h = 1 → 2^16 intervals, h = 1/2^15 = 0.00003 05175 78125 (exact)
With 15 digits precision, and |u| ≤ 1, this implied INTEGRAL internally calculated u's are all exact.

J-F Garnier · 02-05-2020, 08:43 AM

(02-03-2020 11:15 PM)Albert Chan Wrote:
(02-02-2020 03:31 PM)Albert Chan Wrote: From the user's manual, about the algoritm, it uses Romberg's method, after a non-linear transform.

There was an error in the manual. x = ½(3u-u³) substitution *required* integrand limit -1 to 1, not a to b

Thanks, Albert.
Actually, this description originates from a HP-15C manual (text originally written by W. Kahan) ... and was simplified in the HP71 Math ROM manual.
In the HP-15C Advanced Functions Handbook p.41, the integration limits are -1 and 1, and the exact sentence is: "Their spacing can be demonstrated by substituting, say" (emphasis added) to indicate it is an example for illustration only, not an exact detailed description, The important word "say" was removed in the HP75 and HP71 manuals...

Quote:Trivia:
For 2 intervals, h = 1 → 2^16 intervals, h = 1/2^15 = 0.00003 05175 78125 (exact)
With 15 digits precision, and |u| ≤ 1, this implied INTEGRAL internally calculated u's are all exact.

Interesting. The older HP-34C and HP-15C (and HP-41C Advantage ROM) used 13 digits internally, maybe were they more limited in the number of intervals.

J-F

Albert Chan · (This post was last modified: 02-09-2020 11:31 PM by Albert Chan.)

I noticed u-transformation can be applied more than once, accelerate some integrals.

Example: \( \int _{-1} ^1 {dx \over \sqrt{1-x^2}} = 2\;\sin^{-1} 1 = \pi\)

Below, FNB(U) is u-transformed FNA(X), which is more accurate, and run much faster.

Code:

10 INPUT "P=? ";P

20 DEF FNA(X) @ N=N+1 @ FNA=1/SQRT(1-X*X) @ END DEF

30 DEF FNB(U)=1.5*(1-U*U)*FNA(U*(3-U*U)/2)

40 T=TIME @ N=0 @ I=INTEGRAL(-1,1,P,FNA(IVAR))

50 DISP "IA(";N;") =";I,IBOUND,TIME-T;"SEC"

60 T=TIME @ N=0 @ I=INTEGRAL(-1,1,P,FNB(IVAR))

70 DISP "IB(";N;") =";I,IBOUND,TIME-T;"SEC"

>RUN
P=? 1E-4
IA( 8191 ) = 3.14133371246 3.14116984016E-4 6.34 SEC
IB( 31 ) = 3.14159247755 3.14127673304E-4 .03 SEC
>RUN
P=? 1E-5
IA( 65535 ) = 3.14156045534 -3.1415397954E-5 50.97 SEC
IB( 31 ) = 3.14159263336 3.14127673304E-4 .03 SEC

Lets call integrand of IB (= u-transformed FNB) as FNC.
If we plot FNB and FNC, u = -1 to 1, FNB is much flatter, min=FNB(0)=3/2, max=FNB(±1)=√3

Intuition may suggest IA is easier to integrate, but it will be dead wrong.

The code never evaluate the end-points integrand limit, and skip over them.
That's why we only have 2^N-1 samples points, instead of 2^N+1.

→ Trapezoid number, T₁ = 0 (*always*)

FNB, which T₁ should really be 2√3, skewed badly all the romberg numbers.
To confirm, I tried doing by hand, with T₁=2√3, IA(31+2) = 3.14159265428

FNC does not suffer from this, since limit(FNC, U=-1) = limit(FNC, U=+1) = 0

J-F Garnier · 02-05-2020, 06:20 PM

I ran the same tests on the HP75 Math ROM, here is the program adapted from Albert's code for the 75C:

10 A=6371000 @ B=9.46073047258E15 @ C=3.98589196E17 @ A1=A-1
20 DEF FNA(X)
22 N=N+1 @ FNA=C/EXP(X)
25 END DEF
30 DEF FNB(X)
32 N=N+1 @ X=EXP(X) @ FNB=C*X/(X+A1)^2
35 END DEF
40 INPUT "P=?";P
50 T=TIME @ N=0 @ I=INTEGRAL(LOG(A),LOG(B),P,FNA(X))
60 DISP "IA(";N;") =";ROUND(I,0),IBOUND,TIME-T;"SEC"
70 ! T=TIME @ N=0 @ I=INTEGRAL(0,LOG(B-A1),P,FNB(X))
80 ! DISP "IB(";N;") =";ROUND(I,0),IBOUND,TIME-T;"SEC"

The results are a bit different from the HP71 outputs:

P=?2e-2
IA( 63 ) = 62563132323 422758135.629 .11 SEC

P=?1e-3
IA( 127 ) = 62563050559 21501195.6588 .22 SEC

P=?1e-4
IA( 255 ) = 62563050656 2072473.90909 .495 SEC

P=?1e-5
IA( 255 ) = 62563050656 207247.390909 .494 SEC

P=?1e-6
IA( 511 ) = 62563050656 20718.6780627 .989 SEC

For the same target accuracy P, the HP75 uses more samples:
with P=1E-6, the HP71 uses only 127 samples, and the HP75 511 samples.

Comparing with the Excel simulations from Albert for the first two results:

Code:

Intervals       Richardson Extrapolations

64  62401515156    62564066925 62562986272 62563074044 62563117315 62563129267 62563132322

128 62522713769    62563113307 62563049732 62563050740 62563050648 62563050583 62563050564 62563050559

it seems that the HP75 uses the last extrapolated value. In a way it could make sense since each extrapolated value only requires a few operations.

So it turns out that the HP71 and HP75 INTEGRAL are based on the same algorithm but use different criteria for the choice of the number of samples and returned value.

The IBOUND value is still not clear for me, both on the HP71 and HP75. With same number of samples and same extrapolation, why is the IBOUND value different and depends on the user-supplied target accuracy P? For instance in the cases above, IBOUND is exactly divided by 10 when P is changed from 1E-4 to 1E-5 - same samples, same extrapolation, same INTEGRAL result.

J-F

Albert Chan · (This post was last modified: 02-05-2020 08:46 PM by Albert Chan.)

(02-05-2020 06:20 PM)J-F Garnier Wrote: The IBOUND value is still not clear for me, both on the HP71 and HP75. With same number of samples and same extrapolation, why is the IBOUND value different and depends on the user-supplied target accuracy P? For instance in the cases above, IBOUND is exactly divided by 10 when P is changed from 1E-4 to 1E-5 - same samples, same extrapolation, same INTEGRAL result.

I think P is maximum relative error allowed.

Let E = | rough expected integral result | → Maximum absolute error ≈ IBOUND = E × P

Then, it just walk across the romberg numbers, using the difference as a guide, pick the correct gap.

(02-03-2020 01:58 PM)Albert Chan Wrote: Extending rows of 128 intervals, extrapolate to its limits, we have:

62522713769 Δ
62563113307 40399538
62563049732 -63575
62563050740 1007
62563050648 -91
62563050583 -65 <-\
62563050564 -19 <-- extrapolations does not improve estimate, actually made it worse
62563050559 -5 <-/

Say, E = 62563e6

For P=1E-5, 63575 < (IBOUND = 625630) < 40399538 → return 62563050740
For P=1E-6, 1007 < (IBOUND = 62563) < 63575 → return 62563050648

This implied we need at least 4 romberg numbers, or minimum of 7 sample points.

J-F Garnier · (This post was last modified: 02-06-2020 08:57 AM by J-F Garnier.)

Yes, this is the basic idea, I believe.

I found in the HP71 and HP75 user manuals a significant difference in the way IBOUND is computed in the two machines.
This is a simple description, but still helps:
HP71: Convergence is determined using J(k) defined as the kth approximation to the integral of |F| over the same interval of integration.
HP75: Convergence is determined using J(k) defined as the kth approximation to the integral of 10^(int(log|F|)) over the same interval of integration.
In both cases, these estimations are compared to the series of integral extrapolations |I(k,j) - I(k,j-1)| < E.J(k) as you described above, using the user-supplied E value.

Computing 10^(int(log|F|)) may look strange and complicate, but it is not with BCD numbers especially in asm: 10^(int(log|F|)) just means to set the mantissa to 1 and the sign to +, just keeping the exponent. I believe the benefit is to save memory locations, it may come from the original code on the 34C/15C/41C that were quite memory-limited.
From this description, the error estimation is using the weighted samples, before any extrapolation.

This explains the IBOUND results and the difference between the HP71 and HP75, for instance:
HP71:
P=?1e-5
IA( 127 ) = 62563050744 , IBOUND=625265.300862 .
the weighted sample sum was 62526530086, slightly different from the extrapolated value returned by INTEGRAL.

HP75:
P=?1e-3
IA( 127 ) = 62563050559 , IBOUND=21501195.6588
the weighted sample sum (using all mantissa=1) was 21501195659, significantly smaller than the value returned by INTEGRAL.

J-F

Wes Loewer · 02-06-2020, 06:53 PM

(02-05-2020 08:43 AM)J-F Garnier Wrote:
(02-03-2020 11:15 PM)Albert Chan Wrote: There was an error in the manual. x = ½(3u-u³) substitution ...
Actually, this description originates from a HP-15C manual (text originally written by W. Kahan) ... and was simplified in the HP71 Math ROM manual.

The 71b and 15c manuals' descriptions both quote from Kahan's article in the August 1980 Hewlett Packard Journal description of the 34c [https://www.hpl.hp.com/hpjournal/pdfs/Is...980-08.pdf]

A couple of notes on the 3/2*u-1/2*u³ substitution.

1) The manuals and article say something like

Quote:Besides suppressing resonance, the substitution has two additional benefits. First, no sample need be taken from either endpoint of the interval of integration...

If I understand correctly (and please correct me if I'm wrong), it's not this substitution that prevents the end points from being evaluated. The reason the endpoints are not evaluated is because the Romberg variation used is based on rectangle midpoint evaluations rather than something like a trapezoid (trapezium) method which does evaluate endpoints. If you were evaluating at the endpoints of -1 and 1, then after this substitution you would still have the endpoints -1 and 1. (This invariance of the endpoints is part of why this substitution is so useful.)

2) Here's an interesting tidbit. I've never seen this in print, but a few years ago I discovered that the ti-89/92/nspire also use this substitution even though the intervals are already nonuniform beforehand.

The ti-83/84 calculators use the standard 15 point Gauss-Kronrod method which uses nonuniform intervals at the following nodes:

Code:

 0.00000000000000

±0.20778495500790

±0.40584515137740

±0.58608723546769

±0.74153118559939

±0.86486442335977

±0.94910791234276

±0.99145537112081

However, the ti-89/92/nspire numeric integration uses the following nodes which didn't match anything I had ever seen.

Code:

 0.00000000000000

±0.30719191764839

±0.57534429140662

±0.77847088404600

±0.90842445832523

±0.97384146143897

±0.99618089849105

±0.99989079590057

I couldn't figure out where these ti-89 values came from until one day on a whim I plugged the Gauss-Kronrod values into 3/2*u-1/2*u³ and, voila!, out came the ti-89 values.

So then the question was "Why make this substitution if the intervals are already nonuniform?" The hp manuals and article go on to say that this substitution allows you to integrate functions whose slope is infinite at an endpoint. While this is not required since the algorithm is not evaluating at the endpoints anyway, experimentation shows that it indeed gives better results when integrating functions like sqrt(1-u) or sqrt(1-u^2) from -1 to 1 where the slope is infinite at one or more endpoints.

By not using the standard Gauss-Kronrod nodes, the algorithm is not as accurate on higher order polynomials (degree > 14), but is more accurate on functions with square roots. Seems like a good trade off.

Why this happens is illustrated by the semi-circular graph of √(1-X^2) before and after the substitution.

Albert Chan · 02-06-2020, 11:16 PM

(02-06-2020 06:53 PM)Wes Loewer Wrote: A couple of notes on the 3/2*u-1/2*u³ substitution.

If I understand correctly (and please correct me if I'm wrong), it's not this substitution that prevents the end points from being evaluated. The reason the endpoints are not evaluated is because the Romberg variation used is based on rectangle midpoint evaluations rather than something like a trapezoid (trapezium) method which does evaluate endpoints.

A simple test showed that HP71B INTEGRAL do extrapolations from trapezoids, not mid-points rule.

>10 DISP INTEGRAL(-1, 1, 1E-5, 1/SQRT(1-IVAR^2)), IBOUND
>RUN
3.14156045534 -3.1415397954E-5

this failure to converge (65535 sample points !) is due to missing end points evaluation.

\(\large \int _{-1} ^1 f(x) dx
= \int _{-1} ^ 1 {3 \over 2}(1-u^2) f \left({ u (3-u^2) \over 2} \right) du
= \int _{-1} ^ 1 g(u) du \)

\(f(x) = {1\over\sqrt{1-x^2}}\quad → \quad
g(u) = {3 (1-u^2) \over \sqrt{(4-u^2)(1-u^2)^2}}\)

\(\displaystyle{\lim_{u^2 \to 1^-} g(u)} = \displaystyle{\lim_{u^2 \to 1^-}{3\over\sqrt{4-u^2}}} = \sqrt3 ≠ 0\)

Quote:So then the question was "Why make this substitution if the intervals are already nonuniform?"

It does matter, even if gaussian quadrature do not evaluate end points.
u-transformation changes the overall shape of the integrand, affecting all sample points.

try above f(x) and g(u) (with 1-u^2 cancelled)

XCas> f(x) := 1/sqrt(1-x*x)
XCas> g(x) := 3/sqrt(4-x*x)
XCas> time(gaussquad(f(x), x = -1 .. 1)) → 0.0065
XCas> time(gaussquad(g(x), x = -1 .. 1)) → 0.00054

Wes Loewer · (This post was last modified: 02-08-2020 10:35 AM by Wes Loewer.)

(02-06-2020 11:16 PM)Albert Chan Wrote:
(02-06-2020 06:53 PM)Wes Loewer Wrote: If I understand correctly (and please correct me if I'm wrong), it's not this substitution that prevents the end points from being evaluated. The reason the endpoints are not evaluated is because the Romberg variation used is based on rectangle midpoint evaluations rather than something like a trapezoid (trapezium) method which does evaluate endpoints.

A simple test showed that HP71B INTEGRAL do extrapolations from trapezoids, not mid-points rule.

>10 DISP INTEGRAL(-1, 1, 1E-5, 1/SQRT(1-IVAR^2)), IBOUND
>RUN
3.14156045534 -3.1415397954E-5

this failure to converge (65535 sample points !) is due to missing end points evaluation.

But the standard trapezoid method evaluates the endpoints while the rectangle midpoint method does not. As you said, the endpoints are not evaluated, which points to the midpoint method. The 50g definitely uses Romberg extrapolation of midpoints (see edit below) and gives very similar results as above. Setting the number format to FIX 5, you get:

result = 3.14156045554, IERR = -3.14153979546E-5

I could be wrong, but I would be very surprised if they changed up the algorithm for the 71b.

(02-06-2020 11:16 PM)Albert Chan Wrote:
Quote:So then the question was "Why make this substitution if the intervals are already nonuniform?"

It does matter, even if gaussian quadrature do not evaluate end points.
u-transformation changes the overall shape of the integrand, affecting all sample points.

Right, that's just what I was trying to explain and demonstrate in the sample graph.

EDIT: Ack! It seems that my statement "The 50g definitely uses Romberg extrapolation of midpoints" was definitely wrong. :-(
Mea culpa

Albert Chan · (This post was last modified: 02-07-2020 11:09 PM by Albert Chan.)

(02-07-2020 03:49 AM)Wes Loewer Wrote: But the standard trapezoid method evaluates the endpoints while the rectangle midpoint method does not. As you said, the endpoints are not evaluated, which points to the midpoint method. The 50g definitely uses Romberg extrapolation of midpoints and gives very similar results as above. Setting the number format to FIX 5, you get:

result = 3.14156045554, IERR = -3.14153979546E-5

I could be wrong, but I would be very surprised if they changed up the algorithm for the 71b.

Sorry for the confusion.
I thought 65535 sample points failed convergence is enough to show INTEGRAL uses trapezoids, not midpoint rectangles.
Here are the actual numbers, and its extrapolations.

Midpoint numbers, for g(u)=3/sqrt(4-u^2), u = -1 to 1, converged quickly.
Note: table required 127 sample points, since previous iteration numbers cannot be reused.
Note: we can get midpoint numbers from trapezoids: M_n = 2 T_2n - T_n, even if T₁ were wrong.

Code:

1     3                       

2     3.098386677 3.131182236                 

4     3.129937562 3.140454524 3.141072676             

8     3.138610555 3.141501553 3.141571355 3.141579270         

16    3.140842484 3.141586460 3.141592120 3.141592450 3.141592502     

32    3.141404814 3.141592257 3.141592644 3.141592652 3.141592653 3.141592653 

64    3.141545675 3.141592629 3.141592653 3.141592654 3.141592654 3.141592654 3.141592654

Trapezoids, with assumption of T₁ = 0, failed convergence.
Last result = 3.141560455, matching HP71B INTEGRAL result.

Code:

1     0                           

2     1.5         2                       

4     2.299193338 2.565591118 2.603297193                 

8     2.71456545  2.853022821 2.872184934 2.876452994             

16    2.926588003 2.997262187 3.006878145 3.009016132 3.009535987         

32    3.033715243 3.069424323 3.074235132 3.075304291 3.075564245 3.075628788     

64    3.087560028 3.10550829  3.107913888 3.108448471 3.108578449 3.108610721 3.108618775 

128   3.114552852 3.123550459 3.124753271 3.125020563 3.125085551 3.125101687 3.125105714 3.125106721

256   3.12806688  3.132571556 3.133172962 3.133306608 3.133339102 3.13334717  3.133349184 3.133349687

512   3.134828298 3.137082105 3.137382808 3.137449631 3.137465878 3.137469912 3.137470919 3.13747117

1024  3.138210109 3.139337379 3.139487731 3.139521142 3.139529266 3.139531283 3.139531786 3.139531912

2048  3.139901289 3.140465016 3.140540192 3.140556898 3.14056096  3.140561968 3.14056222  3.140562283

4096  3.140746949 3.141028835 3.141066423 3.141074776 3.141076807 3.141077311 3.141077437 3.141077468

8192  3.141169795 3.141310744 3.141329538 3.141333715 3.14133473  3.141334982 3.141335045 3.141335061

16384 3.141381223 3.141451699 3.141461096 3.141463184 3.141463692 3.141463818 3.141463849 3.141463857

32768 3.141486938 3.141522176 3.141526875 3.141527919 3.141528173 3.141528236 3.141528251 3.141528255

65536 3.141539796 3.141557415 3.141559764 3.141560286 3.141560413 3.141560445 3.141560453 3.141560455

The reason for skipping end points is just an optimization, which work for most cases.
If f(x) is finite at the end points, g(u) has end-points = 0, because of (1-u^2) factor.
If f(x) is infinite at the end points, but grow slower than the shrinking (1-u^2), we still have 0 at the end-points.

For other cases, well, at least it returned a negative IBOUND ... Big Grin

J-F Garnier · (This post was last modified: 02-07-2020 08:29 AM by J-F Garnier.)

(02-07-2020 03:49 AM)Wes Loewer Wrote: The 50g definitely uses Romberg extrapolation of midpoints.
...
I could be wrong, but I would be very surprised if they changed up the algorithm for the 71b.

It's easy to check that the HP71 doesn't use the mid points but the trapezoid nodes, by tracing the function sampled points:

10 DEF FNF(X) @ PRINT X; @ FNF=X @ END DEF
20 I=INTEGRAL(-1,1,.00001,FNF(IVAR))

>RUN
0 -.6875 .6875 -.3671875 .3671875 -.9140625 .9140625
the sampled points corresponds to the u values of 0, +/-0.5, +/-0.25, +/-0.75.

I would be interested to know how the 48/49/50 series work. Like the HP71 math ROM, or with further improvements? (I could try but I'm not a RPL guy :-)

J-F

Wes Loewer · 02-07-2020, 01:19 PM

(02-07-2020 08:23 AM)J-F Garnier Wrote:
(02-07-2020 03:49 AM)Wes Loewer Wrote: The 50g definitely uses Romberg extrapolation of midpoints.
...
I could be wrong, but I would be very surprised if they changed up the algorithm for the 71b.

It's easy to check that the HP71 doesn't use the mid points but the trapezoid nodes, by tracing the function sampled points:

10 DEF FNF(X) @ PRINT X; @ FNF=X @ END DEF
20 I=INTEGRAL(-1,1,.00001,FNF(IVAR))

>RUN
0 -.6875 .6875 -.3671875 .3671875 -.9140625 .9140625
the sampled points corresponds to the u values of 0, +/-0.5, +/-0.25, +/-0.75.

Yes, those are the points that you would expect using midpoints. These are the same points used on the 50g.

n=1: the rectangle goes from -1 to 1 with a midpoint at 0
n=2: rectangles go from -1 to 0 and 0 to 1 with midpoints at -0.5, 0.5
n=4: rectangles go from -1 to -0.5, -0.5 to 0, 0 to 0.5, and 0.5 to 1 with midpoints at -0.75, -0.25, 0.25, 0.75

Albert Chan · (This post was last modified: 02-08-2020 02:33 AM by Albert Chan.)

(02-07-2020 01:19 PM)Wes Loewer Wrote: Yes, those are the points that you would expect using midpoints. These are the same points used on the 50g.

Since trapezoids and midpoints are using the same points, I tried to tease out the algorithm, using f(x) = x*log(1+x)

10 DEF FNF(X)=X*LOGP1(X)
20 DEF FNG(U)=1.5*(1-U*U)*FNF(U*(3-U*U)/2)
30 M1=2*FNG(0)
40 M2=FNG(-.5)+FNG(.5)
50 M4=.5*(FNG(-.75)+FNG(-.25)+FNG(.25)+FNG(.75))
60 T1=0
70 T2=(T1+M1)/2
80 T4=(T2+M2)/2
90 T8=(T4+M4)/2
100 DISP "MIDPOINTS=";(M1-20*M2+64*M4)/45
110 DISP "TRAPEZOID=";(-T1+84*T2-1344*T4+4096*T8)/2835
120 DISP "INTEGRAL =";INTEGRAL(-1,1,1,FNF(IVAR))

>RUN
MIDPOINTS= 1.02693666365
TRAPEZOID= .978017069767
INTEGRAL7= .97801706977

For this case, T₁ = 0 is a valid assumption.

>FNG(-.9999), FNG(+.9999)
5.40429438145E-3 2.07933751591E-4

For a detail comparison between these 2 methods: https://math.stackexchange.com/questions...int/674350

J-F Garnier · (This post was last modified: 02-07-2020 05:57 PM by J-F Garnier.)

(02-07-2020 01:19 PM)Wes Loewer Wrote:
(02-07-2020 08:23 AM)J-F Garnier Wrote: the sampled points corresponds to the u values of 0, +/-0.5, +/-0.25, +/-0.75.

Yes, those are the points that you would expect using midpoints. These are the same points used on the 50g.

n=1: the rectangle goes from -1 to 1 with a midpoint at 0
n=2: rectangles go from -1 to 0 and 0 to 1 with midpoints at -0.5, 0.5
n=4: rectangles go from -1 to -0.5, -0.5 to 0, 0 to 0.5, and 0.5 to 1 with midpoints at -0.75, -0.25, 0.25, 0.75

I understand that the 7 values are a single (combined) set of samples, with 8 regions from -1 to -0.75, -0.75 to -0.5 etc.
For these 8 regions, midpoints would be -0.875, -0.625, etc.

The key point may be that from from one iteration to another, the HP71 doesn't compute a new integral forgetting the previous samples, but combines them with new samples (according to the HP articles/manuals).

J-F

Wes Loewer · 02-07-2020, 08:12 PM

(02-06-2020 11:16 PM)Albert Chan Wrote: A simple test showed that HP71B INTEGRAL do extrapolations from trapezoids, not mid-points rule.

>10 DISP INTEGRAL(-1, 1, 1E-5, 1/SQRT(1-IVAR^2)), IBOUND
>RUN
3.14156045534 -3.1415397954E-5

this failure to converge (65535 sample points !) is due to missing end points evaluation.

The 50g also chokes on this integral:
FIX 3: result = 3.13952114413, IERR = 3.13821310172E-3, N = 1023
FIX 4: result = 3.14133371247, IERR = 3.14116984018E-4, N = 8191
FIX 5: bails out after N = 10259

It's very well possible that I have simply believed incorrectly for over a dozen years that the Romberg method hp used involved midpoints rather than trapezoids. I went back and reread the Kahan's original article where he writes, "The simplest method is the midpoint rule, whose nodes all lie in the middles of panels all of the same width" along with the statement that "no sample need be drawn from either end of the interval of integration." I guess I jumped to the conclusion that midpoints were used.

I sure would love to see the source code on this.

Wes Loewer · 02-07-2020, 08:16 PM

(02-07-2020 05:54 PM)J-F Garnier Wrote: The key point may be that from from one iteration to another, the HP71 doesn't compute a new integral forgetting the previous samples, but combines them with new samples (according to the HP articles/manuals).

Good point. I picked up on that when I reread the article.