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