Error formula for Composite Trapezoidal Rule

Question

My textbook gives me the error term for the Composite Trapezoidal Rule as this:

$-\frac{b-a}{12}h^2f''(\mu)$, where $\mu \in(a,b)$ and $f \in C^2 [a,b]$

I am using MatLab to produce approximations with the Composite Trapezoidal Rule for $\int_0^{0.99} \frac{1}{\sqrt{1-x^2}}{\rm d}x$ with the intervals $h = 0.01, 0.005, 0.0025, 0.00125, 0.000625$

Below is my table of the approximations by my code and the absolute error for each interval:

 ....h............S(h)...........abs. err...
 0.010000    1.432052842622   0.002795989152
 0.005000    1.429980957924   0.000724104453
 0.002500    1.429439827337   0.000182973867
 0.001250    1.429302728001   0.000045874530
 0.000625    1.429268330467   0.000011476997

Evaluating the error with the error formula, however, gives me a very different number than what my code is spitting out. For example, evaluating the error term for $h = 0.01, a = 0, b = 0.99$, I end up with $0.437161725$. Should my approximation of the error be that off? Am I not using the error term properly?

Answer 1

You should be careful with this expression:

$$ {\rm err} = -\frac{b-a}{12}h^2f''(\mu) \tag{1} $$

The meaning is: there is a point $\mu \in (a,b)$ such that the error is given by this expression. To show this is true I calculate $S(h)$ for various values of $h$ and the absolute error $\epsilon$. I then find the value of $\mu$ guaranteed by Eq. (1), that is, the value of $\mu$ such that ${\rm err} = \epsilon$

Answer 2

Note that the error is in the derivation the sum over the error of all the sub-intervals. For each of these intervals you get $$ \int_{x_k}^{x_{k+1}}f(s)\,ds-(x_{k+1}-x_k)\frac{f(x_{k+1})+f(x_k)}2=\frac{(x_{k+1}-x_k)^3}{12}f''(μ_k) $$ With constant length of the sub-intervals, $x_{k+1}-x_k=h$, the error is thus a Riemann sum for $$ \frac{h^2}{12}\int_a^bf''(s)\,ds=\frac{h^2}{12}(f'(b)-f'(a)) $$ For $f(x)=(1-x^2)^{-1/2}$ we get $f'(x)=x(1-x^2)^{-3/2}$ so that with $[a,b]=[0, 0.99]$ the constant in the error formula is $C=(f'(b)-f'(a))/12=29.38829..$. The experimental table extended to include estimates of this second order error constant is \begin{array}{r|lllllccc} n&h&S(h)&E(h)=S(h)-I&E(h)/h^2&E(h)/h^2-C\\\hline 99 & 0.01 & 1.4320528426221864 & 0.002795989151717 & 27.95989151717082 & -1.4283996452350003\\ 198 & 0.005 & 1.4299809579235057 & 0.000724104453036 & 28.964178121455575 & -0.42411304095024605\\ 396 & 0.0025 & 1.4294398273373048 & 0.000182973866836 & 29.275818693683675 & -0.11247246872214589\\ 792 & 0.00125 & 1.4293027280009247 & 0.000045874530455 & 29.35969949149353 & -0.028591670912291534\\ 1584 & 0.000625 & 1.4292683304674152 & 0.000011476996946 & 29.381112181567914 & -0.007178980837906579\\ 3168 & 0.0003125 & 1.4292597232453161 & 0.000002869774847 & 29.38649443194663 & -0.0017967304591905986\\ \end{array} which shows that indeed the error behaves like $E(h)=29.38829\cdot h^2+O(h^4)$.