Convergence of trapezoidal integrationΒΆ
We approximate the integral of the function \(f(x)=e^{\sin 7x}\) over the interval \([0,2]\).
using FundamentalsNumericalComputation
f = x -> exp(sin(7*x));
a = 0; b = 2;
In lieu of the exact value, we will use the QuadGK package to find an accurate result.
Q,errest = quadgk(f,a,b,atol=1e-14,rtol=1e-14);
println("Integral = $Q")
Integral = 2.6632197827615394
Here is the error at \(n=40\).
T,t,y = FNC.trapezoid(f,a,b,40)
@show T;
@show err = Q - T;
T = 2.6623029356022876
err = Q - T = 0.0009168471592517768
In order to check the order of accuracy, we double \(n\) a few times and observe how the error decreases.
n = @. 40*2^(0:5)
err = []
for n in n
T,t,y = FNC.trapezoid(f,a,b,n)
push!(err,Q-T)
end
pretty_table((n=n,error=err),backend=:html)
| n | error |
|---|---|
| Int64 | Any |
| 40 | 0.0009168471592517768 |
| 80 | 0.00023006461762786756 |
| 160 | 5.7567555966198114e-5 |
| 320 | 1.4395072795103658e-5 |
| 640 | 3.598966740625542e-6 |
| 1280 | 8.997540876798382e-7 |
Each doubling of \(n\) cuts the error by a factor of about 4, which is consistent with second-order convergence. Another check: the slope on a log-log graph should be \(-2\).
plot(n,abs.(err),m=:o,label="results",
xaxis=(:log10,"n"), yaxis=(:log10,"error"),
title="Convergence of trapezoidal integration")
plot!(n,3e-3*(n/n[1]).^(-2),l=:dash,label="2nd order")