[Math] Richardson’s extrapolation of composite trapezoidal rule

approximation-theoryintegrationnumerical methods

I have applied Richardson's formula to the composite trapezoidal rule, $I_h(f)=\frac{h}{2}(f(a)+\sum_{k=1}^{n-1}f(a+kh)+f(b))$, in an attempt to better approximate the integral $I(f)=\int_0^1 e^{-x}dx$. Here is my work:

Define $N_1(h)=\frac{h}{2}(f(a)+\sum_{k=1}^{n-1}f(a+kh)+f(b))$. Since Composite Trapezoidal rule is of degree of accuracy 1, we want to find $N_2(h)$ to approximate this integral.

$N_2(h)=2N_1(\frac{h}{2})-N_1(h)$
$=\frac{h}{2}[(f(a)+\sum_{k=1}^{n-1}f(a+\frac{kh}{2})+f(b))]-\frac{h}{2}[(f(a)+\sum_{k=1}^{n-1}f(a+kh)+f(b))]$
$=h\sum_{k=1}^{n-1}[f(a+\frac{kh}{2})-f(a+kh)]$.

However, using that the actual value of this integral is $0.63212$, no value of $h$ ranging from $2^{-1}$ to $2^{-8}$ ($n$ from $1$ to $255$, since $h=\frac{b-a}{n}$), results in anything close to this number. In fact, the original composite trapezoidal rule gives much better approximations. I have computed these summations using MatLab. Can anyone help me figure out what I did wrong?

Best Answer

The trapezoidal rule and thus also the composite variant are symmetric. Thus the error order is even, which means that it is $2$. Which implies that the extrapolated value is $$ N_2(h)=\frac13(4N_1(\tfrac h2)-N_1(h)) $$

Next, if $n·h=b-a$, then also $2n·\frac h2=b-a$, so that the explicit formulas are \begin{align} N_1(h)&=h·\left(\tfrac12 f(a)+\sum_{k=1}^{n-1}f(a+kh)+\tfrac12 f(b)\right) \\ N_1(\tfrac h2)&=\tfrac h2·\left(\tfrac12 f(a)+\sum_{k=1}^{2n-1}f(a+k\tfrac h2)+\tfrac12 f(b)\right) \\ N_2(h)&=\tfrac h6·\left(f(a)+4f(a+\tfrac h2)+\sum_{k=1}^{n-1}\left[2f(a+kh)+4f(a+(k+\tfrac12)h)\right]+f(b)\right) \end{align} which is the composite Simpson rule (also symmetric) of error order $4$.

To experiment one could use this python code 

from math import exp

def f(x):
    return exp(-x)

def trapezint(f,a,b,N):
    h=(b-a)/N
    sum=(f(a)+f(b))/2;
    for k in range(1,N):
        sum = sum + f(a+k*h)
    return sum*h

last = 0;
for k in range(10):
    N = 1<<k;
    intval = trapezint(f, 0.0, 1.0, N);
    print "N=%3d, I(N) = %.12f,  N*N*(I(N/2)-I(N))= %f, (4*I(N)-I(N/2))/3 = %.12f" \
    % (N, intval, N*N*(last-intval), (4*intval-last)/3)
    last=intval

with output 

N=  1, I(N) = 0.683939720586,  N*N*(I(N/2)-I(N))= -0.683940, (4*I(N)-I(N/2))/3 = 0.911919627448
N=  2, I(N) = 0.645235190149,  N*N*(I(N/2)-I(N))= 0.154818, (4*I(N)-I(N/2))/3 = 0.632333680004
N=  4, I(N) = 0.635409429028,  N*N*(I(N/2)-I(N))= 0.157212, (4*I(N)-I(N/2))/3 = 0.632134175321
N=  8, I(N) = 0.632943418210,  N*N*(I(N/2)-I(N))= 0.157825, (4*I(N)-I(N/2))/3 = 0.632121414605
N= 16, I(N) = 0.632326313844,  N*N*(I(N/2)-I(N))= 0.157979, (4*I(N)-I(N/2))/3 = 0.632120612389
N= 32, I(N) = 0.632172000094,  N*N*(I(N/2)-I(N))= 0.158017, (4*I(N)-I(N/2))/3 = 0.632120562177
N= 64, I(N) = 0.632133419302,  N*N*(I(N/2)-I(N))= 0.158027, (4*I(N)-I(N/2))/3 = 0.632120559038
N=128, I(N) = 0.632123773957,  N*N*(I(N/2)-I(N))= 0.158029, (4*I(N)-I(N/2))/3 = 0.632120558842
N=256, I(N) = 0.632121362611,  N*N*(I(N/2)-I(N))= 0.158030, (4*I(N)-I(N/2))/3 = 0.632120558829
N=512, I(N) = 0.632120759774,  N*N*(I(N/2)-I(N))= 0.158030, (4*I(N)-I(N/2))/3 = 0.632120558829
Related Question