adaptive Simpson's method¶

basic form¶

The integral I=∫badxf(x)I=∫abdxf(x)

is approximated as

I=S+ES=b−a6(f(a)+4f(m)+f(b))E=−190(b−a2)5f(4)(ξ).I=S+ES=b−a6(f(a)+4f(m)+f(b))E=−190(b−a2)5f(4)(ξ).

for m=a+b2m=a+b2 and some ξξ in within a:b.

adaptive scheme¶

Given the basic form, we can consider two levels of approximations:

S1(a,b)=b−a6(f(a)+4f(m)+f(b))S2(a,b)=S1(a,m)+S1(m,b).S1(a,b)=b−a6(f(a)+4f(m)+f(b))S2(a,b)=S1(a,m)+S1(m,b).

and the respective errors work out to be

E1(a,b)=−190(b−a2)5f(4)(ξ)E2(a,b)=−1902(b−a4)5f(4)(ξ)≈116E1(a,b).E1(a,b)=−190(b−a2)5f(4)(ξ)E2(a,b)=−1902(b−a4)5f(4)(ξ)≈116E1(a,b).

This is expected for a 4-th order scheme. With S1+E1=S2+E2S1+E1=S2+E2 we have a way to construct the error via

E2(a,b)≈115(S2(a,b)−S1(a,b))E2(a,b)≈115(S2(a,b)−S1(a,b))

In [7]:
function adaptive_simpsons(f, a, b; tol=1E-9, imax = 8000)
    
    order = 4
    
    # counter
    i1 = 1
    
    function quad_rule(a1, a2, f1, f2)
        m1 = (a1+a2)/2.0
        fm = f(m1)
        estim = (a2-a1)/6.0 * (f1 + 4.0*fm + f2)
        return m1, fm, estim
    end
    
    function _work(a1, a2, f1, f2, _tol)
        
        # tracker
        i1 += 1
        
        m, fm, estim = quad_rule(a1, a2, f1, f2)
        
        left = quad_rule(a1, m, f1, fm)[3]
        right = quad_rule(m, a2, fm, f2)[3]
        err = left+right - estim
 
        if abs(err) <= (2^order-1) * _tol || i1 > imax
            ret = left+right + err/(2^order-1)
        else
            ret = (
                _work(a1, m, f1, fm, 0.5*_tol) +
                _work(m, a2, fm, f2, 0.5*_tol)
            )
        end
        
        return ret
    end
    
    ret = _work(a, b, f(a), f(b), tol)
    
    return ret
    
end
Out[7]:
adaptive_simpsons (generic function with 1 method)
In [8]:
using BenchmarkTools
println([adaptive_simpsons(x->x*log(1+x), 0.0, 1.0), 0.25])
@benchmark ret1 = adaptive_simpsons(x->x*log(1+x), 0.0, 1.0), 0.25

[0.2500000000002274, 0.25]
Out[8]:
BenchmarkTools.Trial: 10000 samples with 8 evaluations.
 Range (min … max):  3.443 μs …  1.181 ms  ┊ GC (min … max): 0.00% … 99.38%
 Time  (median):     3.787 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   3.991 μs ± 14.325 μs  ┊ GC (mean ± σ):  4.98% ±  1.40%

                      ▂▇█▆▂▁                                  
  ▁▂▁▁▁▁▁▂▄▅▆▅▄▃▂▂▁▁▂▄██████▆▅▄▄▃▃▂▂▂▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  3.44 μs        Histogram: frequency by time        4.33 μs <

 Memory estimate: 4.28 KiB, allocs estimate: 271.
In [ ]: