Geometric representation of Euler-Maclaurin Summation Formula

Solution 1:

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