HP71B Integral Questions

[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