--- jupytext: formats: ipynb,md:myst,py:light text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.13.8 kernelspec: display_name: Python 3 language: python name: python3 --- ```{code-cell} :tags: [hide-cell] import mmf_setup; mmf_setup.nbinit() import os from pathlib import Path FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/' os.makedirs(FIG_DIR, exist_ok=True) import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL) %matplotlib inline import numpy as np, matplotlib.pyplot as plt try: from myst_nb import glue except: glue = None ``` (sec:EulerMaclaurin)= # 12. [Euler-Maclaurin Formula][] :::{margin} The [trapezoidal rule][] is easy to remember: it is just the sum of $f(x_i)$ at all of the equally spaced points with a factor of $1/2$ for the endpoints, with the whole thing multiplied by the lattice spacing $h$ or "$\d{x}$". ::: 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*} :::{margin} In quantum mechanics \begin{gather*} e^{\lambda \op{p}/\I\hbar}\psi(x) = \psi(x+\lambda). \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*} :::{margin} The last formula makes sense since $\op{E}$ and $\op{D}$ commute: $[\op{E}, \op{D}] = 0$. Otherwise, we would have had to be careful about the order: \begin{gather*} (\op{E}^N-1)\op{D}^{-1}f(a). \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*} :::{margin} Hint: recall that \begin{gather*} \sum_{n=0}^{N-1} r^n = \frac{r^{N} - 1}{r-1},\\ \frac{e^t}{e^t-1} = \sum_{n=0}^{\infty} \frac{B_n t^n}{n^!}. \end{gather*} ::: :::{admonition} Do It! Derive this. :class: dropdown \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. ::: *{cite}`Arfken: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 {cite}`Graham: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 {cite}`Graham: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 {ref}`sec:AsymptoticSeries`. [Bernoulli polynomials]: [Bernoulli numbers]: [Asymptotic series]: [Euler-Maclaurin formula]: [Trapezoidal rule]: