adaptive Simpson's method¶

basic form¶

The integral $I = \int_a^b dx f(x)$

is approximated as

$ \begin{aligned} I &= S + E \\ S &= \frac{b-a}{6} \, \left( f(a) + 4 f(m) + f(b) \right) \\ E &= -\frac{1}{90} \left( \frac{b-a}{2} \right)^5 \, f^{(4)}(\xi). \end{aligned} $

for $m = \frac{a+b}{2}$ and some $\xi$ in within a:b.

adaptive scheme¶

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

$ \begin{aligned} S_1(a, b) &= \frac{b-a}{6} \, \left( f(a) + 4 f(m) + f(b) \right) \\ S_2(a, b) &= S_1(a, m) + S_1(m, b). \end{aligned} $

and the respective errors work out to be

$ \begin{aligned} E_1(a, b) &= -\frac{1}{90} \left( \frac{b-a}{2} \right)^5 \, f^{(4)}(\xi) \\ E_2(a, b) &= -\frac{1}{90} \, 2\, \left( \frac{b-a}{4} \right)^5 \, f^{(4)}(\xi) \\ &\approx \frac{1}{16} \, E_1(a, b). \end{aligned} $

This is expected for a 4-th order scheme. With $ S_1 + E_1 = S_2 + E_2 $ we have a way to construct the error via

$ E_2(a, b) \approx \frac{1}{15} \, \left( S_2(a,b)-S_1(a,b) \right) $

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 [ ]: