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.

1. 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.

2. 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$$

3. 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$$

4. 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.

5. 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}$$

6. 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$$

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"}]

• Wow thanks for your contribution to this 3-year-old question of mine that was perhaps conceived out of naivety. Sorry I didn't notice your reply earlier as I haven't been using StackExchange for the past two years. You really did some good work here, and it will take me quite some days to digest it and try your approach out myself. Nov 15, 2021 at 19:42
• @LeeDavidChungLin Thanks for the reply, hope my answer is satisfying. I got the same question with you and searched for it in StackExchange, found your question, realized there were no answers for such a great question. After I got this answer I decided to post for someone in the future searching for an answer to this question. Edited post embedding pictures inline that just became possible by your upvote :)
– 김경민
Nov 16, 2021 at 11:23
• Your approach inspired something interesting that I'm still working on. Meanwhile, I'd like to keep my question open (not accepting your answer just yet) since the direction your analysis is taking seems to be a bit different from what I'd like to see. To be fair, I got feedbacks offline that my original motivation was naive and wrong in some sense. Wish I could upvote your answer twice. Anyway, thanks again. Nov 20, 2021 at 14:02