Euler–Maclaurin formula

Applications

The Basel problem

The Basel problem is to determine the sum $${\displaystyle 1+\frac{1}{4}+\frac{1}{9}+\frac{1}{16}+\frac{1}{25}+\cdots =\sum _{n=1}^{\mathrm{\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 ⁠π²/6⁠, which he proved in the same year.

Sums involving a polynomial

If f is a polynomial and p is big enough, then the remainder term vanishes. For instance, if f(x) = x³, we can choose p = 2 to obtain, after simplification, $${\displaystyle \sum _{i=0}^n{i}^3={\left(\frac{n(n+1)}{2}\right)}^2.}$$

Approximation of integrals

The formula provides a means of approximating a finite integral. Let a < b be the endpoints of the interval of integration. Fix N, the number of points to use in the approximation, and denote the corresponding step size by h = ⁠b − a/N − 1⁠. Set x_i = a + (i − 1)h, so that x₁ = a and x_N = b. Then: $${\displaystyle \begin{array}{cc}I& ={\int }_a^{b}f(x)dx\\ & \sim h\left(\frac{f(x_1)}{2}+f(x_2)+\cdots +f(x_{N-1})+\frac{f(x_N)}{2}\right)+\frac{h^2}{12}[f^{\prime }(x_1)-f^{\prime }(x_N)]-\frac{h^4}{720}[f^{‴}(x_1)-f^{‴}(x_N)]+\cdots \end{array}}$$

This may be viewed as an extension of the trapezoid rule by the inclusion of correction terms. Note that this asymptotic expansion is usually not convergent; there is some p, depending upon f and h, such that the terms past order p increase rapidly. Thus, the remainder term generally demands close attention.

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

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}^{\mathrm{\infty }}\frac{B_{2k}}{(2k)!}\left(f^{(2k-1)}(b)-f^{(2k-1)}(a)\right),}$$

where a and b are integers. Often the expansion remains valid even after taking the limits a → −∞ or b → +∞ 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}^{\mathrm{\infty }}\frac{1}{(z+k{)}^2}\sim \underset{={\displaystyle \frac{1}{z}}}{\underbrace{{\int }_0^{\mathrm{\infty }}\frac{1}{(z+k{)}^2}dk}}+\frac{1}{2z^2}+\sum _{t=1}^{\mathrm{\infty }}\frac{B_{2t}}{z^{2t+1}}.}$$

Here the left-hand side is equal to ψ⁽¹⁾(z), namely the first-order polygamma function defined by

${\displaystyle {\psi }^{(1)}(z)=\frac{d^2}{d{z}^2}\ln \mathrm{\Gamma }(z);}$

the gamma function Γ(z) is equal to (z − 1)! when z is a positive integer. This results in an asymptotic expansion for ψ⁽¹⁾(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

If s is an integer greater than 1 then we have the following asymptotic expansions as ${\displaystyle n\to \mathrm{\infty }}$, using the Riemann zeta function ${\displaystyle \zeta (s)}$:

$${\displaystyle \begin{array}{cc}\sum _{k=1}^n\frac{1}{k^s}& =\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k^s}-\sum _{k=n+1}^{\mathrm{\infty }}\frac{1}{k^s}\\ & =\zeta (s)-\sum _{k=n+1}^{\mathrm{\infty }}\frac{1}{k^s}\\ & \approx \zeta (s)-{\int }_n^{\mathrm{\infty }}\frac{dk}{k^s}+\frac{1}{2n^s}-\sum _{i=1}\frac{B_{2i}}{(2i)!}\frac{(s+2i-2)!}{(s-1)!n^{s+2i-1}}\\ & \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}}.\end{array}}$$

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}+\cdots .}$$

This is related to the trigamma function:

$${\displaystyle \sum _{k=1}^n\frac{1}{k^2}=\frac{{\pi }^2}{6}-{\psi }^{(1)}(n+1)}$$

When s = 1, the sum and the integral go to infinity but the difference between them goes to a limit, the Euler–Mascheroni constant, γ ≈ 0.5772... This gives:

$${\displaystyle \begin{array}{cc}{\int }_1^n\left(\frac{1}{\lfloor x\rfloor }-\frac{1}{x}\right)dx& ={\int }_1^{\mathrm{\infty }}\left(\frac{1}{\lfloor x\rfloor }-\frac{1}{x}\right)dx-{\int }_n^{\mathrm{\infty }}\left(\frac{1}{\lfloor x\rfloor }-\frac{1}{x}\right)dx\\ & =\gamma -{\int }_n^{\mathrm{\infty }}\left(\frac{1}{\lfloor x\rfloor }-\frac{1}{x}\right)dx\\ & \approx \gamma -\frac{1}{2n}-\sum _{k=1}\frac{B_{2k}}{2k{n}^{2k}}\end{array}}$$

from which we have the asymptotic expansion:

$${\displaystyle \sum _{k=1}^n\frac{1}{k}\sim \ln n+\gamma +\frac{1}{2n}-\sum _{k=1}\frac{B_{2k}}{2k{n}^{2k}},}$$

These harmonic numbers are related to the digamma function:

$${\displaystyle \sum _{k=1}^n\frac{1}{k}=\gamma +\psi (n+1)}$$

More generally, $${\displaystyle \sum _{k=1}^n\frac{1}{k^s}=\frac{1}{n^{s-1}}+s{\int }_1^n\frac{x-\lfloor x\rfloor }{x^{s+1}}dx\text{with }s>1}$$

Proofs

Derivation by mathematical induction

We outline the argument given in Apostol.

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{array}{cc}{B}_0(x)& =1,\\ 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},\\ & ⋮\end{array}}$$

The values B_n(1) are the Bernoulli numbers B_n. Notice that for n ≠ 1 we have $${\displaystyle B_n=B_n(1)=B_n(0),}$$ and for n = 1, $${\displaystyle B_1=B_1(1)=-B_1(0).}$$

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\text{for }n\ne 1.}$$

Let k be an integer, and consider the integral $${\displaystyle {\int }_k^{k+1}f(x)dx={\int }_k^{k+1}udv,}$$ where $${\displaystyle \begin{array}{cc}u& =f(x),\\ du& =f^{\prime }(x)dx,\\ dv& =P_0(x)dx& \text{since }{P}_0(x)& =1,\\ v& =P_1(x).\end{array}}$$

Integrating by parts, we get $${\displaystyle \begin{array}{cc}{\int }_k^{k+1}f(x)dx& =[uv{]}_k^{k+1}-{\int }_k^{k+1}vdu\\ & =[f(x)P_1(x){]}_k^{k+1}-{\int }_k^{k+1}{f}^{\prime }(x)P_1(x)dx\\ & =B_1(1)f(k+1)-B_1(0)f(k)-{\int }_k^{k+1}{f}^{\prime }(x)P_1(x)dx.\end{array}}$$

Using B₁(0) = −⁠1/2⁠, B₁(1) = ⁠1/2⁠, and summing the above from k = 0 to k = n − 1, we get $${\displaystyle \begin{array}{cc}{\int }_0^{n}f(x)dx& ={\int }_0^{1}f(x)dx+\cdots +{\int }_{n-1}^{n}f(x)dx\\ & =\frac{f(0)}{2}+f(1)+\cdots +f(n-1)+\frac{f(n)}{2}-{\int }_0^n{f}^{\prime }(x)P_1(x)dx.\end{array}}$$

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+\frac{f(n)-f(0)}{2}+{\int }_0^n{f}^{\prime }(x)P_1(x)dx.}$$

This is the p = 1 case of the summation formula. To continue the induction, we apply integration by parts to the error term: $${\displaystyle {\int }_k^{k+1}{f}^{\prime }(x)P_1(x)dx={\int }_k^{k+1}udv,}$$ where $${\displaystyle \begin{array}{cc}u& =f^{\prime }(x),\\ du& =f^{″}(x)dx,\\ dv& =P_1(x)dx,\\ v& =\frac{1}{2}{P}_2(x).\end{array}}$$

The result of integrating by parts is $${\displaystyle \begin{array}{cc}[uv{]}_k^{k+1}-{\int }_k^{k+1}vdu& ={\left[\frac{f^{\prime }(x)P_2(x)}{2}\right]}_k^{k+1}-\frac{1}{2}{\int }_k^{k+1}{f}^{″}(x)P_2(x)dx\\ & =\frac{B_2}{2}(f^{\prime }(k+1)-f^{\prime }(k))-\frac{1}{2}{\int }_k^{k+1}{f}^{″}(x)P_2(x)dx.\end{array}}$$

Summing from k = 0 to k = n − 1 and substituting this for the lower order error term results in the p = 2 case of the formula, $${\displaystyle \sum _{k=1}^{n}f(k)={\int }_0^{n}f(x)dx+\frac{f(n)-f(0)}{2}+\frac{B_2}{2}(f^{\prime }(n)-f^{\prime }(0))-\frac{1}{2}{\int }_0^n{f}^{″}(x)P_2(x)dx.}$$

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 identities for periodic Bernoulli functions.