Derive the error for the Midpoint rule Problem statement
Let $[a,b]$ be divided into $n$ subintervals by $n+1$ points $x_0 = a, x_1 = a+h , x_2 = a+2h , \dots , x_n = a+nh =b$.

Show that
$$
I- I_M = \frac{1}{24}(b-a)h^2 f''(z)
$$
for some $z \in [a,b]$.

Attempt at solution
$$
I_M = h \sum_{i=0}^{n-1} f\left( a + \left( i + \frac12\right)\right) := h\sum_{i=0}^{n-1}f(x_i)
$$
\begin{align}
I -I_M &= \int_a^b f(x) \mathrm{d}x - h\sum_{i=0}^{n-1}f(x_i)\\
&= \sum_{i=0}^{n-1} \left( \int_{x_i}^{x_{i+1}} f(x) \ \mathrm{d}x - hf(x_i)\right)
\end{align}
Here I am stuck. I know I need to use Taylor's theorem in conjunction with Mean-Value Theorem but I don't know how to proceed.
I suspect that for $x\in [x_i ,x_{i+1}]$ the Taylor expansion is
$$
f(x) = f(x_i) + (x-x_i)f'(x_i) + \frac12 (x-x_i)^2 f''(z_x) \quad \text{for } z_x \in [x_i , x_{i+1}]
$$
but then If I substitute this in the above I get
\begin{align}
&= \sum_{i=0}^{n-1}\left( \int_{x_i}^{x_{i+1}} (f(x_i) + (x-x_i)f'(x_i) + \frac12 (x-x_i)^2 f''(z_x))\mathrm{d}x - hf(x_i)\right)\\
\end{align}
and I don't see what I have accomplished
Does anyone know If I am on the right track and if yes how to continue this derivation ?
 A: There is another proof that avoids double integrals and uses the generalized mean value instead. I have learned this technique from reading @LutzLehmann's answers to very similar questions.

Let $h>0$ be given and consider the problem of computing the integral 
$$I = \int_{-h}^h f(x)dx.$$
The midpoint rule takes the form
$$M_h = 2h f(0).$$
We will now obtain the familiar error formula by studying the auxiliary function $g$ given by
$$ g(x) = \int_{-x}^x f(t) dt - 2x f(0).$$
This function is interesting precisely because
$$ g(h) = I - M_h.$$
Let $F$ denote an anti-derivative of $f$. Then
$$ g(x) = F(x) - F(-x) - 2xf(0).$$
It is clear that $g$ is as smooth as $F$ and that $g(0) = 0$. We have $$g'(x) = f(x) + f(-x) - 2f(0).$$ It follows that $$g'(0) = 0.$$ We have $$g''(x) = f'(x) - f'(-x).$$ It follows that $$g''(0) = 0.$$ We conclude that the graph of $g$ is quite flat near $x=0$. To make further progress we make repeated use of the generalized mean value theorem. There exists at least one $x_1$ between $0$ and $x$ such that
$$ \frac{g(x)}{x^3} = \frac{g(x) - g(0)}{x^3 - 0^3} = \frac{g'(x_1)}{3x_1^2}.$$ There exists at least one $x_2$ between $0$ and $x_1$ such that 
$$ \frac{g'(x_1)}{3x_1^2} = \frac{g'(x_1) - g'(0)}{3x_1^2 - 3\cdot 0^2}= \frac{g''(x_2)}{6x_2}.$$ Similarly, there exists at least one $x_3$ between $0$ and $x_2$ such that 
$$ \frac{g''(x_2)}{6x_2} = \frac{g''(x_2)- g''(0)}{6x_2 - 6 \cdot 0} = \frac{g'''(x_3)}{6} = \frac{f''(x_3) + f''(-x_3)}{6}$$
Finally, by the continuity of $f''$ there is at least one $\xi$ between $-x_3$ and $x_3$ such that $$\frac{g(x)}{x^3} = \frac{f''(\xi)}{3}.$$
In particular, there is at least one $\xi$ in $[-h,h]$ such that $$I - M_h = g(h) = \frac{f''(\xi)}{3} h^3.$$
Replacing $h$ with $h/2$ yields the familiar result, i.e,
$$ \int_{-h/2}^{h/2} f(x) dx - h f(0) = \frac{f''(\xi)}{24} h^3.$$
We finish the proof as suggested by @RRL.
A: You can start by obtaining the rule for a single interval, say
\begin{align*}
\int_{x_i}^{x_{i+1}} f(x) dx -hf(c_i) &= \int_{x_i}^{x_{i+1}}(f(x)-f(c_i))dx\\
& = \int_{x_i}^{x_{i+1}} \left(f(c_i)+f'(c_i)(x-c_i) + \frac{f''(\xi_i(x))}{2}(x-c_i)^2\;-\;f(c_i)\right) dx\\
& = \frac{f''(\xi_i)}{2}\int_{x_i}^{x_{i+1}}(x-c_i)^2 dx = \frac{f''(\xi_i) h^3}{24}
\end{align*}
Now, if you take the sum over the $n$ intervals, you get the error term
$$
\sum_{i=1}^n \frac{f''(\xi_i) h^3}{24} = \frac{h^2(b-a)}{24}\left(\frac 1n \sum_{i=1}^n f''(\xi_i)\right) = \frac{h^2(b-a) f''(\xi)}{24}.
$$

Some notes:

*

*I used the mean value theorem for integrals. Since $(x-c_i)^2$ does not change sign, we know that
$$
\int_{x_i}^{x_{i+1}} \frac{f''(\xi_i(x))}{2}(x-c_i)^2 dx = \frac{f''(\xi_i)}{2} \int_{x_i}^{x_{i+1}} (x-c_i)^2 dx
$$


*If we assume the $f\in C^2$, the term $\frac 1n \sum f''(\xi_i)$ is an average of values of $f''$, hence between the minimum and maximum of $f''$. The intermediate value theorem says that it will correspond to $f''(\xi)$, for some $\xi \in (a,b)$.
A: To make any progress we need to consider the smoothness of $f$. Here I assume $f \in C^2([a,b])$. Denote the partition points as $x_j = a + hj$ and the midpoints as $c_j = a + \frac{2j+1}{2}h$ for $j = 0,1,\ldots, n-1$.
We first find the local error for a subinterval $[x_j,x_{j+1}]$  using the Taylor expansion with integral remainder
$$f(x) = f(c_j) +f'(c_j)(x-c_j) + \int_{c_j}^x(x-t)f''(t)\, dt$$
Integrating over $[x_j,x_{j+1}]$ and using $x_{j+1} - x_j = h$, we get
$$\tag{1}\int_{x_j}^{x_{j+1}}f(x) \, dx = hf(c_j)+ f'(c_j)\int_{x_j}^{x_{j+1}}(x-c_j) \, dx + \int_{x_j}^{x_{j+1}} \int_{c_j}^x(x-t)f''(t)\, dt\, dx$$
The second term on the RHS of (1) is
$$f'(c_j)\int_{x_j}^{x_{j+1}}(x-c_j) \, dx =f'(c_j) \left.\frac{(x-c_j)^2}{2} \right|_{x_j}^{x_{j+1}}= \frac{f'(c_j)}{2} \left[(x_{j+1}-c_j)^2- (x_j - c_j)^2 \right]\\ = \frac{f'(c_j)}{2}\left[\left(\frac{h}{2}\right)^2 - \left(-\frac{h}{2}\right)^2\right]= 0$$
Substituting into (1) and rearranging, we get
$$\tag{2}\int_{x_j}^{x_{j+1}}f(x) \, dx - hf(c_j)=  \int_{x_j}^{x_{j+1}} \int_{c_j}^x(x-t)f''(t)\, dt\, dx$$
In order to apply the mean value theorem to extract the second derivative from the integral on the RHS, some manipulation with indicator functions is required to eliminate the variable integration limit.  Note that
$$\int_{x_j}^{x_{j+1}} \int_{c_j}^x(x-t)f''(t)\, dt\, dx  \\= \int_{c_j}^{x_{j+1}} \int_{c_j}^x(x-t)f''(t)\, dt\, dx + \int_{x_j}^{c_j} \int_{x}^{c_j}(t-x)f''(t)\, dt\, dx \\ =  \int_{c_j}^{x_{j+1}} \int_{c_j}^{x_{j+1}}\mathbf{1}_{}(t \leqslant x)(x-t)f''(t)\, dt\, dx + \int_{x_j}^{c_j} \int_{x_j}^{c_j}\mathbf{1}_{}(t \geqslant x)(t-x)f''(t)\, dt\, dx $$
Now we can apply Fubini's theorem to obtain
$$\tag{3}\int_{x_j}^{x_{j+1}} \int_{c_j}^x(x-t)f''(t)\, dt\, dx \\ =  \int_{c_j}^{x_{j+1}} \int_{c_j}^{x_{j+1}}\mathbf{1}_{(t \leqslant x)}(x-t)f''(t)\, dx\, dt + \int_{x_j}^{c_j} \int_{x_j}^{c_j}\mathbf{1}_{(t \geqslant x)}(t-x)f''(t)\, dx\, dt \\ =  \int_{c_j}^{x_{j+1}} f''(t)\left(\int_{t}^{x_{j+1}}(x-t)\, dx\right)\, dt + \int_{x_j}^{c_j} f''(t)\left(\int_{x_j}^{t}(t-x)\, dx\right)\, dt\\ = \int_{c_j}^{x_{j+1}} f''(t)\frac{(x_{j+1}-t)^2}{2}\, dt + \int_{x_j}^{c_j} f''(t)\frac{(t-x_j)^2}{2}\, dt$$
By the mean value theorem for integrals there are points $z_j^+ \in [c_j,x_{j+1}]$ and $z_j^- \in [x_j,c_j]$ such that
$$\tag{4}\int_{c_j}^{x_{j+1}} f''(t)\frac{(x_{j+1}-t)^2}{2}\, dt + \int_{x_j}^{c_j} f''(t)\frac{(t-x_j)^2}{2}\, dt\\= f''(z_j^+)\int_{c_j}^{x_{j+1}} \frac{(x_{j+1}-t)^2}{2}\, dt +  f''(z_j^-)\int_{x_j}^{c_j} \frac{(t-x_j)^2}{2}\, dt\\ = f''(z_j^+)\frac{h^3}{48} + f''(z_j^-) \frac{h^3}{48} $$
Substituting into (2) with (3) and (4) yields
$$\int_{x_j}^{x_{j+1}}f(x) \, dx - hf(c_j)=\frac{h^3}{48}\left(f''(z_j^+)+f''(z_j^-)\right)$$
Summing from $j=0$ to $j = n-1$, we get
$$I-I_M = \int_a^b f(x) \, dx - h\sum_{j=0}^{n-1}f(c_j)= \frac{h^3}{48}\sum_{j=0}^{n-1}\left(f''(z_j^+)+f''(z_j^-)\right)\\ = \frac{2nh^3}{48}\frac{1}{n}\sum_{j=0}^{n-1}\frac{f''(z_j^+)+f''(z_j^-)}{2},$$
and since $nh = (b-a)$,
$$\tag{5}I - I_M = \frac{(b-a)h^2}{24} \sum_{j=0}^{n-1}\frac{f''(z_j^+)+f''(z_j^-)}{2n}$$
I'll leave it to you to consider if there  exists  $z \in [a,b]$ such that the second derivative at this point equals the weighted average of second derivatives, i.e.,
$$f''(z) =  \sum_{j=0}^{n-1}\frac{f''(z_j^+)+f''(z_j^-)}{2n},$$
and, hence,
$$I - I_M = \frac{(b-a)h^2}{24}f''(z) $$
Nevertheless, if the second derivative is bounded as $|f''(x)| \leqslant \|f''||_{\infty}$ for all $x \in [a,b]$, then we obtain from (5) the well-known midpoint rule error bound
$$|I - I_M| \leqslant \frac{\|f''\|_{\infty}(b-a)h^2}{24} $$
