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.

@show num_nodes = length(t);

plot(f,0,4,color=:black,legend=:none,
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)
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"),
plot!(n,order4,l=:dash,label="\$O(n^{-4})\$")