# 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;

(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/2 + y[n+1]/2) ]
err_2nd = Q .- T

1-element Array{Float64,1}:
6.272367234605447e-5


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/2 + h*sum( f.(t[2:2:n]) ) ]
err_2nd = Q .- T

2-element Array{Float64,1}:
6.272367234605447e-5
1.5367752102146692e-5


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 + h*sum( f.(t[2:2:n]) ) ]
err_2nd = Q .- T

3-element Array{Float64,1}:
6.272367234605447e-5
1.5367752102146692e-5
3.822306969603062e-6


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}:
-4.175546458318191e-7
-2.617474123556285e-8


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 - S) / 15
err_6th = Q .- R

-8.274761431614763e-11


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