Euler–Maclaurin formula

In mathematics, the Euler–Maclaurin formula provides a powerful connection between integrals and sums. It can be used to approximate integrals by finite sums, or conversely to evaluate finite sums and infinite series using integrals and the machinery of calculus. For example, many asymptotic expansions are derived from the formula, and Faulhaber's formula for the sum of powers is an immediate consequence.

The formula was discovered independently by Leonhard Euler and Colin Maclaurin around 1735 (and later generalized as Darboux's formula). Euler needed it to compute slowly converging infinite series while Maclaurin used it to calculate integrals.

The formula[edit]

If ${\displaystyle m}$ and ${\displaystyle n}$ are natural numbers and ${\displaystyle f(x)}$ is a complex or real valued continuous function for real numbers ${\displaystyle x}$ in the interval ${\displaystyle [m,n],}$ then the integral

${\displaystyle I=\int _{m}^{n}f(x)\,dx}$

can be approximated by the sum (or vice versa)

${\displaystyle S=f(m+1)+\cdots +f(n-1)+f(n)}$

(see rectangle method). The Euler–Maclaurin formula provides expressions for the difference between the sum and the integral in terms of the higher derivatives ${\displaystyle f^{(k)}(x)}$ evaluated at the end points of the interval, that is to say when ${\displaystyle x=m}$ and ${\displaystyle x=n.}$

Explicitly, for ${\displaystyle p}$ a positive integer and a function ${\displaystyle f(x)}$ that is ${\displaystyle p}$ times continuously differentiable in the interval ${\displaystyle [m,n],}$ we have

${\displaystyle S-I=\sum _{k=1}^{p}{{\frac {B_{k}}{k!}}(f^{(k-1)}(n)-f^{(k-1)}(m))}+R_{p}}$

where ${\displaystyle B_{k}}$ is the ${\displaystyle k}$th Bernoulli number (with ${\displaystyle B_{1}=1/2}$) and ${\displaystyle R_{p}}$ is an error term which is normally small for suitable values of ${\displaystyle p}$ and depends on ${\displaystyle n,m,p,}$ and ${\displaystyle f.}$

The formula is often written with the subscript taking only even values, since the odd Bernoulli numbers are zero except for ${\displaystyle B_{1},}$ in which case we have^[1][2]

${\displaystyle \sum _{i=m}^{n}f(i)=\int _{m}^{n}f(x)\,dx+{\frac {f(n)+f(m)}{2}}+\sum _{k=1}^{\lfloor p/2\rfloor }{\frac {B_{2k}}{(2k)!}}(f^{(2k-1)}(n)-f^{(2k-1)}(m))+R_{p}}$

or alternatively

${\displaystyle \sum _{i=m+1}^{n}f(i)=\int _{m}^{n}f(x)\,dx+{\frac {f(n)-f(m)}{2}}+\sum _{k=1}^{\lfloor p/2\rfloor }{\frac {B_{2k}}{(2k)!}}(f^{(2k-1)}(n)-f^{(2k-1)}(m))+R_{p}.}$

The first few Bernoulli numbers of even index are:

${\displaystyle B_{2}={\frac {1}{6}},\quad B_{4}=-{\frac {1}{30}},\quad B_{6}={\frac {1}{42}},\quad B_{8}=-{\frac {1}{30}}.}$

We may write the Euler-Maclaurin formula then as

${\displaystyle \sum _{i=m}^{n}f(i)-\int _{m}^{n}f(x)~{\rm {d}}x={\frac {f(m)+f(n)}{2}}+{\frac {1}{6}}{\frac {f'(n)-f'(m)}{2!}}-{\frac {1}{30}}{\frac {f'''(n)-f'''(m)}{4!}}+\cdots +B_{2k}{\frac {f^{(2k-1)}(n)-f^{(2k-1)}(m)}{(2k)!}}+R_{2k}.}$

The Bernoulli polynomials and periodic function[edit]

The formula is derived below using repeated integration by parts applied to successive intervals ${\displaystyle [r,r+1]}$ for integers ${\displaystyle r=m,m+1,\cdots ,n-1.}$ The derivation uses the periodic Bernoulli functions, ${\displaystyle P_{k}(x),}$ which are defined in terms of the Bernoulli polynomials ${\displaystyle B_{k}(x)}$ for ${\displaystyle k=0,1,2\cdots .}$

The Bernoulli polynomials may be defined recursively by

${\displaystyle {\begin{aligned}B_{0}(x)&=1\\B_{k}'(x)&=kB_{k-1}(x){\text{ and }}\int _{0}^{1}B_{k}(x)\,dx=0{\text{ for }}k\geq 1\end{aligned}}}$

and the periodic Bernoulli functions are defined as

${\displaystyle P_{k}(x)=B_{k}\left(x-\lfloor x\rfloor \right)}$

where ${\displaystyle \lfloor x\rfloor }$ denotes the largest integer that is not greater than ${\displaystyle x}$ so that ${\displaystyle x-\lfloor x\rfloor }$ always lies in the interval ${\displaystyle [0,1).}$

It can be shown that ${\displaystyle B_{k}(1)=B_{k}(0)}$ for all ${\displaystyle k\neq 1}$ so that except for ${\displaystyle P_{1}(x),}$ all the periodic Bernoulli functions are continuous. The functions ${\displaystyle P_{k}(x)}$ are sometimes written as ${\displaystyle {\tilde {B}}_{k}(x).}$

The remainder term[edit]

The remainder term ${\displaystyle R_{p}}$ can be written as

${\displaystyle R_{p}=(-1)^{p+1}\int _{m}^{n}f^{(p)}(x){P_{p}(x) \over p!}\,dx.}$

When ${\displaystyle k>0,}$ it can be shown that

${\displaystyle \left|B_{k}(x)\right|\leq {\frac {2\cdot k!}{(2\pi )^{k}}}\zeta (k)}$

where ${\displaystyle \zeta }$ denotes the Riemann zeta function; one approach to prove this inequality is to obtain the Fourier series for the polynomials ${\displaystyle B_{k}(x).}$ The bound is achieved for even ${\displaystyle k}$ when ${\displaystyle x}$ is zero. The term ${\displaystyle \zeta (k)}$ may be omitted for odd ${\displaystyle k}$ but the proof in this case is more complex (see Lehmer).^[3] Using this inequality, the size of the remainder term can be estimated using

${\displaystyle \left|R_{p}\right|\leq {\frac {2\zeta (p)}{(2\pi )^{p}}}\int _{m}^{n}\left|f^{(p)}(x)\right|\ \,dx.}$

Applicable formula[edit]

We can use the formula as a means of approximating a finite integral, with the following simple formula:^[4]

${\displaystyle {\begin{aligned}I&=\int _{x_{0}}^{x_{N}}f(x)\,\mathrm {d} x\\[6pt]&=h\left({\frac {f_{0}}{2}}+f_{1}+f_{2}+\cdots +f_{N-1}+{\frac {f_{N}}{2}}\right)+{\frac {h^{2}}{12}}[f'_{0}-f'_{N}]-{\frac {h^{4}}{720}}[f'''_{0}-f'''_{N}]+\cdots ,\end{aligned}}}$

where ${\displaystyle N}$ is the number of points in the interval of integration from ${\displaystyle x_{0}}$ to ${\displaystyle x_{N}}$ and ${\displaystyle h}$ is the distance between points so that ${\displaystyle h=(x_{N}-x_{0})/N.}$ Note the series above is usually not convergent; indeed, often the terms will increase rapidly after a number of iterations. Thus, attention generally needs to be paid to the remainder term.

This may be viewed as an extension of the trapezoid rule by the inclusion of correction terms.^[4]

Applications[edit]

The Basel problem[edit]

The Basel problem asks to determine the sum

${\displaystyle 1+{\frac {1}{4}}+{\frac {1}{9}}+{\frac {1}{16}}+{\frac {1}{25}}+\cdots =\sum _{n=1}^{\infty }{\frac {1}{n^{2}}}.}$

Euler computed this sum to 20 decimal places with only a few terms of the Euler–Maclaurin formula in 1735. This probably convinced him that the sum equals ${\displaystyle \pi ^{2}/6,}$ which he proved in the same year.^[5] Parseval's identity for the Fourier series of ${\displaystyle f(x)=x}$ gives the same result.

Sums involving a polynomial[edit]

If ${\displaystyle f}$ is a polynomial and ${\displaystyle p}$ is big enough, then the remainder term vanishes. For instance, if ${\displaystyle f(x)=x^{3},}$ we can choose ${\displaystyle p=2}$ to obtain after simplification

${\displaystyle \sum _{i=0}^{n}i^{3}=\left({\frac {n(n+1)}{2}}\right)^{2}}$

(see Faulhaber's formula).

Numerical integration[edit]

The Euler–Maclaurin formula is also used for detailed error analysis in numerical quadrature. It explains the superior performance of the trapezoidal rule on smooth periodic functions and is used in certain extrapolation methods. Clenshaw–Curtis quadrature is essentially a change of variables to cast an arbitrary integral in terms of integrals of periodic functions where the Euler–Maclaurin approach is very accurate (in that particular case the Euler–Maclaurin formula takes the form of a discrete cosine transform). This technique is known as a periodizing transformation.

Asymptotic expansion of sums[edit]

In the context of computing asymptotic expansions of sums and series, usually the most useful form of the Euler–Maclaurin formula is

${\displaystyle \sum _{n=a}^{b}f(n)\sim \int _{a}^{b}f(x)\,dx+{\frac {f(b)+f(a)}{2}}+\sum _{k=1}^{\infty }\,{\frac {B_{2k}}{(2k)!}}(f^{(2k-1)}(b)-f^{(2k-1)}(a))}$

where ${\displaystyle a}$ and ${\displaystyle b}$ are integers.^[6] Often the expansion remains valid even after taking the limits ${\displaystyle a\rightarrow -\infty }$ or ${\displaystyle b\rightarrow +\infty }$ or both. In many cases the integral on the right-hand side can be evaluated in closed form in terms of elementary functions even though the sum on the left-hand side cannot. Then all the terms in the asymptotic series can be expressed in terms of elementary functions. For example,

${\displaystyle \sum _{k=0}^{\infty }{\frac {1}{(z+k)^{2}}}\sim \underbrace {\int _{0}^{\infty }{\frac {1}{(z+k)^{2}}}\,dk} _{=1/z}+{\frac {1}{2z^{2}}}+\sum _{t=1}^{\infty }{\frac {B_{2t}}{z^{2t+1}}}.}$

Here the left-hand side is equal to ${\displaystyle \psi ^{(1)}(z),}$ namely the first-order polygamma function defined through ${\displaystyle \psi ^{(1)}(z)=D^{2}(\log \Gamma (z));}$ the gamma function ${\displaystyle \Gamma (z)}$ is equal to ${\displaystyle (z-1)!}$ if ${\displaystyle z}$ is a positive integer. This results in an asymptotic expansion for ${\displaystyle \psi ^{(1)}(z).}$ That expansion, in turn, serves as the starting point for one of the derivations of precise error estimates for Stirling's approximation of the factorial function.

Examples[edit]

If s is an integer greater than 1 we have:

${\displaystyle \sum _{k=1}^{n}{\frac {1}{k^{s}}}\approx {\frac {1}{s-1}}+{\frac {1}{2}}-{\frac {1}{(s-1)n^{s-1}}}+{\frac {1}{2n^{s}}}+\sum _{i=1}{\frac {B_{2i}}{(2i)!}}\left[{\frac {(s+2i-2)!}{(s-1)!}}-{\frac {(s+2i-2)!}{(s-1)!n^{s+2i-1}}}\right]}$

Collecting the constants into a value of the Riemann zeta function, we can write an asymptotic expansion:

${\displaystyle \sum _{k=1}^{n}{\frac {1}{k^{s}}}\sim \zeta (s)-{\frac {1}{(s-1)n^{s-1}}}+{\frac {1}{2n^{s}}}-\sum _{i=1}{\frac {B_{2i}}{(2i)!}}{\frac {(s+2i-2)!}{(s-1)!n^{s+2i-1}}}}$

For s equal to 2 this simplifies to

${\displaystyle \sum _{k=1}^{n}{\frac {1}{k^{2}}}\sim \zeta (2)-{\frac {1}{n}}+{\frac {1}{2n^{2}}}-\sum _{i=1}{\frac {B_{2i}}{n^{2i+1}}}}$

or

${\displaystyle \sum _{k=1}^{n}{\frac {1}{k^{2}}}\sim {\frac {\pi ^{2}}{6}}-{\frac {1}{n}}+{\frac {1}{2n^{2}}}-{\frac {1}{6n^{3}}}+{\frac {1}{30n^{5}}}-{\frac {1}{42n^{7}}}.}$

We can also derive (from Equation 1 below) the perhaps not-so-useful formula:

${\displaystyle \sum _{k=1}^{n}{\frac {1}{k^{s}}}={\frac {1}{n^{s-1}}}+s\int _{1}^{n}{\frac {\lfloor x\rfloor }{x^{s+1}}}dx\qquad {\text{with }}\quad s\in \mathbb {R} \setminus \{1\}}$

When s is 1, we get the following for the so-called harmonic numbers:

${\displaystyle {\begin{aligned}\sum _{k=1}^{n}{\frac {1}{k}}&=\log n+{\frac {1}{2}}+{\frac {1}{2n}}-\int _{1}^{n}{\frac {x-\lfloor x\rfloor -1/2}{x^{2}}}dx\\&\approx \log n+\gamma +{\frac {1}{2n}}-{\frac {1}{12n^{2}}}+{\frac {1}{120n^{4}}}-{\frac {1}{252n^{6}}}\end{aligned}}}$

where ${\displaystyle \gamma \approx 0.5772156649015328606065}$ is the Euler constant,

Proofs[edit]

Derivation by mathematical induction[edit]

We follow the argument given in Apostol.^[1]

The Bernoulli polynomials B_n(x) and the periodic Bernoulli functions P_n(x) for n = 0, 1, 2, ... were introduced above.

The first several Bernoulli polynomials are

${\displaystyle {\begin{aligned}B_{1}(x)&=x-{\frac {1}{2}}\\B_{2}(x)&=x^{2}-x+{\frac {1}{6}}\\B_{3}(x)&=x^{3}-{\frac {3}{2}}x^{2}+{\frac {1}{2}}x\\B_{4}(x)&=x^{4}-2x^{3}+x^{2}-{\frac {1}{30}}\\&\,\,\,\vdots \end{aligned}}}$

The values B_n(0) are the Bernoulli numbers. Notice that for n ≠ 1 we have

${\displaystyle B_{n}(0)=B_{n}(1)=B_{n}\quad (n{\text{th Bernoulli number}})}$

For n = 1,

${\displaystyle B_{1}(0)=-B_{1}(1)=B_{1}}$

The functions P_n agree with the Bernoulli polynomials on the interval [0, 1] and are periodic with period 1. Furthermore, except when n = 1, they are also continuous. Thus,

${\displaystyle P_{n}(0)=P_{n}(1)=B_{n}\quad (n\neq 1)}$

Let k be an integer, and consider the integral

${\displaystyle \int _{k}^{k+1}f(x)\,dx=\int _{k}^{k+1}u\,dv}$

where

${\displaystyle {\begin{aligned}u&=f(x)\\du&=f'(x)\,dx\\dv&=P_{0}(x)\,dx&&{\text{since }}P_{0}(x)=1\\v&=P_{1}(x)\end{aligned}}}$

Integrating by parts, we get

${\displaystyle {\begin{aligned}\int _{k}^{k+1}f(x)\,dx&={\Big [}uv{\Big ]}_{k}^{k+1}-\int _{k}^{k+1}v\,du\\[8pt]&={\Big [}f(x)P_{1}(x){\Big ]}_{k}^{k+1}-\int _{k}^{k+1}f'(x)P_{1}(x)\,dx\\[8pt]&=B_{1}(1)f(k+1)-B_{1}(0)f(k)-\int _{k}^{k+1}f'(x)P_{1}(x)\,dx\end{aligned}}}$

Using ${\displaystyle B_{1}(0)=-1/2}$, ${\displaystyle B_{1}(1)=1/2}$, and summing the above from k = 0 to k = n − 1, we get

${\displaystyle {\begin{aligned}\int _{0}^{n}f(x)\,dx&=\int _{0}^{1}f(x)\,dx+\dotsb +\int _{n-1}^{n}f(x)\,dx\\[8pt]&={\frac {f(0)}{2}}+f(1)+\dotsb +f(n-1)+{f(n) \over 2}-\int _{0}^{n}f'(x)P_{1}(x)\,dx\end{aligned}}}$

Adding (f(n) − f(0))/2 to both sides and rearranging, we have

${\displaystyle \sum _{k=1}^{n}f(k)=\int _{0}^{n}f(x)\,dx+{f(n)-f(0) \over 2}+\int _{0}^{n}f'(x)P_{1}(x)\,dx\qquad (1)}$

The last two terms therefore give the error when the integral is taken to approximate the sum.

Next, consider

${\displaystyle \int _{k}^{k+1}f'(x)P_{1}(x)\,dx=\int _{k}^{k+1}u\,dv}$

where

${\displaystyle {\begin{aligned}u&=f'(x)\\du&=f''(x)\,dx\\dv&=P_{1}(x)\,dx\\v&={\frac {1}{2}}P_{2}(x)\end{aligned}}}$

Integrating by parts again, we get

${\displaystyle {\begin{aligned}{\Big [}uv{\Big ]}_{k}^{k+1}-\int _{k}^{k+1}v\,du&=\left[{f'(x)P_{2}(x) \over 2}\right]_{k}^{k+1}-{1 \over 2}\int _{k}^{k+1}f''(x)P_{2}(x)\,dx\\&={B_{2} \over 2}(f'(k+1)-f'(k))-{1 \over 2}\int _{k}^{k+1}f''(x)P_{2}(x)\,dx\end{aligned}}}$

Then summing from k = 0 to k = n − 1, and then replacing the last integral in (1) with what we have thus shown to be equal to it, we have

${\displaystyle \sum _{k=1}^{n}f(k)=\int _{0}^{n}f(x)\,dx+{f(n)-f(0) \over 2}+{\frac {B_{2}}{2}}(f'(n)-f'(0))-{1 \over 2}\int _{0}^{n}f''(x)P_{2}(x)\,dx.}$

By now the reader will have guessed that this process can be iterated. In this way we get a proof of the Euler–Maclaurin summation formula which can be formalized by mathematical induction, in which the induction step relies on integration by parts and on the identities for periodic Bernoulli functions.

See also[edit]

• Cesàro summation
• Euler summation
• Gauss–Kronrod quadrature formula
• Darboux's formula
• Euler–Boole summation

Notes[edit]

References[edit]