The Euler-Maclaurin Formula is the exact error \(\epsilon\) when computing a definite
integral with the trapezoidal rule:
\[\begin{gather*}
\int_{a}^{b}f(x)\d{x} = h\left(
\frac{1}{2}f(a) + f(a+h) +\dots + f(b-h) + \frac{1}{2}f(b)
\right) + \epsilon, \\
h = \frac{b-a}{N}.
\end{gather*}\]
To derive the formula, we use property that the derivative \(\op{D} = \partial_x\) generates
translations
\[\begin{gather*}
\op{D}f(x) = f'(x), \qquad
e^{h\op{D}} f(x) = f(x + h).
\end{gather*}\]
This just follows from Taylor expanding the exponential. For the purposes of this
discussion, consider
\[\begin{gather*}
\op{E} = e^{h\op{D}} = \op{1} + h\op{D} + \frac{h^2}{2!}\op{D}^2 + \cdots,
\qquad \op{E}^{n}f(a) = f(a+hn).
\end{gather*}\]
If we let \(\op{D}F(x) = f(x)\) we can express the exact answer
\[\begin{gather*}
\int_{a}^{b}f(x)\d{x} = F(b)-F(a) = F(a+hN) - F(a) = (\op{E}^N - 1)F(a)\\
= \frac{\op{E}^N - 1}{\op{D}}f(a),
\end{gather*}\]
The error can thus be expressed:
\[\begin{gather*}
\epsilon = -\left.\sum_{n=2}^{\infty}\frac{B_n h^{n}f^{(n-1)}(x)}{n!}\right|^{b}_{a}
= h^{2}\frac{f'(a)-f'(b)}{12}
- \left.\sum_{n=4}^{\infty}\frac{B_n h^{n}f^{(n-1)}(x)}{n!}\right|^{b}_{a}\\
\end{gather*}\]
Do It! Derive this.
\[\begin{gather*}
\frac{1}{2}f(a) + f(a+h) +\dots + f(b-h) + \frac{1}{2}f(b) = \\
=
\sum_{n=0}^{N-1} f(a+hn) + \frac{f(b)-f(a)}{2}
=
\left(\sum_{n=0}^{N-1} \op{E}^{n} + \frac{\op{E}^N-1}{2}\right)f(a)\\
=
\left(\frac{\op{E}^{N}-1}{\op{E}-1} + \frac{\op{E}^N-1}{2}\right)f(a)
=
(\op{E}^{N}-1)
\left(\frac{1}{\op{E}-1} + \frac{1}{2}\right)f(a).
\end{gather*}\]
The error is thus
\[\begin{gather*}
\epsilon = \frac{\op{E}^{N}-1}{\op{D}}\left(1 - h\op{D}
\left(\frac{1}{\op{E}-1} + \frac{1}{2}\right)\right)f(a)\\
= \frac{\op{E}^{N}-1}{\op{D}}\left(1 -
\left(\frac{h\op{D}}{\op{E}-1} + \frac{ h\op{D}}{2}\right)\right)f(a)\\
= \frac{\op{E}^{N}-1}{\op{D}}\left(
1 - \frac{ h\op{D}}{2} - \frac{h\op{D}}{e^{h\op{D}}-1}
\right)f(a).
\end{gather*}\]
Before proceeding, think about this equation. We can Taylor-expand the term in
parentheses: Note that the leading orders cancels:
\[\begin{gather*}
\epsilon = \frac{\op{E}^{N}-1}{\op{D}}\left(
1 - \frac{h\op{D}}{2} - \frac{h\op{D}}{h\op{D} +
\frac{1}{2}h^2\op{D}^2 + O(h^3\op{D}^3)}
\right)f(a)\\
= \frac{\op{E}^{N}-1}{\op{D}}\left(
1 - \frac{h\op{D}}{2} - \frac{1}{1 +
\frac{1}{2}h\op{D} + O(h^2\op{D}^2)}
\right)f(a)\\
= \frac{\op{E}^{N}-1}{\op{D}}\left(
1 - \frac{h\op{D}}{2} - \Bigl(
1 - \frac{1}{2}h\op{D} + O(h^2\op{D}^2)\Bigr)
\right)f(a)\\
= \frac{\op{E}^{N}-1}{\op{D}}O(h^2\op{D}^2)f(a)
= (\op{E}^{N}-1)O(h^2\op{D})f(a).
\end{gather*}\]
Note also that the first factor \(\op{E}^{N} - 1\) simply means apply the formula at \(b\)
then subtract at \(a\):
\[\begin{gather*}
\epsilon = O(h^2\op{D})f(a)\Bigr|^{b}_{a}
\end{gather*}\]
Working a little harder one finds the leading order term
\[\begin{gather*}
\epsilon = \frac{-h^2 f'(x)}{12}\Bigr|^{b}_{a} + O(h^3)
= h^2\frac{f'(a)-f'(b)}{12}+ O(h^3).
\end{gather*}\]
However, using the Bernoulli numbers \(B_0=1\), \(B_1=-\tfrac{1}{2}\), \(B_2=\tfrac{1}{6}\),
\(B_3=0\), \(B_4=\tfrac{-1}{30}\), etc. we can express the entire series:
\[\begin{gather*}
\epsilon = \frac{\op{E}^{N}-1}{\op{D}}\left(
1 - \frac{ h\op{D}}{2} - \sum_{n=0}^{\infty}\frac{B_n (h\op{D})^n}{n!}
\right)f(a)\\
= \frac{\op{E}^{N}-1}{\op{D}}\left(
1 - \frac{h\op{D}}{2} - B_0 - B_1 h\op{D} -
\sum_{n=2}^{\infty}\frac{B_n (h\op{D})^n}{n!}
\right)f(a)\\
= -\left.\sum_{n=2}^{\infty}\frac{B_n h^{n}f^{(n-1)}(x)}{n!}\right|^{b}_{a}\\
= h^{2}\frac{f'(a)-f'(b)}{12}
- \left.\sum_{n=4}^{\infty}\frac{B_n h^{n}f^{(n-1)}(x)}{n!}\right|^{b}_{a}\\
\end{gather*}\]
Thus, because \(B_3=0\) we see that the next term is even higher order.
[Arfken et al., 2013] derives the same result using integration by parts and the property
\(B_n'(x) = nB'_{n-1}(x)\) of the Bernoulli polynomials.
This can be rearranged as (see [Graham et al., 1994] Eq. (9.67))
\[\begin{gather*}
\sum_{a\leq n < b} f(n) = \int_{a}^{b}f(x)\d{x}
+ \left.\sum_{n=1}^{m}\frac{B_{n}}{n!}f^{(n-1)}(x)\right|^{b}_{a} + R_{m}\\
R_m = (-1)^{m+1}\int_a^{b}\frac{B_m(\{x\})}{m!}f^{(m)}(x)\d{x}
\end{gather*}\]
where \(a\), \(b\), and \(m\) are integers, \(a\leq b\), and \(1\leq m\). Here \(\{x\}\) means the
fractional part of \(x\). I.e. \(\{\pi\} = .14159\cdots\). This gives an explicit formula
for the remainder and an explicit way to correct the approximation.
The remainder term (courtesy of S.D. Poisson – see [Graham et al., 1994]) is important,
because the sum \(\sum_{n\geq 1} (B_n/n!)f^{(n-1)}(x)|^{b}_{a}\) can diverge (it is often
an asymptotic series). Nevertheless, despite diverging, asymptotic series sometimes
provide spectacular approximations, and often, using a few terms can greatly improve
one’s approximation of the series as we shall now discuss in 12. Asymptotic Series.
