In [1]:
using BenchmarkTools
adaptive Simpson's method¶
basic form¶
The integral $I = \int_a^b dx f(x)$
is approximated as
$$ \begin{align} 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{align} $$
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 [2]:
function adapt_simp(f, a, b, tol=1E-12, max_depth=80)
rule(a1, a2) = (a2-a1) * (f(a1)+f(a2)+4*f((a1+a2)/2))/6
function adapt(a1, a2, tol, depth)
m = (a1+a2)/2
s1 = rule(a1, a2)
s2 = rule(a1, m) + rule(m, a2)
err = abs(s2-s1)
if err < tol && depth < max_depth
return s2
elseif depth >= max_depth
@warn "max. depth $max_depth reached."
return s2
else
tol = max(eps(), tol/2)
return adapt(a1, m, tol, depth+1) + adapt(m, a2, tol, depth+1)
end
end
return adapt(a, b, tol, 0)
end
Out[2]:
adapt_simp (generic function with 3 methods)
In [3]:
println([adapt_simp(x->tan(x), 0.0, 1.0), -log(cos(1))])
println([adapt_simp(x->tanh(x), 0.0, 1.0), log(cosh(1))])
[0.6156264703860375, 0.6156264703860141] [0.433780830483044, 0.4337808304830271]
In [4]:
@benchmark adapt_simp(x->tanh(x), 0.0, 1.0)
Out[4]:
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample. Range (min … max): 24.791 μs … 1.766 ms ┊ GC (min … max): 0.00% … 97.60% Time (median): 25.291 μs ┊ GC (median): 0.00% Time (mean ± σ): 26.453 μs ± 18.672 μs ┊ GC (mean ± σ): 1.21% ± 2.91% ▅██▇▅▂▁ ▁▁▂▁▁ ▁▃▅▆▅▃▁ ▂ ▇█████████████████████████████▅▅▄▄▄▆▅▆▆▅▅▅▃▄▆▅▅▆▇▆▄▅▃▄▅▁▄▄▄ █ 24.8 μs Histogram: log(frequency) by time 32.1 μs < Memory estimate: 30.61 KiB, allocs estimate: 1958.
In [5]:
function fast_simp(f, a, b, tol=1e-12, max_depth=80)
# not that much faster
function rule(a1, a2, a3, f1, f2, f3)
return (a3-a1)*(f1+f3+4*f2)/6
end
function adapt(a1, a2, a3, f1, f2, f3, tol, depth)
m12 = (a1+a2)/2; fm12 = f(m12)
m23 = (a2+a3)/2; fm23 = f(m23)
s1 = rule(a1, a2, a3, f1, f2, f3)
s2 = rule(a1, m12, a2, f1, fm12, f2)
s2 += rule(a2, m23, a3, f2, fm23, f3)
err = abs(s2 - s1)
ret = s2
if err < tol || depth >= max_depth
depth >= max_depth && @warn "max depth $max_depth reached"
return ret
else
tol = max(tol/2, eps())
return adapt(a1, m12, a2, f1, fm12, f2, tol, depth+1) +
adapt(a2, m23, a3, f2, fm23, f3, tol, depth+1)
end
end
mid = (a+b)/2
fa = f(a)
fb = f(b)
fm = f(mid)
return adapt(a, mid, b, fa, fm, fb, tol, 0)
end
Out[5]:
fast_simp (generic function with 3 methods)
In [6]:
function integrate(f, a, b, tol=1e-12, max_depth=80)
function rule(a1, a2, a3, f1, f2, f3)
ret1 = (a3-a1)*(f1+f3+4*f2)/6
ret2 = (a2-a1)*(f1+f2+4*f((a1+a2)/2))/6
ret2 += (a3-a2)*(f2+f3+4*f((a2+a3)/2))/6
return (16*ret2 - ret1)/15
end
function adapt(a1, a2, a3, f1, f2, f3, tol, depth)
m12 = (a1+a2)/2
m23 = (a2+a3)/2
fm12 = f(m12)
fm23 = f(m23)
s1 = rule(a1, a2, a3, f1, f2, f3)
s2 = rule(a1, m12, a2, f1, fm12, f2)
s2 += rule(a2, m23, a3, f2, fm23, f3)
err = abs(s2 - s1)
# ret = s2
# further improvement
ret = (32*s2-s1)/31
if err < tol || depth >= max_depth
depth >= max_depth && @warn "max depth $max_depth reached"
return ret
else
tol = max(tol/2, eps())
return adapt(a1, m12, a2, f1, fm12, f2, tol, depth+1) +
adapt(a2, m23, a3, f2, fm23, f3, tol, depth+1)
end
end
m = (a+b)/2
fa = f(a)
fb = f(b)
fm = f(m)
return adapt(a, m, b, fa, fm, fb, tol, 0)
end
Out[6]:
integrate (generic function with 3 methods)
In [7]:
println([fast_simp(x->tanh(x), 0.0, 1.0), integrate(x->tanh(x), 0.0, 1.0), log(cosh(1))])
println([fast_simp(x->tan(x), 0.0, 1.0), integrate(x->tan(x), 0.0, 1.0), -log(cos(1))])
[0.433780830483044, 0.4337808304830294, 0.4337808304830271] [0.6156264703860375, 0.6156264703860105, 0.6156264703860141]
In [8]:
@benchmark fast_simp(x->tanh(x), 0.0, 1.0)
Out[8]:
BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample. Range (min … max): 22.083 μs … 92.583 μs ┊ GC (min … max): 0.00% … 51.17% Time (median): 23.208 μs ┊ GC (median): 0.00% Time (mean ± σ): 24.071 μs ± 3.889 μs ┊ GC (mean ± σ): 0.93% ± 4.25% ▄██▅▂ ▂▄███████▇▇▆▅▅▄▄▄▃▃▃▃▃▃▅▆▇▇▆▅▅▄▄▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▂▁▂ ▃ 22.1 μs Histogram: frequency by time 29.5 μs < Memory estimate: 61.17 KiB, allocs estimate: 3914.
In [9]:
@benchmark integrate(x->tanh(x), 0.0, 1.0)
Out[9]:
BenchmarkTools.Trial: 10000 samples with 8 evaluations per sample. Range (min … max): 3.438 μs … 12.177 μs ┊ GC (min … max): 0.00% … 66.90% Time (median): 3.510 μs ┊ GC (median): 0.00% Time (mean ± σ): 3.594 μs ± 442.488 ns ┊ GC (mean ± σ): 0.59% ± 3.46% ▅█▇▆▅▃▂▁ ▃▄▂ ▂ ▇█████████▇▇▇▇▆▆▆▇▆▇██████▇▅▆▅▄▄▄▂▄▄▄▄▃▄▄▄▄▄▄▄▄▂▃▃▃▄▅▅▅▂▄▂▃ █ 3.44 μs Histogram: log(frequency) by time 4.67 μs < Memory estimate: 5.42 KiB, allocs estimate: 346.
compare to Romberg Integration¶
In [10]:
function romberg(f, a, b; Max=25, tol=1E-12)
Rmat = Matrix{Float64}(undef, Max, Max) # Avoids zero-initialization
Rmat[1, 1] = 0.5 * (b - a) * (f(a) + f(b))
estim = Rmat[1, 1]
for mm in 2:Max
# Compute next trapezoidal estimate without redundant calculations
dx = (b - a) / 2^(mm - 1)
inside = dx * sum(f(a + (2 * kk - 1) * dx) for kk in 1:2^(mm - 2))
Rmat[mm, 1] = 0.5 * Rmat[mm - 1, 1] + inside
# Perform Richardson extrapolation
for nn in 2:mm
Rmat[mm, nn] = (4^(nn - 1) * Rmat[mm, nn - 1] - Rmat[mm - 1, nn - 1]) / (4^(nn - 1) - 1)
end
# Convergence check
if abs(Rmat[mm, mm] - Rmat[mm - 1, mm - 1]) < tol * abs(Rmat[mm, mm]+1)
return Rmat[mm, mm]
end
end
@warn "Warning: Maximum iterations ($Max) reached. Result may not be accurate."
return Rmat[Max, Max]
end
Out[10]:
romberg (generic function with 1 method)
In [11]:
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[11]:
BenchmarkTools.Trial: 10000 samples with 252 evaluations per sample. Range (min … max): 301.754 ns … 743.552 ns ┊ GC (min … max): 0.00% … 42.63% Time (median): 319.940 ns ┊ GC (median): 0.00% Time (mean ± σ): 352.902 ns ± 87.169 ns ┊ GC (mean ± σ): 8.11% ± 13.62% █▇▆▅▅▄▆▇▄▃▂▂▁▁▂▁ ▁▃▃▂▂▁▁▁▁ ▁▁▁▁ ▂ ████████████████▆▆▆▄▅▄▃▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄███████████████ █ 302 ns Histogram: log(frequency) by time 635 ns < Memory estimate: 5.08 KiB, allocs estimate: 3.
In [ ]: