Proving error bound of Simpson's rule

I feel like writing about some techniques I found online on how people prove the error bound of Simpson's rule for numerical integration.

The Simpson's rule:

\( \int_a^b f(x) dx \) can be approximated as \(\frac{b-a}{6}(f(a) + 4 f(\frac{a+b}{2}) + f(b))\).

You can find more info about this rule on Wikipedia.

The topic of today is: how do you prove the error bound for the Simpson's rule, which is given by:

\[ -\frac{1}{90} \left( \frac{a-b}{2} \right)^5 f^{(4)}(c) \]

for some constant \(c\).

And I think this is the best proof I have found out there, which might probably be the standard treatment of proving such error bounds: https://math.stackexchange.com/questions/1759326/proving-error-bound-on-simpsons-rule-numerical-integration

Since the answer there is quite succinct, here I'm going to expand it for my own understanding and future reference.

Also, something that I found quite interesting, turns out MVT is quite useful in coming up with error bounds such as this. I want to write it down here:

Mean value theorem (this is not the most accurate version):

If \(f\) is continuous and differentiable function, then there exists \(c \in (a, b)\) such that \( \frac{f(b)-f(a)}{b-a} = f'(c) \).

Also there's this "extended" mean value theorem, also known as the "Cauchy's mean value theorem" which roughly says:

If \(f, g\) are continuous and differentiable, then there exists \( c \in (a, b) \) such that \( \frac{f(b)-f(a)}{g(b)-g(a)} = \frac{f'(c)}{g'(c)} \).

Both can be proved by Rolle's theorem, which states that:

If \(f\) is continuous in \([a,b]\) and differentiable in \((a,b)\), then if \(f(a) = f(b)\), there must exists \(c \in (a,b)\) such that \(f'(c) = 0\).

The proofs for all of those are available in Wikipedia.

Ok let's dive into the proof to the error bound. The technique described here can also be applied to proving the rectangle rule and trapezoid rule, I think (I have tried it on the rectangle rule, I haven't tried it on the trapezoid rule).

What do we want?

We want to establish the bound of the following:

\[ \int_{a}^{b} f(x) dx - \frac{b-a}{6}( f(a) + 4f((a+b)/2) + f(b) ) \].

Simplify upper and lower bound

Instead of dealing with \(a\) and \(b\), the standard treatment I think is to rewrite it as integration from \(-t\) to \(t\). Let me show you how we can do that.

Starting with

\[ \int_a^b f(x) dx \]

Rewrite \( x = y + \frac{a+b}{2} \) to obtain:

\[ \int_{(a-b)/2}^{(b-a)/2} f\left(y + \frac{a+b}{2}\right) dy \]

Then, define \(g(y) = f(y + \frac{a+b}{2}) \) to obtain:

\[ \int_{(a-b)/2}^{(b-a)/2} g(y) dy \]

Finally define \(t = \frac{b-a}{2} \), then we basically are looking at:

\[ \int_{-t}^{t} g(y) dy \].

Also, 

\[\begin{align} \frac{b-a}{6}(f(a) + 4 f((a+b)/2) + f(b)) & = \frac{b-a}{6}(g((a-b)/2) + 4g(0) + g((b-a)/2)) \\ &= \frac{t}{3}(g(-t) + 4g(0) + g(t)) \end{align}\].

So indeed what we are interested in is the following:

\[ I(t) = \int_{-t}^t g(x) dx - \frac{t}{3}(g(-t) + 4g(0) + g(t)) \]

By the fundamental theorem of calculus, there exists \(G(x)\) such that \(G'(x) = g(x)\) which means \( \int_{-t}^{t} g(x) dx = G(t) - G(-t) \). Hence we can finally derive the final form of \(I(t)\):

\[ I(t) = G(t) - G(-t) - \frac{t}{3}(g(-t) + 4g(0) + g(t)) \]

Note that

\[\begin{align} I'(t) &= g(t) + g(-t) - \frac{1}{3}(g(t) + g(-t) + 4g(0)) - \frac{t}{3}(g'(t) - g'(-t))) \\  I''(t) &= \frac{1}{3}(g'(t) - g'(-t)) - \frac{t}{3}(g''(t) + g''(-t)) \\  I'''(t) &= -\frac{t}{3}(g'''(t)-g'''(-t) ) \end{align}\]

And more importantly, that \(I(0) = I'(0) = I''(0) = I'''(0) = 0\).

MVT Chain

This is the meat of the technique, which is to apply Cauchy's MVT several time. Let me show you in detail the argument for the first iteration, and then I'll show you the chain.

Pick \(h(x) = x^5\) which is continuous and differentiable in the range. Then by Cauchy's MVT we have \(t_1 \in (0, t) \) such that

\[ \frac{I(t) - I(0)}{h(t) - h(0)} = \frac{I'(t_1)}{h'(t_1)} = \frac{I'(t_1)}{5{t_1}^4}\].

But since \(I(0) = 0\) and \(h(0)\), really what we have is:

\[ \frac{I(t)}{t^5} = \frac{I'(t_1)}{5{t_1}^4} \].

Repeat this argument in a chain, we have:

\[ \frac{I(t)}{t^5} =  \frac{I'(t_1)}{5{t_1}^4} =  \frac{I''(t_2)}{20{t_2}^3} =  \frac{I'''(t_3)}{60{t_3}^2} \].

Where \(t > t_1 > t_2 > t_3 > 0 \).

Also notice that

\[\begin{align} \frac{I'''(t_3)}{60{t_3}^2} &= -\frac{1}{90}\frac{g'''(t_3) - g'''(-t_3)}{{t_3} - (-t_3)} \\    &= -\frac{1}{90} g^{(4)}(t_4) \end{align}\]

by MVT where \( t_4 \in (-t_3, t_3)\).

Hence we have

\[ I(t) = -\frac{1}{90}g^{(4)}(t_4)t^5 \].

Substituting \(t = (b-a)/2\) and \(f\) back, we get the desired results. \( \square\)