HP71B Integral Questions

[Wes Loewer]

I'm still wrapping my head around my very long-term misunderstanding of hp's implementation of the Romberg method. So how does this work to use trapezoids but not evaluate the endpoints? Do you just leave out the two outermost trapezoids?

(If I'm not mistaken, I think you can also reuse previous samples for midpoints if you trisect the intervals.)

[J-F Garnier]

(02-08-2020 10:46 AM)Wes Loewer Wrote:  I'm still wrapping my head around my very long-term misunderstanding of hp's implementation of the Romberg method. So how does this work to use trapezoids but not evaluate the endpoints? Do you just leave out the two outermost trapezoids?

Actually, the endpoints are not really left out, just it's useless to evaluate the function at the endpoints since their contribution to the sum is null due to the (1-u²) factor. So nothing is missing in the sum.

J-F

[Wes Loewer]

(02-08-2020 10:59 AM)J-F Garnier Wrote:  Actually, the endpoints are not really left out, just it's useless to evaluate the function at the endpoints since their contribution to the sum is null due to the (1-u²) factor. So nothing is missing in the sum.

If the outside trapezoids don't contribute to the sum, isn't that equivalent to leaving them out? Wouldn't the approximate area of the outer intervals be missing from the sum?

[Albert Chan]

(02-07-2020 05:08 PM)Albert Chan Wrote:  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

Besides showing INTEGRAL algorithm using romberg's method with trapezoids, I noticed a few things.

• Integrand that evaluated to zero for end points meant the outside trapezoids is approximated as triangles.
  Area is 0 if there no sample points in between. Thus trapezoid area, T₁ = 0
• Even though the midpoint numbers does not incorprate previous sample points, it will, when we do the extrapolation.
  For 7 sample points, midpoints extrapolated estimate = (M1-20*M2+64*M4)/45
• For this example, extrapolated trapezoid numbers seems to converge with half the error.
  We may improve the estimate wih a 2:1 weighted average of both approaches.
  (unfortunately, in general, we have no idea what the weight is)
  
  
  Code:
  

  points  t=trapezoids    m=midpoints     (t+m)/2         (2t+m)/3
7       0.9780170698    1.0269366637    1.0024768667    0.9943236011
15      0.9945010276    1.0108562045    1.0026786161    0.9999527533
31      0.9986339048    1.0027587099    1.0006963073    1.0000088398
63      0.9996599242    1.0006854427    1.0001726835    1.0000017637
127     0.9999151794    1.0001704034    1.0000427914    1.0000002541
255     0.9999788206    1.0000424598    1.0000106402    1.0000000336
511     0.9999947084    1.0000105961    1.0000026523    1.0000000043
1023    0.9999986775    1.0000026466    1.0000006621    1.0000000005
2047    0.9999996694    1.0000006613    1.0000001654    1.0000000001

Note: \(\int _{-1}^1 x \log(1+x)\;dx = 1\)
We can prove using integration by parts, but is easier with symmetry

\(\log(1+x) = \int _0 ^x {dt\over1+t} = \int _0 ^x (1-t+t^2-t^3+t^4+\;...) dt
= {x\over1} - {x^2\over2} + {x^3\over3} - {x^4\over4} +\; ... \)

\(\int _{-1}^1 x \log(1+x)\;dx
= \int _{-1}^1 ({x^2\over1}-{x^3\over2}+{x^4\over3}-{x^5\over4}+\; ... )\; dx
= 2 \int _0 ^1 ({x^2\over1} + {x^4\over3} + {x^6\over5} +\;...) dx
\)

Summing integrated terms, we have a telescoping series

\({2\over1×3} + {2\over3×5} + {2\over5×7} + \;... = {1\over1}-{1\over3} + {1\over3}-{1\over5}+{1\over5}-{1\over7}+\;... = 1\)

[Wes Loewer]

Yes, I see what you mean. I get the same values as your tables.

(02-09-2020 01:43 PM)Albert Chan Wrote:  [*] Even though the midpoint numbers does not incorprate previous sample points, it will, when we do the extrapolation.

True, but the fact that for the midpoint method each new row has to start over from scratch is the very thing that makes it not as good as the trapezoid. If the previous points were reusable, then all midpoint values on your Romberg table would shift up by one row, making it the better option. This is the issue that I was missing in my thinking.

Quote:For this example, extrapolated trapezoid numbers seems to converge with half the error.

This makes perfect sense. The error bounds for midpoint and trapezoid methods are: |Em | < k(b-a)^3/(24n^2) and |Et | < k(b-a)^3/(12n^2)
So for a given number of intervals (n) the midpoint method has about half the error of the trapezoid method. But for a given number of evaluation points, the trapezoid will have twice as many intervals, so now trapezoid will have about half the error.

Just for fun, I wanted to see what would happen with Romberg using midpoints with trisected intervals, thereby allowing midpoints to reuse previous points. (If I did my Romberg calculations right, the only difference was changing 4^p-1 in the denominator to 9^p-1.) Here's what I get:

Code:

points  midpoints with trisection
1       0.000000000
3       1.34471664
9       1.01187573
27      1.00152749
81      1.00016874
243     1.00001866
729     1.00000207
2187    1.00000023

This is only slightly better than the trapezoid method for this example. But for the previously mentioned problematic 1/sqrt(1-x^2), it makes a huge difference

Code:

points   midpoints with trisection
1       3.000000000000
3       3.136485386483
9       3.141533702802
27      3.141592493620
81      3.141592652733

[Wes Loewer]

(02-08-2020 03:04 PM)Wes Loewer Wrote:  

(02-08-2020 10:59 AM)J-F Garnier Wrote:  Actually, the endpoints are not really left out, just it's useless to evaluate the function at the endpoints since their contribution to the sum is null due to the (1-u²) factor. So nothing is missing in the sum.

If the outside trapezoids don't contribute to the sum, isn't that equivalent to leaving them out? Wouldn't the approximate area of the outer intervals be missing from the sum?

I see now why 1/sqrt(1-x^2) has such a difficult time converging with this algorithm. Every time the algorithm goes to a smaller interval, the total width being considered keeps growing. It's like trying to hit a moving target.

For most functions, the cubic transformation will cause the endpoints to approach 0, so the added area has a shrinking height and width. But in this case, the endpoints approach a finite value (~1.73), so while the width is shrinking, the height is not.

This has all been very educational. :-)

[Albert Chan]

(02-09-2020 08:33 PM)Wes Loewer Wrote:  

Quote:For this example, extrapolated trapezoid numbers seems to converge with half the error.

This makes perfect sense.
The error bounds for midpoint and trapezoid methods are: |Em | < k(b-a)^3/(24n^2) and |Et | < k(b-a)^3/(12n^2)
So for a given number of intervals (n) the midpoint method has about half the error of the trapezoid method. But for a given number of evaluation points, the trapezoid will have twice as many intervals, so now trapezoid will have about half the error.

The constant k of |Em| and |Et| are not the same k, but depends on actual curve.
Even if the k's are the same, the formulas only give the upper bound, not equality.

Also, INTEGRAL never return raw trapezoid numbers, but a minimum of 3 extrapolations.
(even more extrapolations on top with more than 7 sample points).

When I try compare errors of 2 methods, error ratios goes all over the place.
The two schemes likely bracketed the true area (but, I also found exceptions)

Example, from Kahan's article, Handheld Calculator Evaluates Integrals, page 30

f(x) = sqrt(x)/(x-1) - 1/log(x), x=0 to 1        → Et/Em ≈ -10
u-transformed f(x), u=-1 to 1                    → Et/Em ≈ -0.4
u-transformed f(x), 2 times                      → Et/Em ≈ -0.06

Equivalent formula, with x=w^2 substitution:

g(w) = 2*w^2/((w+1)*(w-1)) - w/log(w), w=0 to 1  → Et/Em ≈ -0.4
u-transformed g(w), u=-1 to 1                    → Et/Em ≈ -0.05
u-transformed g(w), 2 times                      → Et/Em ≈ -0.02 to +0.02

Note: error ratios only consider final few iterations, initial numbers can go wild.

u-transformation make the curve "bumpy", creating area strips errors of opposite signs, cancelling each other.
If u-transformed more than once, area converge faster (at the cost of more time to evaluate functions).

[Albert Chan]

Trapezoid and midpoint extrapolated numbers are using same points.
Doing both schemes together in parallel are very cheap.

Code:

-- integ.lua

local function tm(f, a, b, t)
    local n, h = 2, (b-a)*0.5
    t = t or (f(a)+f(b))*h  -- trapezoid, T1
    local function iter()       
        while true do
            local m = 0
            for i = 1, n, 2 do m = m + f(a + i*h) end
            m = m*h         -- midpoint area / 2
            t = t*0.5 + m   -- trapezoid area
            coroutine.yield(t, m+m)
            n, h = n*2, h*0.5
        end
    end
    return coroutine.wrap(iter), t
end

local function simpson(f, a, b, t)
    local gen, t0 = tm(f, a, b, t)
    local t1, m1 = gen()    -- line up sequences
    local m0 = m1           -- disable first correction
    local function iter()
        while true do
            coroutine.yield(t1+(t1-t0)/3, m1+(m1-m0)/3)
            t0, m0, t1, m1 = t1, m1, gen()
        end
    end
    return coroutine.wrap(iter)
end

local function romberg(gen, k, t0, m0)
    local function iter()
        for t1, m1 in gen do
            coroutine.yield(t1+(t1-t0)/k, m1+(m1-m0)/k)
            t0, m0 = t1, m1
        end
    end
    return coroutine.wrap(iter)
end

local 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

local function integ(gen, limit)
    limit = limit or 1024       -- default max intervals
    local k, n, t1 = 15, 4
    local t0, m0 = gen()        -- first m0 discarded
    repeat
        t1, m0 = gen()
        t0 = t1 + (t1-t0)/k
        if limit>0 then print(n, t0, m0) end
        gen = romberg(gen, k, t1, m0)
        n, k = 2*n, 4*k + 3
    until limit*limit < k
    return t0, m0
end

return {integ=integ, tm=tm, simpson=simpson, u=u}

Note: midpoints "first" simpson numbers were a placeholder, to line-up the 2 sequences.
(1 sample point cannot produce simpson number, which required 3+ points)

lua> I = require'integ'
lua> log1p = require'mathx'.log1p
lua> function f(x) return x*log1p(x) end
lua> s = I.simpson(I.u(f,-1,1),-1,1,0)
lua> a = I.integ(s)

4     0.9275194244636257   1.7390989208692982
8     0.9780170697704409   1.0269366636614183
16    0.9945010276157439   1.0108562045403804
32    0.9986339047563956   1.002758709871382
64    0.9996599242359795   1.0006854427294893
128   0.9999151793902487   1.0001704033854413
256   0.9999788205588988   1.0000424597853743
512   0.9999947084064521   1.0000105961327912
1024  0.9999986775131146   1.000002646612206

We don't know the error ratios of 2 methods.
But, we can apply Aitken's Extrapolation for the better convergence sequence.
In other words, extrapolate along the romberg table diagonals.

For u-transfomed f(x), trapezoid numbers is better. ("error" between iterations is smaller)
Aitkens extrapolation of last 5 trapezoid numbers gives area = 1.00000 000003

Code:

64    0.99965992424
128   0.99991517939
256   0.99997882056  0.99999995784
512   0.99999470841  0.99999999440
1024  0.99999867751  0.99999999928  1.00000000003

For comparison, INTEGRAL (with built-in u-transform), and maximum 65536 intervals.

>INTEGRAL(-1,1,0,IVAR*LOGP1(IVAR))
.999999 999665

[Albert Chan]

Amazingly, Aitken's extrapolation trick work even for the troublesome 1/sqrt(1-x^2)

lua> I=require'integ'
lua> function f(x) return 1/sqrt(1-x*x) end
lua> u = I.u(f,-1,1)
lua> a = I.integ(I.simpson(u,-1,1,0),2^10) -- with bad assumption, T1 = 0

4     2.6032971925091086   3.131182235954578
8     2.876452993708263    3.141072676119944
16    3.009535987413591    3.1415792702305954
32    3.075628788308323    3.1415925017013078
64    3.1086187748017564   3.141592652903347
128   3.1251067205395437   3.1415926535886425
256   3.1333498128442754   3.141592653589793
512   3.1374712489390655   3.141592653589794
1024  3.139531953229666    3.141592653589784

Code:

Aitken Extrapolation on Trapezoid numbers
64    3.10861877480
128   3.12510672054
256   3.13334981284  3.14159114440
512   3.13747124894  3.14159246493
1024  3.13953195323  3.14159263001  3.14159265359

As expected, with correct T1 value, trapezoid numbers converge better.

lua> I.integ(I.simpson(u,-1,1,2*sqrt(3)),2^10)

4     3.1421574437527586   3.131182235954578
8     3.141606450669424    3.141072676119944
16    3.1415928071550314   3.1415792702305954
32    3.141592654278877    3.1415925017013078
64    3.141592653590944    3.141592652903347
128   3.1415926535897936   3.1415926535886425
256   3.1415926535897936   3.141592653589793
512   3.1415926535897944   3.141592653589794
1024  3.1415926535897887   3.141592653589784

[Albert Chan]

Out of curiosity, I tried using the actual Gamma function integral

\(\Gamma(x+1) = \int _0 ^1 (-\ln(t))^x dt \)

This integral, for any big x, choke on HP71B INTEGRAL

Instead, I tried letting \(t = s^{x+1}\), turning the extreme spike into a bell-shaped curve

\(\Gamma(x+1) = \int _0^1 (x+1) s^x (-(x+1)\ln(s))^x ds = (x+1)^{x+1} \int _0^1 (-s\ln(s))^x ds\)

Scale RHS 2 terms to avoid hitting overflow limits, we have:

Code:

10 P=1/10^8
20 DEF FNA(X)=INTEGRAL(0,1,P,(-LN(IVAR))^X)
30 DEF FNB(X)=100*((X+1)/100)^(X+1)*INTEGRAL(0,1,P,(-100*IVAR*LN(IVAR))^X)
40 INPUT "X=?";X
50 T=TIME @ DISP FNA(X),IBOUND,TIME-T
60 T=TIME @ DISP FNB(X),IBOUND,TIME-T
70 DISP GAMMA(X+1),"GAMMA(X+1)"

>RUN
X=?10
3623769.01194     -3.61915267163E-2     40.36
3628800.00001     12719498.961          .09
3628800           GAMMA(X+1)
>RUN
X=?20
1.66431877492E18  -15440997253.2        40.47
2.4329020082E18   4.16433765069E22      .11
2.43290200818E18  GAMMA(X+1)
>RUN
X=?30
1.39229859877E31  -1.13405116532E23     40.31
2.65252859812E32  1.55401137325E38      .19
2.65252859812E32  GAMMA(X+1)