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?
$\endgroup$1 Answer
$\begingroup$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=intvalwith 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
