Adaptive integration¶
using FundamentalsNumericalComputation
f = x -> (x+1)^2*cos((2*x+1)/(x-4.3));
@show exact,errest = quadgk(f,0,4,atol=1e-14,rtol=1e-14); # 'exact' value
(exact, errest) = quadgk(f, 0, 4, atol = 1.0e-14, rtol = 1.0e-14) = (-2.825533373437446, 2.1510571102112408e-14)
We perform the integration and show the nodes selected underneath the curve.
Q,t = FNC.intadapt(f,0,4,0.001)
@show num_nodes = length(t);
plot(f,0,4,color=:black,legend=:none,
xlabel="x", ylabel="f(x)", title="Adaptive node selection")
plot!(t,f.(t),seriestype=:sticks,m=(:o,2))
num_nodes = length(t) = 69
The error turns out to be a bit more than we requested. It’s only an estimate, not a guarantee.
@show err = exact - Q;
err = exact - Q = -0.022002813037627078
Let’s see how the number of integrand evaluations and the error vary with the requested tolerance.
table = (tol=[],err=[],n=[])
for tol in 10.0.^(-4:-1:-14)
Q,t = FNC.intadapt(f,0,4,tol)
push!(table.tol,tol)
push!(table.err,exact-Q)
push!(table.n,length(t))
end
pretty_table(table,["tolerance","error","f-evaluations"],backend=:html)
| tolerance | error | f-evaluations |
|---|---|---|
| 0.0001 | -0.00041946881777477074 | 113 |
| 1.0e-5 | 4.7897662537543795e-5 | 181 |
| 1.0e-6 | 6.314383849126415e-6 | 297 |
| 1.0e-7 | -6.639249763296107e-7 | 489 |
| 1.0e-8 | 7.18080701567203e-8 | 757 |
| 1.0e-9 | 1.2652385805722588e-8 | 1193 |
| 1.0e-10 | -8.441261023506286e-10 | 2009 |
| 1.0e-11 | 2.6127100483108734e-11 | 3157 |
| 1.0e-12 | 4.044631296551415e-11 | 4797 |
| 1.0e-13 | -1.937117133365973e-12 | 7997 |
| 1.0e-14 | 1.6298074001497298e-13 | 12609 |
As you can see, even though the errors are not less than the estimates, the two columns decrease in tandem. If we consider now the convergence not in \(h\) (which is poorly defined) but in the number of nodes actually chosen, we come close to the fourth order accuracy of the underlying Simpson scheme.
n = table.n
plot(n,abs.(table.err),m=:o,label="results",
xaxis=(:log10,"number of nodes"), yaxis=(:log10,"error"),
title="Convergence of adaptive quadrature")
order4 = @. 0.01*(n/n[1])^(-4)
plot!(n,order4,l=:dash,label="\$O(n^{-4})\$")