How to prove the equation $\int_{a}^{b} f(x) \mathrm{d} x = \frac{b-a}{2}[f(a)+f(b)]-\frac{(b-a)^{3}}{12} f^{\prime \prime}(\eta)$? In numerical analysis, we know that the following equation is true:
$$\int_{a}^{b} f(x) \mathrm{d} x = \frac{b-a}{2}[f(a)+f(b)]-\frac{(b-a)^{3}}{12} f^{\prime \prime}(\eta), \quad \eta \in(a, b)$$
I wonder if there is any profound and concise way to prove this conclusion. It is better to use Taylor series to prove this conclusion.

The link Trapezoid rule shows the error term is bounded but is a little unsatisfactory since it does not prove equality in this particular case.
Simplifying to the interval [0,1], If we assume $f$ is twice continuously differentiable on $[0,1]$ then, using integration by parts we can proceed thus,
\begin{aligned}
\int_0^1 f(x)~dx &= \Big[ \big(x-\tfrac{1}{2}\big) f(x) \Big]_0^1 - \int_0^1\big(x-\tfrac{1}{2}\big) f'(x) ~ dx \\
&=\frac{f(1)+f(0)}{2}-\Big[\big(\tfrac{1}{2}x^2-\tfrac{1}{2}x\big)f'(x)\Big]_0^1+\int_0^1 \big(\tfrac{1}{2}x^2-\tfrac{1}{2}x\big)f''(x)~dx \\
&= \frac{f(1)+f(0)}{2}+\int_0^1 \big(\tfrac{1}{2}x^2-\tfrac{1}{2}x\big)f''(x)~dx 
\end{aligned}
Now, $f''(x)$ is continious and therefore bounded, say $m \leqslant f’’(x) \leqslant M$.  The expression $\tfrac{1}{2}x^2-\tfrac{1}{2}x \leqslant 0 $ throughout $[0,1]$, and therefore
\begin{aligned}
M\int_0^1 \big(\tfrac{1}{2}x^2-\tfrac{1}{2}x \big)~dx \leqslant \int_0^1 \big(\tfrac{1}{2}x^2-\tfrac{1}{2}x \big)f’’(x) ~dx \leqslant m \int_0^1 \big(\tfrac{1}{2}x^2-\tfrac{1}{2}x \big)~dx.
\end{aligned}
where the inequalities can be taken as strict unless $f’’$ is constant.  This becomes,
\begin{aligned}
-\frac{M}{12} \leqslant \int_0^1 \big(\tfrac{1}{2}x^2-\tfrac{1}{2}x\big)f’’(x)~dx \leqslant -\frac{m}{12}
\end{aligned}
Moreover, being continuous, $f’’$ takes every value between $m$ and $M$ so there exists $\xi\in(0,1) ^1$ such that,
\begin{aligned}
\int_0^1 \big(\tfrac{1}{2}x^2-\tfrac{1}{2}x\big)f’’(x)~dx = -\frac{1}{12}f’’(\xi)
\end{aligned}
which completes the proof, that there exists $\xi\in(0,1)$ such that:
\begin{aligned}
\int_0^1 f(x) ~dx = \frac{f’’(0)+f’’(1)}{2} - \frac{1}{12} f''(\xi).
\end{aligned}
The result for the case $[a,b]$ can be derived through a change of variable.
Corrected to refer to $f’’$ consistently

Footnote: if $f’’$ is not constant we can choose $\xi$ in the interior of $[0,1]$ because the inequalities are then strict;  if $f’’$ is constant any point $\xi$ will serve.
