View source for Euler–Maclaurin formula

{{Use American English|date = March 2019}}
{{Short description|Summation formula}}
In [[mathematics]], the '''Euler–Maclaurin formula''' is a formula for the difference between an [[integral]] and a closely related sum.  It can be used to approximate integrals by finite sums, or conversely to evaluate finite sums and [[series (mathematics)|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.  Euler needed it to compute slowly converging infinite series while Maclaurin used it to calculate integrals.  It was later generalized to [[Darboux's formula]].

==The formula==
If <math>m</math> and <math>n</math> are [[natural number]]s and <math>f(x)</math> is a complex or real valued [[continuous function]] for [[real number]]s <math>x</math> in the interval <math>[m,n]</math>, then the integral

:<math>I = \int_m^n f(x)\,dx</math>

can be approximated by the sum (or vice versa)

:<math>S = f(m + 1) + \cdots + f(n - 1) + f(n)</math>

(see [[Riemann sum#Right Riemann sum|rectangle method]]). The Euler–Maclaurin formula provides expressions for the difference between the sum and the integral in terms of the higher derivatives <math>f^{(k)}(x)</math> evaluated at the end points of the interval, that is to say when <math>x=m</math> and <math>x=n</math>.

Explicitly, for <math>p</math> a positive integer and a function <math>f(x)</math> that is <math>p</math> times continuously differentiable in the interval <math>[m, n]</math>, we have

:<math> S - I = \sum_{k=1}^p {\frac{B_k}{k!} (f^{(k - 1)}(n) - f^{(k - 1)}(m))} + R_p,</math>

where <math>B_k</math> is the <math>k</math>th [[Bernoulli numbers|Bernoulli number]] (with <math>B_1=1/2</math>) and <math>R_p</math> is an [[numerical integration|error term]] which depends on <math>n</math>, <math>m</math>, <math>p</math>, and <math>f</math> and is usually small for suitable values of <math>p</math>.

The formula is often written with the subscript taking only even values, since the odd Bernoulli numbers are zero except for <math>B_1</math>.  In this case we have<ref name=":0" /><ref name="DLMF">{{cite web|url=http://dlmf.nist.gov/2.10|title=Digital Library of Mathematical Functions: Sums and Sequences|publisher=[[National Institute of Standards and Technology]]}}</ref>

:<math>\sum_{i=m}^n f(i) =
    \int^n_m 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,
</math>

or alternatively

:<math>\sum_{i=m+1}^n f(i) =
    \int^n_m 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.
</math>

===The remainder term===
{{seealso|Bernoulli polynomials}}

The remainder term arises because the integral is usually not exactly equal to the sum.  The formula may be derived by applying repeated [[integration by parts]] to successive intervals <math>[r, r+1]</math> for <math>r = m, m + 1, \ldots, n - 1</math>.  The boundary terms in these integrations lead to the main terms of the formula, and the leftover integrals form the remainder term.

The remainder term has an exact expression in terms of the periodized Bernoulli functions <math>P_k(x)</math>.  The Bernoulli polynomials may be defined recursively by <math>B_0(x) = 1</math> and, for <math>k \ge 1</math>,
:<math>\begin{align}
  B_k'(x) &= kB_{k - 1}(x), \\
  \int_0^1 B_k(x)\,dx &= 0.
\end{align}</math>
The periodized Bernoulli functions are defined as
:<math>P_k(x) = B_k(x - \lfloor x\rfloor),</math>
where <math>\lfloor x \rfloor</math> denotes the largest integer less than or equal to <math>x</math> (so that <math>x - \lfloor x \rfloor</math> always lies in the interval <math>[0,1)</math>).

With this notation, the remainder term <math>R_p</math> equals
:<math>R_p = (-1)^{p+1}\int_m^n f^{(p)}(x) \frac{P_p(x)}{p!}\,dx. </math>

When <math>k>0</math>, it can be shown that

:<math>|B_k(x)| \le \frac{2 \cdot k!}{(2\pi)^k}\zeta(k),</math>

where <math>\zeta</math>  denotes the [[Riemann zeta function]]; one approach to prove this inequality is to obtain the Fourier series for the polynomials <math>B_k(x)</math>.   The bound is achieved for even <math>k</math> when <math>x</math> is zero.  The term <math>\zeta(k)</math> may be omitted for odd <math>k</math> but the proof in this case is more complex (see Lehmer).<ref name="Lehmer">{{cite journal|last1=Lehmer|first1=D. H. |title=On the maxima and minima of Bernoulli polynomials | date=1940 | journal=The American Mathematical Monthly|volume=47|issue=8|pages=533–538 |doi=10.2307/2303833}}</ref> Using this inequality, the size of the remainder term can be estimated as

:<math>|R_p| \leq \frac{2 \zeta(p)}{(2\pi)^p}\int_m^n |f^{(p)}(x)|\,dx.</math>

===Low-order cases===

The Bernoulli numbers from <math>B_1</math> to <math>B_7</math> are <math>\tfrac{1}{2}, \tfrac{1}{6}, 0, -\tfrac{1}{30}, 0, \tfrac{1}{42}, 0.</math> Therefore the low-order cases of the Euler-Maclaurin formula are:

:<math>\begin{align}
\sum_{i=m}^n f(i) - \int_m^n f(x)\,dx &= \frac{f(m)+f(n)}{2} + \int_m^n f'(x)P_1(x)\,dx \\
&=\frac{f(m)+f(n)}{2} + \frac{1}{6}\frac{f'(n) - f'(m)}{2!} - \int_m^n f''(x)P_2(x)\,dx \\
&=\frac{f(m)+f(n)}{2} + \frac{1}{6}\frac{f'(n) - f'(m)}{2!} + \int_m^n f'''(x)P_3(x)\,dx \\
&=\frac{f(m)+f(n)}{2} + \frac{1}{6}\frac{f'(n) - f'(m)}{2!} - \frac{1}{30}\frac{f'''(n) - f'''(m)}{4!}-\int_m^n f^{(4)}(x) P_4(x)\, dx \\
&=\frac{f(m)+f(n)}{2} + \frac{1}{6}\frac{f'(n) - f'(m)}{2!} - \frac{1}{30}\frac{f'''(n) - f'''(m)}{4!} + \int_m^n f^{(5)}(x)P_5(x)\,dx \\
&=\frac{f(m)+f(n)}{2} + \frac{1}{6}\frac{f'(n) - f'(m)}{2!} - \frac{1}{30}\frac{f'''(n) - f'''(m)}{4!} + \frac{1}{42}\frac{f^{(5)}(n) - f^{(5)}(m)}{6!} - \int_m^n f^{(6)}(x)P_6(x)\,dx \\
&=\frac{f(m)+f(n)}{2} + \frac{1}{6}\frac{f'(n) - f'(m)}{2!} - \frac{1}{30}\frac{f'''(n) - f'''(m)}{4!} + \frac{1}{42}\frac{f^{(5)}(n) - f^{(5)}(m)}{6!} + \int_m^n f^{(7)}(x)P_7(x)\,dx.
\end{align}</math>

==Applications==

===The Basel problem===
The [[Basel problem]] is to determine the sum

:<math> 1 + \frac14 + \frac19 + \frac1{16} + \frac1{25} + \cdots = \sum_{n=1}^\infty \frac{1}{n^2}. </math>

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 <math>\pi^2/6</math>, which he proved in the same year.<ref>Pengelley, David J. [http://www.math.nmsu.edu/~davidp/euler2k2.pdf "Dances between continuous and discrete: Euler's summation formula"], in: Robert Bradley and Ed Sandifer (Eds), ''Proceedings, Euler 2K+2 Conference (Rumford, Maine, 2002)'', [[Euler Society]], 2003.</ref> 

===Sums involving a polynomial===
{{seealso|Faulhaber's formula}}
If <math>f</math> is a [[polynomial]] and <math>p</math> is big enough, then the remainder term vanishes.  For instance, if <math>f(x)=x^3</math>, we can choose <math>p=2</math> to obtain, after simplification,

:<math>\sum_{i=0}^n i^3 = \left(\frac{n(n + 1)}{2}\right)^2.</math>

===Approximation of integrals===
The formula provides a means of approximating a finite integral.  Let <math>a < b</math> be the endpoints of the interval of integration.  Fix <math>N</math>, the number of points to use in the approximation, and denote the corresponding step size by <math>h = (b - a) / (N - 1)</math>.  Set <math>x_i = a + (i-1)h</math>, so that <math>x_1 = a</math> and <math>x_N = b</math>.  Then:<ref name="Devries">{{cite book |last1=Devries |first1=Paul L. |last2=Hasbrun |first2=Javier E. |title=A first course in computational physics. |edition=2nd |publisher=Jones and Bartlett Publishers |year=2011 |page=156}}</ref>

:<math>
\begin{align}
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'(x_1) - f'(x_N)] - \frac{h^4}{720}[f'''(x_1) - f'''(x_N)] + \cdots.
\end{align}
</math>

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 <math>p</math>, depending upon <math>f</math> and <math>h</math>, such that the terms past order <math>p</math> increase rapidly.  Thus, the remainder term generally demands close attention.<ref name="Devries"/>

The Euler–Maclaurin formula is also used for detailed [[error analysis (mathematics)|error analysis]] in [[numerical quadrature]]. It explains the superior performance of the [[trapezoidal rule]] on smooth [[periodic function]]s and is used in certain [[Series acceleration|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 expansion]]s of sums and [[Series (mathematics)|series]], usually the most useful form of the Euler–Maclaurin formula is

:<math>\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)!} \left(f^{(2k - 1)}(b) - f^{(2k - 1)}(a)\right),</math>

where <math>a</math> and <math>b</math> are integers.<ref>{{harvtxt|Abramowitz|Stegun|1972}}, 23.1.30</ref> Often the expansion remains valid even after taking the limits <math>a\rightarrow-\infty</math> or <math>b\rightarrow +\infty</math> or both.  In many cases the integral on the right-hand side can be evaluated in [[Differential Galois theory|closed form]] in terms of [[elementary function]]s 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,

:<math>\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}}.</math>

Here the left-hand side is equal to <math>\psi^{(1)}(z)</math>, namely the first-order [[Polygamma function#Series representation|polygamma function]] defined by <math>\psi^{(1)}(z) = (d/dz)^2(\log \Gamma(z))</math>; the [[gamma function]] <math>\Gamma(z) </math> is equal to <math>(z-1)! </math> if <math>z </math> is a [[positive integer]]. This results in an asymptotic expansion for <math>\psi^{(1)}(z)</math>.  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 {{mvar|s}} is an integer greater than 1 we have:
:<math>\sum_{k=1}^n \frac{1}{k^s} \approx \frac 1{s-1}+\frac 12-\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].</math>

Collecting the constants into a value of the [[Riemann zeta function]], we can write an asymptotic expansion:
:<math>\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}}.</math>

For {{mvar|s}} equal to 2 this simplifies to
:<math>\sum_{k=1}^n \frac{1}{k^2} \sim\zeta(2)-\frac 1n+\frac 1{2n^2}-\sum_{i=1}\frac{B_{2i}}{n^{2i+1}},</math>
or
:<math>\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.</math>

When {{math|s {{=}} 1}}, the corresponding technique gives an asymptotic expansion for the [[harmonic number]]s:
:<math>\sum_{k=1}^n \frac{1}{k} \sim \gamma + \log n + \frac{1}{2n} - \sum_{k=1}^\infty \frac{B_{2k}}{2kn^{2k}},</math>
where <math>\gamma \approx 0.5772156649015328606065</math> is the [[Euler–Mascheroni constant]].

== Proofs ==

=== Derivation by mathematical induction ===
We outline the argument given in Apostol.<ref name=":0">{{Cite journal| first = T. M. | title = An Elementary View of Euler's Summation Formula | journal = The [[American Mathematical Monthly]] | volume = 106 | date = 1 May 1999 | issue = 5 | pages = 409–418  | issn = 0002-9890| last = Apostol | doi = 10.2307/2589145| authorlink= Tom M. Apostol| publisher =    Mathematical Association of America| jstor =    2589145}}</ref>

The [[Bernoulli polynomials]] {{math|''B<sub>n</sub>''(''x'')}} and the periodic Bernoulli functions {{math|''P<sub>n</sub>''(''x'')}} for {{math|''n'' {{=}} 0, 1, 2, ...}} were introduced above.

The first several Bernoulli polynomials are

:<math>\begin{align}
  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{align}</math>

The values {{math|''B<sub>n</sub>''(0)}} are the [[Bernoulli numbers]] {{math|''B''<sub>''n''</sub>}}.  Notice that for {{math|''n''&nbsp;≠&nbsp;1}} we have

:<math>B_n = B_n(0) = B_n(1),</math>

and for {{math|''n''&nbsp;{{=}}&nbsp;1}},

:<math>B_1 = B_1(0) = -B_1(1).</math>

The functions {{math|''P''<sub>''n''</sub>}} agree with the Bernoulli polynomials on the interval {{math|[0,&nbsp;1]}}  and are [[periodic function|periodic]] with period 1.  Furthermore, except when {{math|''n'' {{=}} 1}}, they are also continuous. Thus,

:<math> P_n(0) = P_n(1) = B_n \quad (n \neq 1)</math>

Let {{math|''k''}} be an integer, and consider the integral

:<math> \int_k^{k + 1} f(x)\,dx = \int_k^{k + 1} u\,dv,</math>

where

:<math>\begin{align}
   u &= f(x), \\
  du &= f'(x)\,dx, \\
  dv &= P_0(x)\,dx && (\text{since }P_0(x) = 1), \\
   v &= P_1(x).
\end{align}</math>

[[integration by parts|Integrating by parts]], we get

:<math>\begin{align}
 \int_k^{k + 1} f(x)\,dx &= \bigl[uv\bigr]_k^{k + 1} - \int_k^{k + 1} v\,du \\
                         &= \bigl[f(x)P_1(x)\bigr]_k^{k + 1} - \int_k^{k+1} f'(x)P_1(x)\,dx \\
                         &= B_1(1)f(k+1)-B_1(0)f(k) - \int_k^{k+1} f'(x)P_1(x)\,dx.
\end{align}</math>

Using <math>B_1(0)=-1/2</math>, <math>B_1(1)=1/2</math>, and summing the above from {{math|''k'' {{=}} 0}} to {{math|''k'' {{=}} ''n'' − 1}}, we get

:<math>\begin{align}
\int_0^n f(x)\, dx &= \int_0^1 f(x)\,dx + \dotsb + \int_{n-1}^n f(x)\,dx  \\
&= \frac{f(0)}{2}+ f(1) + \dotsb + f(n-1) + \frac{f(n)}{2} - \int_0^n f'(x) P_1(x)\,dx.
\end{align}</math>

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

:<math> \sum_{k=1}^n f(k) = \int_0^n f(x)\,dx + \frac{f(n) - f(0)}{2} + \int_0^n f'(x) P_1(x)\,dx.</math>

This is the {{math|''p'' {{=}} 1}} case of the summation formula.  To continue the induction, we apply integration by parts to the error term:

:<math>\int_k^{k+1} f'(x)P_1(x)\,dx = \int_k^{k + 1} u\,dv,</math>

where

:<math>\begin{align}
   u &= f'(x), \\
  du &= f''(x)\,dx, \\
  dv &= P_1(x)\,dx, \\
   v &= \frac{1}{2}P_2(x).
\end{align}</math>

The result of integrating by parts is

:<math>\begin{align}
  \bigl[uv\bigr]_k^{k + 1} - \int_k^{k + 1} v\,du &= \left[\frac{f'(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'(k + 1) - f'(k)) - \frac{1}{2}\int_k^{k + 1} f''(x)P_2(x)\,dx.
\end{align}</math>

Summing from {{math|''k'' {{=}} 0}} to {{math|''k'' {{=}} ''n'' &minus; 1}} and substituting this for the lower order error term results in the {{math|''p'' {{=}} 2}} case of the formula,

:<math>\sum_{k=1}^n f(k) = \int_0^n f(x)\,dx + \frac{f(n) - f(0)}{2} + \frac{B_2}{2}(f'(n) - f'(0)) - \frac{1}{2}\int_0^n f''(x)P_2(x)\,dx.</math>

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.

==See also==
* [[Cesàro summation]]
* [[Euler summation]]
* [[Gauss–Kronrod quadrature formula]]
* [[Darboux's formula]]
* [[Euler–Boole summation]]

==Notes==
{{reflist}}

==References==
{{refbegin|2}}
* {{Cite book | editor1-last=Abramowitz | editor1-first=Milton | editor1-link=Milton Abramowitz | editor2-last=Stegun | editor2-first=Irene A. | editor2-link=Irene Stegun | title=[[Abramowitz and Stegun|Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables]] | publisher=[[Dover Publications]] | location=New York | isbn=978-0-486-61272-0 | year=1972 | id=tenth printing | ref=harv }}, pp.&nbsp;16, 806, 886
* {{MathWorld|title=Euler–Maclaurin Integration Formulas|urlname=Euler-MaclaurinIntegrationFormulas}}
* Gourdon, Xavier; Sebah, Pascal ''[http://numbers.computation.free.fr/Constants/Miscellaneous/bernoulli.html Introduction on Bernoulli's numbers]'', (2002)
* {{cite book | first1=Hugh L. |last=Montgomery | authorlink=Hugh Montgomery (mathematician) |first2=Robert C. |last2=Vaughan |authorlink2=Robert Charles Vaughan (mathematician)  | title=Multiplicative number theory I. Classical theory | series=Cambridge tracts in advanced mathematics | volume=97 | year=2007 | isbn=0-521-84903-9 | pages=495–519}}
{{refend}}

{{DEFAULTSORT:Euler-Maclaurin Formula}}
[[Category:Approximation theory]]
[[Category:Asymptotic analysis]]
[[Category:Hilbert space]]
[[Category:Numerical integration (quadrature)]]
[[Category:Articles containing proofs]]
[[Category:Theorems in analysis]]
[[Category:Summability methods]]
[[Category:Leonhard Euler]]