Café Math : Euler-MacLaurin Formula

Hi, today I want to talk about one of my favorite formula in Mathematics, called the Euler-MacLaurin Formula. $$ \frac{1}{2}\left[\, Q(n) + Q(0) \,\right] + \sum_{k = 1}^{n-1} Q(k) = \int_0^n Q(x)\,dx + R_n, $$ with, $$ R_n = \sum_{k = 1}^\infty \frac{B_{2k}}{(2k)!} \,\left[\, Q^{(2k-1)}(n) - Q^{(2k-1)}(0) \,\right]. $$ It can be interpreted as a series expression of the reminder term $R_n$ of the $n$-step integration using trapezoid method.

The formula is exact for polynomials $Q$ of any degree. Giving it meaning for more general functions other than polynomials involves tricky convergence issues. Sometimes, it doesn't even converge but that's where the fun begins.

There are a lot of interesting things to be said about the Taylor formula,
$$
f(x + h) = \sum_{n \ge 0} f^{(n)}(x)\,\frac{h^n}{n!}.
$$
Of course it gives useful approximations of the function $f$ by polynomials in the vicinity of $x$ as we all learned, but can we go beyond that? The answer is yes... way beyond. First, it gives the exponential generating series of the successive derivative functions of $f$,
$$
f^{(n)}(x) \overset{def.}{=} \frac{d^n}{dx^n}f(x).
$$
Can we think of anything else ? Well yes, and the following interpretation is a gem of Mathematics.
Let $D$ be the differentiation operator $\frac{d}{dx}$, then rewriting the taylor formula as,
$$
e^{hD} f(x) = f(x+h).
$$
gives a way to express the *shift operator*, that sends the function $f(x)$ to the shifted version $f(x + h)$, as the exponential of the derivation $hD$.

Digression. It is hard to overestimate the importance and value of that simple formula. It is one starting point of umbral calculus and what is known as Sheffer sequences, Appel sequences etc... We'll talk of all that, but not today. For the moment let's just put a name on the ingredients of this formula to point out the connection with the relevant part of differential geometry.

- We have the one parameter group $G$ of translations by $h$ acting on functions.
- The Lie algebra $\mathfrak{L}_G$ of this group, spanned as a vector space by the derivation $D = \frac{d}{dx}$.
- The exponential transform,
$$
\exp : \mathfrak{L}_G \rightarrow G,
$$
that maps the
*infinitesimal generator*$hD \in \mathfrak{L}_G$ to the corresponding group element which is the translation by $h$.

The Bernoulli numbers $B_n$ are the coefficients of the following exponential generating series, $$ \frac{z}{e^z- 1} = \sum_{k \ge 0} B_k \,\frac{z^k}{k!} . $$ As the following function, $$ \frac{z}{2} + \frac{z}{e^z - 1} = \frac{z}{2}\,\frac{e^z+1}{e^z-1} = -\frac{z}{2}\,\frac{e^{-z}+1}{e^{-z}-1}, $$ is an even function, every odd numbered Bernoulli number is zero apart from $B_1 = -\frac{1}{2}$. This means that, $$ \frac{z}{e^z- 1} = 1 + -\frac{z}{2}+\sum_{k \ge 1} \frac{B_{2k}}{(2k)!}\,z^{2k}. $$

Table 1. The first few other terms are, \begin{align*} B_0 &= 1, & B_1 &= -\frac{1}{2}, & B_2 &= \frac{1}{6}, & B_{{4}}&=-\frac{1}{30}, \\ \\ B_{{6}}&=\frac{1}{42}, & B_{{8}}&=-\frac{1}{30}, &B_{{10}}&={\frac {5}{66}}, &B_{{12}}&=-{\frac {691}{2730}}, \\ \\ B_{{14}}&=\frac{7}{6}, & B_{{16}}&=-{\frac {3617}{510}}, & B_{{18}}&={\frac {43867}{798}}, &B_{{20}}&=-{\frac {174611}{330}}, \dots \end{align*}

We give here a simple demonstration of the Euler-MacLaurin formula. Let's start with the following equivalent relations, \begin{align*} P(x) &= \int_x^{x+1} Q(x)\,dx, \\ P'(x) &= Q(x+1) - Q(x). \end{align*} Expressing the second relation using the operators we get, $ D P(x) = (e^D - 1)\,Q(x) $ or again, $$ \frac{D}{e^D - 1} \, P(x) = Q(x). $$

Using what we know about Bernoulli numbers, $$ \frac{D}{e^D-1}\,P(x) = P(x) - \frac{1}{2}\,P'(x) + \sum_{k = 1}^{\infty} \frac{B_{2k}}{(2k)!} \, D^{2k} P(x). $$ We replace each parts of this expression so as to remove $P$ from the landscape, $$ \frac{1}{2}\,\left[\,Q(x+1) + Q(x)\,\right] = \int_x^{x+1} Q(x)\,dx + \sum_{k = 1}^{\infty} \frac{B_{2k}}{(2k)!} \, \left[\,Q^{(2k-1)}(x+1) - Q^{(2k-1)}(x)\,\right]. $$ Finally, we sum this equation for $x = 0, ..., n-1$ and we get the Euler-MacLaurin formula, $$ \frac{1}{2}\left[\, Q(n) + Q(0) \,\right] + \sum_{k = 1}^{n-1} Q(k) = \int_0^n Q(x)\,dx + \sum_{k = 1}^\infty \frac{B_{2k}}{(2k)!} \,\left[\, Q^{(2k-1)}(n) - Q^{(2k-1)}(0) \,\right]. $$

There is more to the topic than just Bernoulli numbers, there are also polynomials. Bernoulli polynomials $B_n(x)$ are defined by the following exponential generating series, $$ \frac{z}{e^z-1}\,e^{zx} = \sum_{n = 0}^{\infty} B_n(x) \, \frac{z^n}{n!} $$

The polynomials are obviously related to the numbers. If you take $x=0$ in that generating series you recover the generating series of Bernoulli numbers, $$ \left. \frac{z}{e^z-1}\,e^{zx}\right|_{x=0} = \frac{z}{e^z-1}. $$ This means that we have $B_n(0) = B_n$. Bernoulli numbers are the constant terms of the Bernoulli polynomials. There is also a formula for the other way around, which express Bernoulli polynomials from Bernoulli numbers.

Theorem 1.*
For all $n \ge 0$,
$$
B_n(x) = \sum_{k = 0}^n \left(n \atop k \right) B_k \,x^{n-k}.
$$
*

Table 2. \begin{align*} B_{{0}} (x) &=1\\ B_{{1}} (x) &=x-1/2\\ B_{{2}} (x) &=1/6+{x}^{2}-x\\ B_{{3}} (x) &=1/2\,x+{x}^{3}-3/2\,{x}^{2}\\ B_{{4}} (x) &=-1/30+{x}^{2}+{x}^{4}-2\,{x}^{3}\\ B_{{5}} (x) &=-1/6\,x+5/3\,{x}^{3}+{x}^{5}-5/2\,{x}^{4}\\ B_{{6}} (x) &=1/42-1/2\,{x}^{2}+5/2\,{x}^{4}+{x}^{6}-3\,{x}^{5}\\ B_{{7}} (x) &=1/6\,x-7/6\,{x}^{3}+7/2\,{x}^{5}+{x}^{7}-7/2\,{x}^{6}\\ B_{{8}} (x) &=-1/30+2/3\,{x}^{2}-7/3\,{x}^{4}+14/3\,{x}^{6}+{x}^{8}-4\,{x}^{7}\\ B_{{9}} (x) &=-3/10\,x+2\,{x}^{3}-21/5\,{x}^{5}+6\,{x}^{7}+{x}^{9}-9/2\,{x}^{8}\\ B_{{10}} (x) &=5/66-3/2\,{x}^{2}+5\,{x}^{4}-7\,{x}^{6}+15/2\,{x}^{8}+{x}^{10}-5\,{x}^{9} \end{align*}

This last theorem is important because it reveals the *umbral* nature of the Bernoulli polynomials.
Umbral calculus is the *black magic* of nineteenth century. It has been put on clear rigorous ground by Steven Roman and Gian-Carlo Rota in the years 1970.

We shall use umbral calculus to prove the following translation theorem,
$$
B_n(x+y) = \sum_{k=0}^n \left(n \atop k\right) B_k(x)\, y^{n-k}
$$
The resemblance of this identity with the Newton binomial theorem is striking.
What if we introduce a letter, say $b$ and that we write,
$$
(b+x)^n = \sum_{k=0}^n \left(n \atop k\right) b^k\, x^{n-k}
$$
If we consider the linear operator $L$ on the space of polynomials that sends $b^k$ to $B_k$ and we apply it to that binomial formula we get,
\begin{align*}
L\left[\, (b+x)^n \,\right] &= L\left[\, \sum_{k=0}^n \left(n \atop k\right) b^k\, x^{n-k} \,\right] \\
&= \sum_{k=0}^n \left(n \atop k\right) L[\,b^k\,]\, x^{n-k} \\
&= \sum_{k=0}^n \left(n \atop k\right) B_k\, x^{n-k} = B_k(x)
\end{align*}
The following computation leads to the translation theorem,
\begin{align*}
B_n(x+y) &= L\left[\, (b+x+y)^n \,\right] \\
&= L\left[\, \sum_{k=0}^{n} \left(n \atop k\right) (b+x)^k\, y^{n-k} \,\right] \\
&= \sum_{k=0}^{n} \left(n \atop k\right) L\left[\, (b+x)^k \,\right] y^{n-k} \\
&= \sum_{k=0}^{n} \left(n \atop k\right) B_k(x)\, y^{n-k}
\end{align*}
In that computation, the $L$ comes from Rota 1970, but the subject was known well before Rota and comes from the work of Sylvester, Cayley etc... Imagine for a moment that you take the previous computation and remove the $L$. By doing that, you're doing umbral calculus in the old fashioned way. No doubt how *shadowy* the whole subject used to look. By the way, the word umbral actually *means* shadowy.
This is why I call it black magic of nineteenth century.

Hmm there is a famous formula for the sum of powers, namely the sum, $$ \sum_{k=0}^{n} k^s = \frac{B_{s+1}(n+1)-B_{s+1}(0)}{s+1}, $$ where $B_n(x)$ are the Bernoulli polynomials. It is called the Faulhaber's formula, after Johann Faulhaber who devoted half of his life to the explicit computation of the first polynomials $s = 0$ to $17$. In most textbooks it is derived as a consequence of Euler-MacLaurin formula. The two subjects are closely connected for sure but I personally don't like the usual proof. My favorite proof of this statement use operators and generating functions.

Ok the starting point is to realize that applying the derivation $D = \frac{d}{dx}$ to the exponential $e^{xz}$ is the same as multiplying it by $z$. Pretty dull yes, but look at what follows. Repeated application of $D$ to $e^{xz}$ multiplies it with $z$ each time, $$ D^n \, e^{xz} = z^n \, e^{xz} $$ Applying a series $G(D)$ of $D$ to $e^{xz}$ likewise results in applying it by $G(z)$ by extended linearity. What do we get if we use the generating series of Bernoulli numbers, $$ \frac{D}{e^D - 1} \, e^{xz} = \frac{z}{e^z - 1} \, e^{xz} $$ which means, $$ \sum_{n\ge 0} \left(\frac{D}{e^D - 1}\, x^n\right)\frac{z^n}{n!} = \sum_{n \ge 0} B_n(x)\frac{z^n}{n!} $$ and so, $$ \frac{D}{e^D - 1}\, x^n = B_n(x) $$ Or equivalently, $$ D\, x^n = (e^D - 1) \, B_n(x) $$

The operator $\Delta = e^D - 1$ is called the *forward difference operator*. It sends a function $f(x)$ to $f(x+1) - f(x)$. There is a world of relations and formulae with this operator.
From that, we see that the derivation of ordinary monomials equals the forward difference operator applied to the Bernoulli polynomials.
$$
D\, x^n = \Delta \, B_n(x)
$$
Explicitely, we have,
$$
n\, x^{n-1} = B_n(x+1) - B_n(x)
$$
Putting $n = s+1$ and summing from $x = 0$ to $n$ we finally get,
$$
\sum_{k=0}^{n} k^s = \frac{B_{s+1}(n+1)-B_{s+1}(0)}{s+1}.
$$

Table 3. \begin{align*} \sum_{k=0}^{n}k&=\frac{{n}^{2}+n}{2} \\ \sum_{k=0}^{n}{k}^{2}&=\frac{2\,n^3+3\,n^2+n}{6} \\ \sum_{k=0}^{n}{k}^{3}&=\frac{n^4 + 2\,n^3+n^2}{4} \\ \sum_{k=0}^{n}{k}^{4}&=\frac{6\,n^5+15\,n^4+10\,n^3-n}{30} \\ \sum_{k=0}^{n}{k}^{5}&=\frac{2\,n^6+6\,n^5+5\,n^4-n^2}{12} \\ \sum_{k=0}^{n}{k}^{6}&=\frac{6\,n^7+21\,n^6+21\,n^5-7\,n^3+n}{42} \\ \sum_{k=0}^{n}{k}^{7}&=\frac{3\,n^8+12\,n^7+14\,n^6-7\,n^4+2\,n^2}{24} \\ \end{align*}

Some people told me that they have trouble understanding the part on Taylor formula.

One frequently asked question was. I always saw the Taylor formula written like that, $$ f(x) = \sum_{n = 0}^\infty \frac{f^{(n)}(x_0)}{n!} (x-x_0)^n $$

Answer. Well, you are totally right. Just replace all at once $x$ by $x + h$ and $x_0$ by $x$ and you get. $$ f(x + h) = \sum_{n = 0}^\infty f^{(n)}(x) \frac{h^n}{n!} $$ This is the same formula.

Another question was. I don't get what is the exponential of $D$.

Answer. No problem, just consider, $$ e^{hD} = \sum_{n = 0}^{\infty} \frac{(hD)^n}{n!} = \sum_{n=0}^{\infty} D^n \frac{h^n}{n!} $$ Ok now apply this operator to the function $f(x)$, $$ e^{hD} f(x) = \sum_{n=0}^{\infty} D^n f(x) \frac{h^n}{n!} = \sum_{n=0}^{\infty} f^{(n)}(x) \frac{h^n}{n!} $$ And by Taylor formula this is, $$ e^{hD} f(x) = f(x + h) $$

Hope this helps ;-)