When would you absolutely need recursion that couldn't be fulfilled with a while loop?
It's a nice to have that makes your life easier
What you do is, you do a three-point approximation and a five-point approximation. The difference between the two gives you a fairly good estimate of the error. If the difference is too high, you cut the region in half, and recursively call the same function on each half.
That calling twice is what makes it hard for a while loop. I mean, yes, you could do it with a work queue of intervals or something, but it would be much less straightforward than a recursive call.
local function q(f, start, fstart, mid, fmid, stop, fstop, ε)
-- Scaled approximation of integral on 3 points using Simpson’s
-- rule.
-- <https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas#Closed_Newton%E2%80%93Cotes_formulas>
local Q3 = (1/6)*(fstart + 4*fmid + fstop)
-- Scaled approximation of integral on 5 points.
local mid1, mid2 = (1/2)*(start + mid), (1/2)*(mid + stop)
local fmid1, fmid2 = f(mid1), f(mid2)
-- Boole’s rule, giving exact results for up to cubic polynomials
-- and an O(h⁷) error. This means that every additional level of
-- recursion, adding 1 bit of precision to the independent
-- variable, gives us 7 more bits of precision for smooth
-- functions! So we need like 256 evaluations to get results to
-- 53-bit machine precision?
local Q5 = (7/90)*(fstart+fstop) + (32/90)*(fmid1+fmid2) + (12/90)*fmid
--print(mid, Q3, Q5)
-- Approximate the error by comparing the two.
if (stop - start)*math.abs(Q5 - Q3) <= ε then return (stop - start) * Q5 end
-- Recursively subdivide the interval.
return
q(f, start, fstart, mid1, fmid1, mid, fmid, ε) +
q(f, mid, fmid, mid2, fmid2, stop, fstop, ε)
end
function adaquad(f, start, stop, ε)
local mid = (1/2)*(start + stop)
-- In some quick testing, ε/2 seems to work okay.
return q(f, start, f(start), mid, f(mid), stop, f(stop), (1/2)*ε)
end
return adaquad
It's kind of insane that 16 lines of code can do this!