Geometric representation of Euler-Maclaurin Summation Formula In Tom Apostol's expository article (here's a free link), upon seeing the figure below (or this from the Wolfram project) I was expecting more diagrams to come to continue the error decomposition of the shaded regions in the shape of "curved triangular pieces".
For each horizontal unit interval, it seems that the hypotenuse of an implicit triangle (within shaded region) is the linear correction. Then, there's a "cap" that fills in the gap between the hypotenuse and the actual curve. In this figure, each cap is attached to the implicit hypotenuse from below due to the curve being convex.

It was a disappointment not seeing further decomposition of the "cap" into regions described by parabolic, cubic, then quartic curves etc. Somehow none of the textbooks and other materials I've read have such diagrams for higher orders.
Question Statement
What is the correct analysis to demonstrate the 'geometry' of the successive orders of the Euler-Maclaurin formula?
Is there any existing source with such diagrams (visualization) similar to the above regarding higher order?
General Remarks
I'm actually not sure whether my proposal is a valid idea (that such a demonstration is possible). Below is a more detailed description (repeating things said above) of what I mean by the tentative decomposition of pieces as the terms of Euler-Maclaurin formula.
Basically it is approximating the definite integral  $F(b)-F(a)$ where $b>a$, unfolding in the direction "opposite" to Euler-Maclaurin formula itself.
$$
\int_{a}^b f(x) \,\mathrm{d}x \approx F_0 + F_1 + F_2 + \cdots
$$
The approximation is done by columns of unit width (summation of  $f(k)$, discretely sampled points) plus corrections to get close to the curve (via derivatives of end points).
$$
\int_{a}^b f(x) \,\mathrm{d}x \approx \overbrace{ \underbrace{ \sum_{k = a}^b f(k) }_{ \textbf{unit columns} }  - \frac{f(a)+f(b)}2 }^{ \textbf{unit columns centered} } - \underbrace{ \frac{f'(b) - f'(a)}{12} - 0 }_{ \textstyle\genfrac{}{}{0pt}{}{\textbf{net result of} }{ \textbf{triangular & parabolic cap} } } - \overbrace{ \frac{f'''(b) - f'''(a)}{720} - 0}^{ \textbf{cubic & guartic cap} }  - \cdots$$
The correct analysis would have to account for both the emergence of Bernoulli numbers and the individual orders of correction ... and the remainder if possible.
Elaboration on the Ideas
For simplicity, consider $(a,b) \in \mathbb{N}^2$ and let the summation be in unit steps as usual, then pretend that things will scale nicely.

*

*The zero-th order approximation for the integral is the columns of unit width $$F_0 = \sum_{k=a}^b f(k)$$ Note that each column of height $f(k)$ is left-aligned. For example, the first column for $k = a$ occupies $x \in [a, a+1]$. This means the last column at $k = b$ is entirely outside of the range of integration.


*The $2$nd term (still zero-th order) shifts all the columns by $\frac12$ to make them centered. That is, each column of height $f(k)$ now resides at $x \in [a-\frac12, a+\frac12]$. The last column is now half-outside, and so is the first column. Thus we remove half of each of them. $$ F_1 = -\frac{f(a) + f(b)}2$$


*The $3$rd term ($1$st order, linear) correction is supposed to shave off a triangular top from each column (note that the outermost are half-columns at $k = b$ and $k = a$). The height of each triangle is $f'(k)$ and the width is unity. $\color{brown}{\textit{Somehow}}$ there should be a telescoping of terms (combine or cancellation) with the $3$rd order correction and we end up with $$F_2 = - \frac16 \frac{f'(b) - f'(a)}2$$


*The $4$th term ($2$nd order, parabolic) correction is supposed to be something involving both $f''(k)$ and $f'(k)$. It should be a parabolic piece on top of the triangular (linear) correction of the previous order, similar in spirit** to what one sees in the quadrature of parabola. Note that there's a fixed ratio of $\frac13$ (or $\frac23$) between the parabolic area and the rectangle it spans. $\color{brown}{\textit{Somehow}}$ there should be telescoping with the $2$nd as well as the $4$th order so that $$F_3 = 0~.$$
**Here one uses parabolic pieces to fill the gaps for one iteration then uses cubic pieces for the next iteration and so on, while in ancient quadrature one keeps using triangular pieces that are linear.


*The $5$th term ($3$rd order, cubic) is supposed to some kind of lune (moon crest) on top of the parabolic correction, and $\color{brown}{\textit{somehow}}$ after cancellation it will end up being $$F_4 = -\frac{f'''(b) - f'''(a)}{720}$$


*So on and so forth ....
In the end, one shall keep the summation $F_0$ on one side of the equation by itself, moving the terms $F_1$, $F_2$, $F_3$ etc to the side of the integral, to have the Euler-Maclaurin summation formula.
$$F_0 \approx \int_{a}^b f(x) \,\mathrm{d}x - F_1 - F_2 - \cdots$$
 A: Let’s think were going to approximate $\int_0^1 f(x) \, dx $ in a certain way, then the integral we want to approximate can be expressed as below, then can be approximated.
$$
∫_a^b f(x) \, dx = \sum_{k=a}^{b-1} { \int_k^{k+1} f(x) \, dx} = \sum_{k=a}^{b-1} {\int_0^1 f(x+k) \, dx} 
$$
Now, we’re going to set a polynomial $p_n (x)$ which has degree lower or equal than $2n+1$ that satisfies
$$
p_n(0) = f(0) , \; p_n(1) = f(1) , \; p^{(2 j - 1)}_n(0) = f^{(2j-1)}(0) , \; p^{(2j-1)}_n(1) = f^{(2j-1)}(1) \quad (j = 1,2,\cdots n)
$$
When $n=0$, which means there are no derivative conditions, $p_0(x)$ simply becomes a line that connects $(0,f(0))$ and $(1,f(1))$. Roughly speaking, as $n$ increases, even though we’re not using derivatives of even orders, $p_n(x)$ would look more similar with $f(x)$ in $[0,1]$ since $p_n(x)$ is unique and $f(x)$ can be considered as a infinite degree polynomial by taylor expansion, that satisfies all the condition $p_n(x)$ has to satisfy. Being a bit more rigorous knowing the result, if the remainder term of the Euler-Maclaurin formula for $f(x)$ in $[0,1]$ goes to $0$, $p_n(x)$ would converge to $f(x)$.
Beyond topic, the unique existence of this polynomial $p_n(x)$ can be prooved. There are $2n+2$ conditions, and the number of coefficients in a polynomial of degree lower or equal to $2n+1$ is $2n+2$, and the equations for the coefficients would be a linear problem. $(2n+2) \times (2n+2)$ matrix of this linear problem can't have linearly dependent rows because then it would generate a general relation between $g(0),g(1),g'(0),g'(1),\cdots g^{(2n-1)}(0),g^{(2n-1)}(1)$ that works for all $g(x)$, polynomial of degree lower or equal than $2n+1$.
Now we simply approximate as below, expressing error as $R_n$.
$$
\int_0^1 f(x) \, dx = \int_0^1 p_n (x) \, dx + R_n, \quad R_n =: \int_0^1 \left(f(x) - p_n (x) \right) \, dx 
$$
The proof that this approximation gives the Euler-Maclaurin formula can be prooved by itself, the Euler-Maclaurin formula. Since $\deg( p_n (x)) \leq 2n+1 $, derivatives of order larger than $2n+1$ will be $0$. Also, derivate of order $2n+1$ would be a constant. Therefore,
$$
\begin{split}
\int_0^1 p_n (x) \, dx &= \frac{p_n(0)+p_n(1)}{2} - \sum_{j=1}^{n}{ \frac{B_{2j}}{(2j)!} \left( p^{(2j-1)}_n(1)-p^{(2j-1)}_n(0) \right)} \\ &= \frac{f(0)+f(1)}{2} - \sum_{j=1}^{n}{ \frac{B_{2j}}{(2j)!} \left( f^{(2j-1)}(1)-f^{(2j-1)}(0) \right)}
\end{split}
$$
Now, let’s use $f(x+k)$ instead of $f(x)$, and say the approximated polynomial is $p_{n,k} (x)$, and the error term is $R_{n,k}$. After telescoping, we can get the Euler-Maclaurin formula written for approximating the integral.
$$
\begin{split}
\int_a^b f(x) \, dx &= \sum_{k=a}^{b-1} \left( \int_0^1 p_{n,k}(x) \, dx + R_{n,k} \right) \\ 
&= \sum_{k=a}^b f(k) - \frac{f(a)+f(b)}{2} - \sum_{j=1}^{n}{ \frac{B_{2j}}{(2j)!} \left( f^{(2j-1)}(b)-f^{(2j-1)}(a) \right)} + \sum_{k=a}^{b-1} R_{n,k}
\end{split}
$$
You can easily see how this can be changed to a geometrical representation of the Euler-Maclurin formula. Simpily draw a function $P_n(x)$, which is $p_{n,k}(x-k)$ for $x \in [k,k+1]$ for $k=a,a+1,\cdots b-1 $, and the Euler-Maclurin formula now means that $\int_a^b f(x) \, dx$ can be approximated to $\int_a^b P_n(x) dx$.
The reason for not using the derivative of even orders can be explained intuitively. When we are summing, telescoping should work and cancel out terms so only the derivatives at $x=a$ or $x=b$ are used. But when you think about the relation of $\int_0^1 f(x) \, dx$ and the value $f^{(k)} (0)$, $f^{(k)}(1)$, you can see that when $k$ is even, $f^{(k)}(0)$ or $f^{(k)}(1)$ increasing would lead $\int_0^1 f(x)\, dx$ increase, and it won’t be able to do telescoping. On the other hand, when $k$ is odd, $f^{(k)}(0)$ increasing or $f^{(k)}(1)$ decreasing would lead $\int_0^1 f(x)\, dx$ increase and the possiblity of telescoping is alive.
This is a picture of visualization for $f(x)=\cos(5x)$. You can see that as $n$ increases, $P_n(x)$ gets close to $f(x)$.
$P_0(x), P_1(x), P_2(x), P_3(x)$ for $f(x)=\cos(5x)$" />
And this a picture that visualizes the error for $n=4$, positive area is red and negative area is blue.
$n=4$, $f(x)=\cos(5x)$" />
As I mentioned earlier, $P_n(x)$ can diverge if the remainder term of the Euler-Maclaurin formula diverges. You can check this by an example, $f(x)=\cos(7x)$. The graph is in below.
$P_0(x), P_1(x), P_2(x), P_3(x)$ for $f(x)=\cos(7x)$" />
These are the Mathematica codes that I used for drawing these graphs. You can try it for yourself with different $f(x),a,b,n$ if you have Mathematica.
f[x_] := Cos[5 x];
a = 0; b = 5;

getsolp[n0_, k0_, x_] := Module[{n = n0, k = k0},
  g[zeta_] := Sum[c[i] * zeta^i, {i, 1, 2 n + 1}] + c[0];
  cond = Join[{g[k] == f[k]}, {g[k + 1] == f[k + 1]}, 
    Table[(D[g[zeta], {zeta, 2 j - 1}] /. 
        zeta -> k) == (D[f[zeta], {zeta, 2 j - 1}] /. zeta -> k), {j, 
      1, n}], Table[(D[g[zeta], {zeta, 2 j - 1}] /. 
        zeta -> k + 1) == (D[f[zeta], {zeta, 2 j - 1}] /. 
        zeta -> k + 1), {j, 1, n}]];
  sol = NSolve[cond, Table[c[i], {i, 0, 2 n + 1}]];
  g[x] /. sol
  ]

P[n_, x_] := getsolp[n, Floor[x], x]

Plot[{f[x], P[0, x], P[1, x], P[2, x], P[3, x]}, {x, a, b}, 
 PlotLegends -> {"f(x)", Subscript[P, 0][x], Subscript[P, 1][x], 
   Subscript[P, 2][x], Subscript[P, 3][x]}, AxesLabel -> {"x", "y"}]

Plot[{f[x], P[4, x]}, {x, a, b}, 
 Filling -> {1 -> {{2}, {{Opacity[0.2], Red}, {Opacity[0.2], Blue}}}},
  PlotLegends -> {"f(x)", Subscript[P, 4][x]}, 
 AxesLabel -> {"x", "y"}]

