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 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 m and n are natural numbers and f(x) is a real or complex valued continuous function for real numbers x in the interval [m,n], then the integral
I = ∫ m n f ( x ) d x {\displaystyle I=\int _{m}^{n}f(x)\,dx}
can be approximated by the sum (or vice versa)
S = f ( m + 1 ) + ⋯ + f ( n − 1 ) + f ( n ) {\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 f(k)(x) evaluated at the endpoints of the interval, that is to say x = m and x = n. Explicitly, for p a positive integer and a function f(x) that is p times continuously differentiable on the interval [m,n], we have
S − I = ∑ k = 1 p B k k ! ( f ( k − 1 ) ( n ) − f ( k − 1 ) ( m ) ) + R p , {\displaystyle S-I=\sum _{k=1}^{p}{{\frac {B_{k}}{k!}}\left(f^{(k-1)}(n)-f^{(k-1)}(m)\right)}+R_{p},}
where Bk is the kth Bernoulli number (with B1 = 1/2) and Rp is an error term which depends on n, m, p, and f and is usually small for suitable values of p. The formula is often written with the subscript taking only even values, since the odd Bernoulli numbers are zero except for B1. In this case we have
∑ i = m n f ( i ) = ∫ m n f ( x ) d x + f ( n ) + f ( m ) 2 + ∑ k = 1 ⌊ p 2 ⌋ B 2 k ( 2 k ) ! ( f ( 2 k − 1 ) ( n ) − f ( 2 k − 1 ) ( m ) ) + R p , {\displaystyle \sum _{i=m}^{n}f(i)=\int _{m}^{n}f(x)\,dx+{\frac {f(n)+f(m)}{2}}+\sum _{k=1}^{\left\lfloor {\frac {p}{2}}\right\rfloor }{\frac {B_{2k}}{(2k)!}}\left(f^{(2k-1)}(n)-f^{(2k-1)}(m)\right)+R_{p},}
or alternatively
∑ i = m + 1 n f ( i ) = ∫ m n f ( x ) d x + f ( n ) − f ( m ) 2 + ∑ k = 1 ⌊ p 2 ⌋ B 2 k ( 2 k ) ! ( f ( 2 k − 1 ) ( n ) − f ( 2 k − 1 ) ( m ) ) + R p . {\displaystyle \sum _{i=m+1}^{n}f(i)=\int _{m}^{n}f(x)\,dx+{\frac {f(n)-f(m)}{2}}+\sum _{k=1}^{\left\lfloor {\frac {p}{2}}\right\rfloor }{\frac {B_{2k}}{(2k)!}}\left(f^{(2k-1)}(n)-f^{(2k-1)}(m)\right)+R_{p}.}
The remainder term
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 [r, r + 1] for r = m, m + 1, …, n − 1. 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 Pk(x). The Bernoulli polynomials may be defined recursively by B0(x) = 1 and, for k ≥ 1,
B k ′ ( x ) = k B k − 1 ( x ) , ∫ 0 1 B k ( x ) d x = 0. {\displaystyle B_{k}'(x)=kB_{k-1}(x),\qquad \int _{0}^{1}B_{k}(x)\,dx=0.}
The periodized Bernoulli functions are defined as
P k ( x ) = B k ( x − ⌊ x ⌋ ) , {\displaystyle P_{k}(x)=B_{k}{\bigl (}x-\lfloor x\rfloor {\bigr )},}
where ⌊x⌋ denotes the largest integer less than or equal to x, so that x − ⌊x⌋ always lies in the interval [0,1). With this notation, the remainder term Rp equals
R p = ( − 1 ) p + 1 ∫ m n f ( p ) ( x ) P p ( x ) p ! d x . {\displaystyle R_{p}=(-1)^{p+1}\int _{m}^{n}f^{(p)}(x){\frac {P_{p}(x)}{p!}}\,dx.}
When k > 0, it can be shown that for 0 ≤ x ≤ 1,
| B k ( x ) | ≤ 2 ⋅ k ! ( 2 π ) k ζ ( k ) , {\displaystyle {\bigl |}B_{k}(x){\bigr |}\leq {\frac {2\cdot k!}{(2\pi )^{k}}}\zeta (k),}
where ζ denotes the Riemann zeta function; one approach to prove this inequality is to obtain the Fourier series for the polynomials Bk(x). The bound is achieved for even k when x is zero. The term ζ(k) may be omitted for odd k but the proof in this case is more complex (see Lehmer). Using this inequality, the size of the remainder term can be estimated as
| R p | ≤ 2 ζ ( p ) ( 2 π ) p ∫ m n | f ( p ) ( x ) | d x . {\displaystyle \left|R_{p}\right|\leq {\frac {2\zeta (p)}{(2\pi )^{p}}}\int _{m}^{n}\left|f^{(p)}(x)\right|\,dx.}
Low-order cases The Bernoulli numbers from B1 to B7 are 1/2, 1/6, 0, −1/30, 0, 1/42, 0. Therefore, the low-order cases of the Euler–Maclaurin formula are:
∑ i = m n f ( i ) − ∫ m n f ( x ) d x = f ( m ) + f ( n ) 2 + ∫ m n f ′ ( x ) P 1 ( x ) d x = f ( m ) + f ( n ) 2 + 1 6 f ′ ( n ) − f ′ ( m ) 2 ! − ∫ m n f ″ ( x ) P 2 ( x ) 2 ! d x = f ( m ) + f ( n ) 2 + 1 6 f ′ ( n ) − f ′ ( m ) 2 ! + ∫ m n f ‴ ( x ) P 3 ( x ) 3 ! d x = f ( m ) + f ( n ) 2 + 1 6 f ′ ( n ) − f ′ ( m ) 2 ! − 1 30 f ‴ ( n ) − f ‴ ( m ) 4 ! − ∫ m n f ( 4 ) ( x ) P 4 ( x ) 4 ! d x = f ( m ) + f ( n ) 2 + 1 6 f ′ ( n ) − f ′ ( m ) 2 ! − 1 30 f ‴ ( n ) − f ‴ ( m ) 4 ! + ∫ m n f ( 5 ) ( x ) P 5 ( x ) 5 ! d x = f ( m ) + f ( n ) 2 + 1 6 f ′ ( n ) − f ′ ( m ) 2 ! − 1 30 f ‴ ( n ) − f ‴ ( m ) 4 ! + 1 42 f ( 5 ) ( n ) − f ( 5 ) ( m ) 6 ! − ∫ m n f ( 6 ) ( x ) P 6 ( x ) 6 ! d x = f ( m ) + f ( n ) 2 + 1 6 f ′ ( n ) − f ′ ( m ) 2 ! − 1 30 f ‴ ( n ) − f ‴ ( m ) 4 ! + 1 42 f ( 5 ) ( n ) − f ( 5 ) ( m ) 6 ! + ∫ m n f ( 7 ) ( x ) P 7 ( x ) 7 ! d x . {\displaystyle {\begin{aligned}\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){\frac {P_{2}(x)}{2!}}\,dx\\&={\frac {f(m)+f(n)}{2}}+{\frac {1}{6}}{\frac {f'(n)-f'(m)}{2!}}+\int _{m}^{n}f'''(x){\frac {P_{3}(x)}{3!}}\,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){\frac {P_{4}(x)}{4!}}\,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){\frac {P_{5}(x)}{5!}}\,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){\frac {P_{6}(x)}{6!}}\,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){\frac {P_{7}(x)}{7!}}\,dx.\end{aligned}}}
Applications
The Basel problem The Basel problem is to determine the sum
1 + 1 4 + 1 9 + 1 16 + 1 25 + ⋯ = ∑ n = 1 ∞ 1 n 2 . {\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 π2/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) = x3, we can choose p = 2 to obtain, after simplification,
∑ i = 0 n i 3 = ( n ( n + 1 ) 2 ) 2 . {\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 xi = a + (i − 1)h, so that x1 = a and xN = b. Then:
I = ∫ a b f ( x ) d x ∼ h ( f ( x 1 ) 2 + f ( x 2 ) + ⋯ + f ( x N − 1 ) + f ( x N ) 2 ) + h 2 12 [ f ′ ( x 1 ) − f ′ ( x N ) ] − h 4 720 [ f ‴ ( x 1 ) − f ‴ ( x N ) ] + ⋯ {\displaystyle {\begin{aligned}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}}{\bigl [}f'(x_{1})-f'(x_{N}){\bigr ]}-{\frac {h^{4}}{720}}{\bigl [}f'''(x_{1})-f'''(x_{N}){\bigr ]}+\cdots \end{aligned}}}
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
∑ n = a b f ( n ) ∼ ∫ a b f ( x ) d x + f ( b ) + f ( a