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}:
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[1]/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]/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[2] - S[1]) / 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 |