Extrapolation for numerical integration

We estimate \(\displaystyle\int_0^2 x^2 e^{-2x}\, dx\) using extrapolation.

using FundamentalsNumericalComputation
f = x -> x^2*exp(-2*x);
a = 0;  b = 2; 
Q,errest = quadgk(f,a,b,atol=1e-14,rtol=1e-14)
(0.1904741736116139, 1.4432899320127035e-15)

We start with the trapezoid formula on \(n=N\) nodes.

N = 20;       # the coarsest formula
n = N;  h = (b-a)/n;
t = h*(0:n);   y = f.(t);

We can now apply weights to get the estimate \(T_f(N)\).

T = [ h*(sum(y[2:n]) + y[1]/2 + y[n+1]/2) ]
err_2nd = Q .- T
1-element Array{Float64,1}:

Now we double to \(n=2N\), but we only need to evaluate \(f\) at every other interior node.

n = 2n;  h = h/2;  t = h*(0:n);
T = [ T; T[1]/2 + h*sum( f.(t[2:2:n]) ) ]
err_2nd = Q .- T
2-element Array{Float64,1}:

As expected for a second-order estimate, the error went down by a factor of about 4. We can repeat the same code to double \(n\) again.

n = 2n;  h = h/2;  t = h*(0:n);
T = [ T; T[2]/2 + h*sum( f.(t[2:2:n]) ) ]
err_2nd = Q .- T
3-element Array{Float64,1}:

Let us now do the first level of extrapolation to get results from Simpson’s formula. We combine the elements T[i] and T[i+1] the same way for \(i=1\) and \(i=2\).

S = [ (4*T[i+1]-T[i])/3 for i in 1:2 ]
err_4th = Q .- S
2-element Array{Float64,1}:

With the two Simpson values \(S_f(N)\) and \(S_f(2N)\) in hand, we can do one more level of extrapolation to get a 6th-order accurate result.

R = (16*S[2] - S[1]) / 15
err_6th = Q .- R

If we consider the computational time to be dominated by evaluations of \(f\), then we have obtained a result with twice as many accurate digits as the best trapezoid result, at virtually no extra cost.

Here are tables of the values and errors. Each time we add a column we have to skip one more estimated value, so we get a triangular structure.

table = (n=[N,2N,4N],trap=T,simp=[missing,S...],sixth=[missing,missing,R])
pretty_table(table,["n","Trapezoid","Simpson","6th order"],backend=:html)
n Trapezoid Simpson 6th order
20 0.19041144993926784 missing missing
40 0.19045880585951175 0.19047459116625973 missing
80 0.1904703513046443 0.19047419978635513 0.1904741736943615
table = (n=[N,2N,4N],trap=Q.-T,simp=[missing,Q.-S...],sixth=[missing,missing,Q-R])
pretty_table(table,["n","2nd order","4th order","6th order"],backend=:html)
n 2nd order 4th order 6th order
20 6.272367234605447e-5 missing missing
40 1.5367752102146692e-5 -4.175546458318191e-7 missing
80 3.822306969603062e-6 -2.617474123556285e-8 -8.274761431614763e-11