Proving error bound for Simpson's rule

Question

The Simpson's rule can be stated as follows:

$$\int\limits_{x_0}^{x_2}f(x)dx\approx \frac{h}3\left[f(x_0)+4f(x_1)+f(x_2)\right]$$

The way I'm trying to find the error bound for the Simpson's rule is as follows:

• Taylor-expand $f(x)$ about $x_0$, $x_1$ and $x_2$ up to and including the 4th derivative:

$$f(x) = f(x_i)+(x-x_i)f'(x)+(x-x_i)^2\frac{f''(x_i)}2 +(x-x_i)^3\frac{f^{(3)}(x_i)}6+(x-x_i)^4\frac{f^{(4)}(\xi_i)}{24}$$

• Add the three expansions of $f(x)$ as in the Simpson's rule and rearrange the terms and coefficients accordingly, so that $f(x)$ is LHS, and the rest is RHS.

• Take the integrals of both sides, namely integrate the RHS to find the error. And here's the part I'm having an issue with. Upon integrating I'm getting \

$$\int\limits_{x_0}^{x_2} f(x) dx = \frac{h}3 (f(x_0)+4f(x_1)+f(x_2)) + \frac12 h^2f'(x_0)-\frac12 f'(x_2)+\frac43 h^3 f''(x_0)$$ $$+\frac43 h^3 f''(x_1)+\frac43 h^3 f''(x_2)+\frac23 h^4f^{(3)}(x_0)-\frac23 h^4 f^{(3)}(x_2)$$ $$+\frac{4}{15}h^5 f^{(5)}(\xi_0)+\frac{1}{15}h^5 f^{(5)}(\xi_1) + \frac{4}{15}h^5 f^{(5)}(\xi_2)$$

It looks like I might be on the right track, but my derivative terms preceding the fourth derivative don't cancel out. What can I do to fix this?

Update: MVT might possibly help here, I'll try it out in a few hours and will report of the results.

Answer 1

Summary answer: No, this isn't gonna fly.

Long answer: There is a conceptual problem invalidating your first step. You seem to be saying that $$\begin{align}\int_{x_0}^{x_2}\frac1{24}(x-x_0)^4f^{(4)}(\xi_0)dx&=\frac1{24}f^{(4)}(\xi_0)\int_{x_0}^{x_2}(x-x_0)^4dx=\left.\frac1{120}f^{(4)}(\xi_0)(x-x_0)^5\right|_{x_0}^{x_2}\\ &=\frac1{120}f^{(4)}(\xi_0)(x_2-x_0)^5=\frac4{15}h^5f^{(4)}(\xi_0).\end{align}$$ But that's not really the case because the parameter $\xi_0$ is a function of $x$ so it can't just be pulled out of the integral like this.

Now, we can fix this problem by taking the Taylor series for the primitive $F^{\prime}(x)=f(x)$ at the $3$ points and taking linear combinations: $$\begin{align}\int_{x_0}^{x_2}f(x)dx&=a\left(F(x_0)+2hf(x_0)+2h^2f^{\prime}(x_0)+\frac43h^3f^{\prime\prime}(x_0)+\frac23h^4f^{\prime\prime\prime}(x_0)+\frac4{15}h^5f^{(4)}(\xi_1)\right)\\ &+b\left(F(x_1)+hf(x_1)+\frac12h^2f^{\prime}(x_1)+\frac16h^3f^{\prime\prime}(x_1)+\frac1{24}h^4f^{\prime\prime\prime}(x_1)+\frac1{120}h^5f^{(4)}(\xi_2)\right)\\ &+cF(x_2)\\ &-dF(x_0)\\ &-e\left(F(x_1)-hf(x_1)+\frac12h^2f^{\prime}(x_1)-\frac16h^3f^{\prime\prime}(x_1)+\frac1{24}h^4f^{\prime\prime\prime}(x_1)-\frac1{120}h^5f^{(4)}(\xi_3)\right)\\ &-g\left(F(x_2)-2hf(x_2)+2h^2f^{\prime}(x_2)-\frac43h^3f^{\prime\prime}(x_2)+\frac23h^4f^{\prime\prime\prime}(x_2)-\frac4{15}h^5f^{(4)}(\xi_4)\right).\end{align}$$ This will be valid provided $a+b+c=d+e+g=1$ because in that case we end up with a fancy expression for $F(x_2)-F(x_0)$. But we want to zap the terms involving the primitive and want the zero-order terms to amount to Simpson's rule so we require $$\begin{align}a-d&=0\\ b-e&=0\\ c-g&=0\\ 2a&=1/3\\ b+e&=4/3\\ 2g&=1/3.\end{align}$$ We can solve to get $$\begin{align}a&=d=1/6\\ b&=e=2/3\\ c&=g=1/6.\end{align}$$ And we observe that in fact $a+b+c=d+e+g=1$ as required. Using these results we make forward progress to $$\begin{align}\int_{x_0}^{x_2}f(x)dx&=\frac h3\left(f(x_0)+4f(x_1)+f(x_2)\right)+\frac13h^2\left(f^{\prime}(x_0)-f^{\prime}(x_2)\right)\\&+\frac29h^3\left(f^{\prime\prime}(x_0)+f^{\prime\prime}(x_1)+f^{\prime\prime}(x_2)\right)\\ &+\frac19h^4\left(f^{\prime\prime\prime}(x_0)-f^{\prime\prime\prime}(x_2)\right)+\frac1{180}h^5\left(8f^{(4)}(\xi_1)+f^{(4)}(\xi_2)+f^{(4)}(\xi_3)+8f^{(4)}(\xi_4)\right).\end{align}$$ Comparing this to your expression we can see that you have made many calculational errors and typos. But the method that allowed us to get to here also provides a hint as to the next step to advance our derivation. We are going to take the Taylor series for $f^{\prime}(x)$ at the $3$ points and take more linear combinations: $$\begin{align}f^{\prime}(x_0)-f^{\prime}(x_2)&=if^{\prime}(x_0)\\ &+j\left(f^{\prime}(x_1)-hf^{\prime\prime}(x_1)+\frac12h^2f^{\prime\prime\prime}(x_1)-\frac16h^3f^{(4)}(\xi_5)\right)\\ &+k\left(f^{\prime}(x_2)-2hf^{\prime\prime}(x_2)+2h^2f^{\prime\prime\prime}(x_2)-\frac43h^3f^{(4)}(\xi_6)\right)\\ &-\ell\left(f^{\prime}(x_0)+2hf^{\prime\prime}(x_0)+2h^2f^{\prime\prime\prime}(x_0)+\frac43h^3f^{(4)}(\xi_7)\right)\\ &-m\left(f^{\prime}(x_1)+hf^{\prime\prime}(x_1)+\frac12h^2f^{\prime\prime\prime}(x_1)+\frac16h^3f^{(4)}(\xi_8)\right)\\ &-nf^{\prime}(x_2).\end{align}$$ Again this will be an identity if $i+j+k=\ell+m+n=1$ and we want to zap the terms involving the first and second derivatives in our current expression for the integral so we must solve $$\begin{align}i-\ell&=0\\ j-m&=0\\ k-n&=0\\ -2\ell/3+2/9&=0\\ (-j-m)/3+2/9&=0\\ -2k/3+2/9&=0.\end{align}$$ And our solution $$\begin{align}i&=\ell=1/3\\ j&=m=1/3\\ k&=n=1/3.\end{align}$$ Does indeed satisfy $i+j+k=\ell+m+n=1$ so we can write down our new expression for the integral $$\begin{align}\int_{x_0}^{x_2}f(x)dx&=\frac h3\left(f(x_0)+4f(x_1)+f(x_2)\right)+\frac19h^4\left(f^{\prime\prime\prime}(x_2)-f^{\prime\prime\prime}(x_0)\right)\\ &+\frac1{180}h^5\left(8f^{(4)}(\xi_1)+f^{(4)}(\xi_2)+f^{(4)}(\xi_3)+8f^{(4)}(\xi_4)\right)\\ &-\frac1{54}h^5\left(f^{(4)}(\xi_5)+8f^{(4)}(\xi_6)+8f^{(4)}(\xi_7)+f^{(4)}(\xi_8)\right).\end{align}$$ At this point it should be clear that we want to expand the third derivative as a Taylor series about the $3$ points and take a linear combination to get $$\begin{align}f^{\prime\prime\prime}(x_2)-f^{\prime\prime\prime}(x_0)&=p\left(f^{\prime\prime\prime}(x_0)+2hf^{(4)}(\xi_9)\right)\\ &+q\left(f^{\prime\prime\prime}(x_1)+hf^{(4)}(\xi_{10})\right)\\ &+rf^{\prime\prime\prime}(x_2)\\ &-sf^{\prime\prime\prime}(x_0)\\ &-t\left(f^{\prime\prime\prime}(x_1)-hf^{(4)}(\xi_{11})\right)\\ &-u\left(f^{\prime\prime\prime}(x_2)-2hf^{(4)}(\xi_{12})\right).\end{align}$$ We require $p+q+r=s+t+u=1$ for this to be an identity and also $p-s=q-t=r-u=0$ to zap the third derivative terms. The system is underdetermined so we must choose some solution such as $$\begin{align}p&=s=1/2\\ q&=t=0\\ r&=u=1/2\end{align}$$ Now we have eliminated the lower derivative terms and arrived at $$\begin{align}\int_{x_0}^{x_2}f(x)dx&=\frac h3\left(f(x_0)+4f(x_1)+f(x_2)\right)\\ &+\frac1{180}h^5\left(8f^{(4)}(\xi_1)+f^{(4)}(\xi_2)+f^{(4)}(\xi_3)+8f^{(4)}(\xi_4)\right)\\ &-\frac1{54}h^5\left(f^{(4)}(\xi_5)+8f^{(4)}(\xi_6)+8f^{(4)}(\xi_7)+f^{(4)}(\xi_8)\right)\\ &+\frac19h^5\left(f^{(4)}(\xi_9)+f^{(4)}(\xi_{12})\right).\end{align}$$ Now if all those fourth derivatives could be considered to be $f^{(4)}(\xi)$, then we would get the right answer
$$\int_{x_0}^{x_2}f(x)dx=\frac h3\left(f(x_0)+4f(x_1)+f(x_2)\right)-\frac1{90}h^5f^{(4)}(\xi).$$ But that isn't valid and I provided a counterexample just to show how tricky this is. The problem with the current approach is that after the first step the error of $\frac1{10}h^5f^{(4)}(\xi)$ is a factor of $9$ too large and has the wrong sign. We can't easily make this smaller.

I have posted a straightforward proof that adapts the proof of the error term for Taylor's series rather than its value so it doesn't run into issues of trying to justify canceling out large error terms to arrive at a smaller error.

Answer 2

There is a nicely conceptual proof on the order of the error term in Simpson's Rule involving Hermite interpolating polynomials. Let $[a,b]$ be a closed interval in $\mathbf R$ with the midpoint $m,$ let $h=(b-a)/2,$ and let $f \in C^4([a,b]).$ Consider the third Hermite polynomial $H(x)$ which meets the data \begin{align*} & f(a), \\ & f(m), f'(m), \\ & f(b). \end{align*} The polynomial $H(x)$ is of degree at most three, and Simpson's rule is exact for polynomials of degree at most three, whence we get $$ \int_a^b f(x) \approx \int_a^b H(x)\,dx = \frac h3 (f(a)+4f(m)+f(b)). $$ To complete the proof, one can use the following result on the error term for the approximation of $f(x)$ by $H(x)$ stating that for every $x \in [a,b],$ $$ f(x) = H(x) + \frac{f^{(4)}(\xi(x))}{4!}(x-a)(x-m)^2 (x-b), $$ where $\xi =\xi(x) \in [a,b];$ we'll need the Weighted Mean Value theorem for integrals and the formula $$ \int_a^b (x-a) (x-m)^2 (x-b)\,dx = -\frac{(b-a)^5}{120}= -\frac{4h^5}{15}. $$