In [62]:
# useful integration routines

Romberg Integration¶

In [11]:
# Romberg Integration

function romberg(f, a, b; Max=20, tol=1E-10)
    
    function interior(n)
        # for internal points
        dx = (b - a) / 2^(n-1)
        return dx * sum([f(a + (2*i1-1)*dx) for i1 in 1:2^(n-2)])
    end
        
    Rmat = zeros(Max, Max)
    # initial result
    Rmat[1, 1] = 0.5 * (b-a) * (f(a)+f(b))
    estim = Rmat[1, 1]
      
    for i1 in 2:Max
        
        # without duplicated calculations
        Rmat[i1, 1] = 0.5 * Rmat[i1-1, 1] + interior(i1)
        
        for i2 in 2:i1
            # Rmat[i1, i2] = Rmat[i1, i2-1] + (Rmat[i1, i2-1]-Rmat[i1-1, i2-1])/(4^(i2-1)-1)
            Rmat[i1, i2] = (4^(i2-1)*Rmat[i1, i2-1]-Rmat[i1-1, i2-1])/(4^(i2-1)-1)
        end
 
        if abs(estim-Rmat[i1,i1])*(4^(i1-1)-1) < tol || i1 == Max
            i1 == Max ? println("max iteration reached...") : nothing
            break
        else
            estim = Rmat[i1, i1]
        end
            
    end
    
    return estim
    
end
Out[11]:
romberg (generic function with 1 method)
In [7]:
using BenchmarkTools
println([romberg(x->x*log(1+x), 0.0, 1.0), 0.25])
@benchmark ret1 = romberg(x->x*log(1+x), 0.0, 1.0), 0.25

[0.25000000000000006, 0.25]
Out[7]:
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.033 μs … 710.696 μs  ┊ GC (min … max): 0.00% … 99.51%
 Time  (median):     1.279 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.378 μs ±   7.095 μs  ┊ GC (mean ± σ):  5.13% ±  1.00%

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

 Memory estimate: 4.69 KiB, allocs estimate: 8.

Adaptive Simpsons¶

In [8]:
function adaptive_simpsons(f, a, b; tol=1E-10, imax = 2500)
    
    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[8]:
adaptive_simpsons (generic function with 1 method)
In [10]:
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.2500000000000068, 0.25]
Out[10]:
BenchmarkTools.Trial: 10000 samples with 5 evaluations.
 Range (min … max):  6.158 μs …  1.820 ms  ┊ GC (min … max): 0.00% … 99.29%
 Time  (median):     6.717 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   7.078 μs ± 22.320 μs  ┊ GC (mean ± σ):  4.38% ±  1.40%

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

 Memory estimate: 7.53 KiB, allocs estimate: 479.
In [ ]: