Fact-checked by Grok 2 weeks ago

Euler–Maclaurin formula

The Euler–Maclaurin formula is a fundamental result in mathematical analysis that provides an asymptotic expansion relating the sum of a smooth function f over consecutive integers to an integral of f, augmented by correction terms involving the function's derivatives at the endpoints and the Bernoulli numbers.^[1] Developed independently by Leonhard Euler in 1732 and Colin Maclaurin in 1742, the formula expresses \sum_{k=a}^{b} f(k) = \int_a^b f(x) \, dx + \frac{f(a) + f(b)}{2} + \sum_{m=1}^p \frac{B_{2m}}{(2m)!} \left( f^{(2m-1)}(b) - f^{(2m-1)}(a) \right) + R_p, where B_{2m} are the even Bernoulli numbers (with B_0 = 1, B_2 = 1/6, B_4 = -1/30, etc.) and R_p is a remainder term that can be controlled for sufficiently smooth f.^[2] The Bernoulli numbers, central to the formula, arise from the generating function \frac{x}{e^x - 1} = \sum_{n=0}^\infty B_n \frac{x^n}{n!}, and the formula's validity requires f to be sufficiently differentiable.^[1] This summation technique bridges discrete and continuous mathematics, enabling precise approximations for large sums by converting them into more tractable integrals while accounting for boundary effects through the endpoint corrections.^[3] For polynomials, the formula yields exact finite sums, as higher derivatives vanish beyond the degree, recovering identities like \sum_{k=1}^n k^3 = \left( \frac{n(n+1)}{2} \right)^2.^[2] In asymptotic analysis, it facilitates expansions for divergent series via smoothing, such as interpreting \sum_{n=1}^\infty n = -\frac{1}{12} in regularized senses relevant to physics and number theory.^[3] Key applications include deriving Stirling's approximation for factorials, n! \sim \sqrt{2\pi n} \left( \frac{n}{e} \right)^n as n \to \infty, by applying the formula to \ln n!; estimating the Euler-Mascheroni constant \gamma = \lim_{n \to \infty} \left( \sum_{k=1}^n \frac{1}{k} - \ln n \right) \approx 0.57721; and providing analytic continuations for the Riemann zeta function \zeta(s) for \Re(s) > 0 via smoothed sums.^[1] The formula also extends to more general settings, such as sums over lattice points in polytopes or in the context of Poisson summation, underscoring its versatility in analytic number theory and computational mathematics.^[3]

Introduction

Historical background

The Euler–Maclaurin formula has roots in earlier 17th-century work on finite differences and interpolation for summation. James Gregory developed a general interpolation formula for equidistant data in 1670, which relates sums to differences and laid groundwork for approximating series through polynomial interpolation.^[4] Independently, Isaac Newton advanced these ideas in the 1670s, formulating divided difference interpolation for arbitrary intervals and applying finite differences to approximate partial sums of series in his Principia (1687).^[4] These contributions provided initial motivations for connecting discrete sums to continuous integrals via differences. Leonhard Euler discovered the formula around 1732 during his studies on approximating slowly converging infinite series, motivated by the need to sum progressions more efficiently.^[5] He first outlined it in correspondence with James Stirling in 1736 but published the full result in 1738 in the Commentarii Academiae Scientiarum Petropolitanae, where it appeared as a general method for summing series using integrals and derivatives.^[5] Euler's approach built on finite difference techniques, extending them to higher-order corrections for better approximations of sums like those involving reciprocals. Independently, Colin Maclaurin derived a similar formula around 1742, publishing it in his Treatise of Fluxions, where he applied it to numerical computation of definite integrals through fluxional calculus.^[5] Maclaurin's work emphasized geometric and integral aspects, complementing Euler's series-focused perspective, and the formula became known by their names due to these parallel discoveries in the mid-18th century.^[5]

Overview and significance

The Euler–Maclaurin formula serves as a powerful bridge between discrete sums and continuous integrals, expressing a finite sum of a smooth function over integer points as its integral plus a series of correction terms involving the function's derivatives evaluated at the endpoints and coefficients known as Bernoulli numbers. This conceptual framework allows mathematicians and scientists to approximate complex summations by leveraging the more tractable tools of calculus, refining the basic trapezoidal rule approximation through these endpoint adjustments. The formula's structure highlights how sums can be viewed as Riemann sums with systematic error corrections, providing an intuitive way to understand the interplay between discrete and continuous mathematics.^[1] In numerical analysis and asymptotic methods, the formula holds significant importance by enabling high-accuracy approximations of sums without exhaustive computation, particularly for large ranges where direct evaluation is inefficient. It underpins the development of advanced quadrature techniques and asymptotic expansions, such as those for special functions like the gamma function via Stirling's approximation, and facilitates the analysis of divergent series by revealing their regularized behaviors. This has broad implications for deriving identities in analytic number theory and optimizing computational estimates in scientific simulations.^[1]^[6] The formula's relevance extends to theoretical physics, where it approximates partition functions in statistical mechanics by converting sums over energy levels into manageable integrals, aiding calculations of thermodynamic properties for systems like ideal gases or oscillators. In quantum mechanics, it contributes to semiclassical trace formulas, such as those in quantum chaos, by relating quantum spectral densities to classical periodic orbits through integral approximations that quantify chaotic enhancements in accuracy. Modern applications in computer science leverage the formula for efficient summation algorithms in numerical libraries and high-performance computing, supporting tasks in data analysis and simulation where precise integral-sum conversions are essential.^[7]^[8]^[1]

Mathematical formulation

Statement of the formula

The Euler–Maclaurin formula relates the sum of a function over integers to its integral plus correction terms derived from the function's derivatives at the endpoints. For integers a \leq b and a function f that is sufficiently smooth, the formula states

\sum_{k=a}^{b} f(k) = \int_{a}^{b} f(x) \, dx + \frac{f(a) + f(b)}{2} + \sum_{k=1}^{m} \frac{B_{2k}}{(2k)!} \left( f^{(2k-1)}(b) - f^{(2k-1)}(a) \right) + R_{m},

where R_{m} denotes the remainder term.^[9] Here, f is assumed to be C^{2m} (continuously differentiable up to order $2m) on the closed interval [a, b], ensuring the existence of the required derivatives.^[9] The notation B_{2k} refers to the Bernoulli numbers, a sequence of rational numbers defined by the generating function \frac{t}{e^t - 1} = \sum_{k=0}^{\infty} B_k \frac{t^k}{k!}, with B_{2k+1} = 0 for all k \geq 1, and f^{(n)} denotes the nth derivative of f.^[10] In cases where the remainder R_m can be made negligible or is omitted, the formula extends to an infinite series form:

\sum_{k=a}^{b} f(k) = \int_{a}^{b} f(x) \, dx + \frac{f(a) + f(b)}{2} + \sum_{k=1}^{\infty} \frac{B_{2k}}{(2k)!} \left( f^{(2k-1)}(b) - f^{(2k-1)}(a) \right).

This series is typically divergent for entire functions due to the rapid factorial-like growth of the Bernoulli numbers, specifically |B_{2k}| \sim 2 (2k)! / (2\pi)^{2k}, rendering it an asymptotic expansion rather than a convergent series.^[3]^[10]

Remainder term

The remainder term in the Euler–Maclaurin formula quantifies the error incurred when truncating the expansion after a finite number of correction terms involving Bernoulli numbers. For the approximation up to the term with B_{2m}, the remainder R_{m+1}(f; a, b) takes the explicit integral form

R_{m+1}(f; a, b) = \int_a^b \frac{\bar{B}_{2m+1}(x)}{(2m+1)!} f^{(2m+1)}(x) \, dx,

where \bar{B}_{2m+1}(x) = B_{2m+1}(\{x\}) denotes the periodic Bernoulli polynomial of degree $2m+1, obtained by periodically extending the Bernoulli polynomial B_{2m+1}(x) with period 1, and \{x\} = x - \lfloor x \rfloor is the fractional part of x.^[1]^[11] This expression arises naturally from repeated applications of integration by parts in the derivation, with the periodic extension ensuring the formula holds for non-integer limits a and b.^[1] To assess the magnitude of R_{m+1} on a finite interval [a, b], a standard bound is

|R_{m+1}(f; a, b)| \leq \frac{b - a}{(2m+1)!} \cdot \max_{x \in [a, b]} |f^{(2m+1)}(x)| \cdot \max_{y \in [0, 1]} |\bar{B}_{2m+1}(y)|.

The supremum norm of the periodic Bernoulli polynomial satisfies \max_{y \in [0, 1]} |\bar{B}_{2m+1}(y)| \leq 2 \zeta(2m+1) / (2\pi)^{2m+1}, where \zeta(s) is the Riemann zeta function; this estimate follows from the Fourier series expansion of \bar{B}_{2m+1}(y).^[11] For practical purposes, especially when boundary terms vanish (e.g., for integer a and b), further integration by parts shifts the derivative order, yielding an alternative bound

|R_{m+1}(f; a, b)| \leq \frac{(b - a)^{2m+2}}{(2m+2)!} \cdot \max_{x \in [a, b]} |f^{(2m+2)}(x)|,

which leverages a crude estimate via Taylor's theorem applied to f^{(2m+1)}(x) within the interval.^[1] As m increases, the asymptotic behavior of the remainder depends on the decay of the higher derivatives of f. If |f^{(k)}(x)| decreases sufficiently rapidly with k (faster than the growth of (k!) and the Bernoulli norms), then |R_{m+1}| diminishes, allowing for improved approximations with higher truncation orders.^[1] However, due to the factorial growth in the denominators and the super-exponential growth of the underlying Bernoulli numbers (|B_{2m}| \sim 4 \sqrt{\pi m} (m / ( \pi e))^{2m}), the full series is typically divergent and serves as an asymptotic expansion rather than a convergent one.^[1] Convergence of the infinite Euler–Maclaurin series occurs under restrictive conditions on f, such as when f is analytic in a horizontal strip |\Im z| < d containing the real interval [a, b], with suitable decay of f and its derivatives at infinity (e.g., f(x) = O(|x|^{-s}) for some s > 1 as |x| \to \infty). In such cases, the remainder can be analyzed via contour integration or resurgence techniques, ensuring the series sums to the exact difference between the sum and integral.^[12] For typical smooth functions without such analytic continuation, the optimal truncation point balances the decreasing early terms against the eventual divergence, minimizing the remainder.^[12]

Special cases

Low-order approximations

The zeroth-order approximation in the Euler–Maclaurin formula replaces the sum \sum_{k=a}^{b} f(k) with the integral \int_{a}^{b} f(x) \, dx, providing a basic estimate that overlooks endpoint contributions and higher derivatives, akin to a precursor of the rectangle rule in numerical integration.^[13] The first-order approximation improves this by adding the endpoint correction \frac{f(a) + f(b)}{2}, yielding \sum_{k=a}^{b} f(k) \approx \int_{a}^{b} f(x) \, dx + \frac{f(a) + f(b)}{2}. This form corresponds directly to the trapezoidal rule for approximating integrals, where the sum over function values at integer points serves as a discrete analog of the continuous integral with linear interpolation at the boundaries.^[13]^[1] For the second-order approximation, the formula incorporates the next correction term involving the second Bernoulli number B_2 = \frac{1}{6}: \sum_{k=a}^{b} f(k) \approx \int_{a}^{b} f(x) \, dx + \frac{f(a) + f(b)}{2} + \frac{B_2}{2!} \left( f'(b) - f'(a) \right) = \int_{a}^{b} f(x) \, dx + \frac{f(a) + f(b)}{2} + \frac{1}{12} \left( f'(b) - f'(a) \right). This adjustment accounts for the linear variation in the function's derivative, enhancing accuracy for smoother functions.^[13]^[14] Truncating the Euler–Maclaurin expansion at low orders introduces an error primarily captured by the remainder term, which for practical numerical integration behaves asymptotically like the first omitted correction; for instance, stopping at the first order yields an error on the order of the second-order term \frac{1}{12} (f'(b) - f'(a)), while second-order truncation shifts the leading error to higher derivatives involving B_4. This truncation strategy is particularly useful in computations where higher derivatives are costly to evaluate, balancing precision and efficiency in approximating sums by integrals.^[13] A representative example is the approximation of the nth harmonic number H_n = \sum_{k=1}^{n} \frac{1}{k}, where the low-order Euler–Maclaurin expansion gives H_n \approx \ln n + \gamma + \frac{1}{2n} - \frac{1}{12n^2}, with \gamma \approx 0.57721 denoting the Euler–Mascheroni constant; here, the zeroth-order term is \ln n + \gamma, the first-order adds \frac{1}{2n}, and the second-order subtracts \frac{1}{12n^2}, demonstrating rapid convergence for large n in series summation tasks.^[14]

Involvement of Bernoulli numbers

The Bernoulli numbers B_n appear in the Euler–Maclaurin formula as coefficients in the asymptotic expansion that corrects the integral approximation of sums, specifically in the terms involving even derivatives of the function.^[15] They are defined through the exponential generating function

\frac{x}{e^x - 1} = \sum_{n=0}^\infty B_n \frac{x^n}{n!},

valid for |x| < 2\pi, where the numbers are rational with B_0 = 1.^[14] The convention adopted in this context sets B_1 = -\frac{1}{2}, ensuring consistency with the formula's standard form, though an alternative B_1 = +\frac{1}{2} appears in some historical derivations.^[15] A key property facilitating their computation is the recurrence relation

\sum_{j=0}^m \binom{m+1}{j} B_j = 0

for all integers m \geq 1, which allows solving sequentially for higher B_m once lower terms are known.^[15] This relation stems directly from the generating function and underpins efficient recursive evaluation.^[16] In the Euler–Maclaurin formula, only even-indexed Bernoulli numbers contribute to the expansion beyond the first-order term, because B_{2k+1} = 0 for all integers k \geq 1.^[14] This vanishing of odd-indexed terms (except B_1) simplifies the series to involve solely B_{2k}, reflecting the symmetry in the periodic extension of the Bernoulli polynomials underlying the formula.^[15] Explicit values for small indices illustrate their role; for instance, B_2 = \frac{1}{6} and B_4 = -\frac{1}{30}, with signs alternating as (-1)^{k-1} B_{2k} > 0 for k \geq 1.^[14] These can be computed via the recurrence, starting from B_0 and B_1.^[15] The asymptotic growth of the Bernoulli numbers, given by

|B_{2n}| \sim 4 \sqrt{\pi n} \left( \frac{n}{\pi e} \right)^{2n}

as n \to \infty, combined with the (2k)! denominators in the formula's terms, renders the series divergent but useful as an asymptotic expansion, where truncating at optimal order minimizes error.^[17] This factorial growth ensures the terms eventually decrease before increasing, a hallmark of such approximations.^[15]

Proofs

Derivation via integration by parts

To derive the Euler–Maclaurin formula, begin by considering the sum \sum_{k=a}^{b} f(k) where a and b are integers with a < b, and f is a sufficiently smooth function. Express the integral \int_a^b f(x) \, dx as a sum of integrals over unit intervals: \int_a^b f(x) \, dx = \sum_{k=a}^{b-1} \int_k^{k+1} f(x) \, dx.^[18] The difference between the sum and the integral arises from approximating each integral \int_k^{k+1} f(x) \, dx by values of f at the endpoints. A precise first-order relation is the trapezoidal rule, obtained via integration by parts. Let F(x) be an antiderivative of f(x), so F'(x) = f(x). For each interval, use a change of variable s = x - k for x \in [k, k+1], so \int_k^{k+1} f(x) \, dx = \int_0^1 f(s + k) \, ds. Integrate by parts with respect to s using the function P_1(s) = s - \frac{1}{2}, which satisfies \int_0^1 P_1(s) \, ds = 0 and P_1'(s) = 1.^[1]^[18] Specifically, \int_0^1 f(s + k) \, ds = \int_0^1 P_1'(s) F(s + k) \, ds = [P_1(s) F(s + k)]_0^1 - \int_0^1 P_1(s) f'(s + k) \, ds. Here, P_1(1) = \frac{1}{2} and P_1(0) = -\frac{1}{2}, so the boundary term is P_1(1) F(k+1) - P_1(0) F(k) = \frac{1}{2} F(k+1) + \frac{1}{2} F(k) = \frac{f(k) + f(k+1)}{2}, since F(k+1) - F(k) = \int_k^{k+1} f(x) \, dx but the evaluation gives the average directly in the limit for point values. Rearranging yields the trapezoidal approximation:

\int_k^{k+1} f(x) \, dx = \frac{f(k) + f(k+1)}{2} + \int_k^{k+1} P_1(x - k) f'(x) \, dx.

Summing over k from a to b-1, the trapezoidal terms telescope: \sum_{k=a}^{b-1} \frac{f(k) + f(k+1)}{2} = \sum_{k=a}^b f(k) - \frac{f(a) + f(b)}{2}. Thus,

\int_a^b f(x) \, dx = \sum_{k=a}^b f(k) - \frac{f(a) + f(b)}{2} + \int_a^b \bar{P}_1(x) f'(x) \, dx,

where \bar{P}_1(x) = P_1(\{x\}) is the periodic extension of P_1 with period 1, and \{x\} = x - \lfloor x \rfloor. Rearranging gives:

\sum_{k=a}^{b} f(k) = \int_a^b f(x) \, dx + \frac{f(a) + f(b)}{2} - \int_a^b \bar{P}_1(x) f'(x) \, dx.

This establishes the initial correction (noting the sign convention for the remainder).^[18] To refine this, repeat integration by parts on the remainder integral \int_a^b \bar{P}_1(x) f'(x) \, dx. Introduce the next periodic function \bar{P}_2(x) such that \bar{P}_2'(x) = \bar{P}_1(x), defined by \bar{P}_2(x) = \int_0^x \bar{P}_1(t) \, dt + c_2, where the constant c_2 ensures periodicity: \int_0^1 \bar{P}_2(x) \, dx = 0. Iteratively, define \bar{P}_{m+1}(x) = \int_0^x \bar{P}_m(t) \, dt + c_{m+1} for m \geq 1, with c_{m+1} chosen for zero mean over [0,1]. These \bar{P}_m(x) are related to the periodic Bernoulli polynomials \bar{B}_m(x) = B_m(\{x\}), where B_m(y) are the Bernoulli polynomials satisfying B_m'(y) = m B_{m-1}(y) and B_m(y+1) - B_m(y) = m y^{m-1}. Specifically, \bar{P}_m(x) = \frac{\bar{B}_m(x)}{m!}.^[1]^[18] Applying integration by parts m times yields endpoint contributions from the derivatives of f. The boundary terms accumulate as differences of odd-order derivatives at the endpoints, scaled by the values of the Bernoulli polynomials at integers. Since \bar{B}_{2k}(n) = B_{2k} (the Bernoulli numbers) for integer n, and odd Bernoulli numbers B_{2k+1} = 0 for k \geq 1, only even-indexed terms survive in the series. After m steps, the formula becomes:

\sum_{k=a}^{b} f(k) = \int_a^b f(x) \, dx + \frac{f(b) + f(a)}{2} + \sum_{j=1}^m \frac{B_{2j}}{(2j)!} \left[ f^{(2j-1)}(b) - f^{(2j-1)}(a) \right] + R_m,

where the remainder R_m = (-1)^{m+1} \int_a^b \frac{\bar{B}_{2m+1}(x)}{(2m+1)!} f^{(2m+1)}(x) \, dx (adjusting sign for consistency). The emergence of Bernoulli numbers B_{2j} stems from the Taylor expansion of the periodic Bernoulli polynomials around integers, where \bar{B}_{2j}(x) = B_{2j} + higher terms vanishing at integers. As m \to \infty, under suitable decay conditions on f^{(2m+1)}, the series converges to the exact relation.^[1]^[18]

Proof by mathematical induction

The Euler–Maclaurin formula can be proved by mathematical induction on the order of the approximation, particularly in the context of finite differences, where the formula provides an exact representation for polynomials of sufficiently low degree. The proof leverages the forward difference operator to establish the inductive step, connecting the discrete sum to the integral through correction terms involving Bernoulli numbers. This approach is especially suited for exact finite sums, as higher-order differences vanish for polynomials. The forward difference operator is defined as \Delta f(x) = f(x+1) - f(x), with higher-order differences given by \Delta^n f(x) = \sum_{k=0}^n (-1)^{n-k} \binom{n}{k} f(x+k). These operators capture the discrete analog of derivatives and play a key role in the proof, as the inductive step involves applying \Delta to both sides of the assumed formula to generate the next order. The Bernoulli numbers enter through their generating function \frac{t}{e^t - 1} = \sum_{k=0}^\infty \frac{B_k t^k}{k!}, which provides the expansion for the antidifference operator \Delta^{-1}, linking finite differences to the continuous integral and endpoint corrections.^[19]^[20] For the base case with order m=0, the formula reduces to the trapezoidal approximation:

\sum_{n=a+1}^b f(n) \approx \int_a^b f(t) \, dt + \frac{f(a) + f(b)}{2}.

This can be verified directly by considering the identity for the first Bernoulli polynomial B_1(x) = x - \frac{1}{2}, where the endpoint average corrects the integral for the linear interpolation error over each interval. For a linear function f(t) = t, the sum and integral match exactly with the endpoint term, confirming the base case holds without remainder.^[19] Assume the formula holds for order $2k-1:

\sum_{n=a+1}^b f(n) = \int_a^b f(t) \, dt + \sum_{j=1}^{k} \frac{B_{2j}}{(2j)!} \left[ f^{(2j-1)}(b) - f^{(2j-1)}(a) \right] + R_{2k-1},

where the remainder R_{2k-1} involves the next Bernoulli polynomial. To prove the step for order $2k+1, apply the forward difference operator \Delta to both sides. The left side becomes \Delta \sum f(n) = f(b+1) - f(a+1), while the right side transforms as \Delta \int f = \int \Delta f + boundary adjustments from the Leibniz rule for differences. The endpoint terms telescope under \Delta, and the correction terms shift to produce the next Bernoulli coefficient via the relation \Delta B_m(x) = m x^{m-1}, which holds for Bernoulli polynomials and induces the inductive extension. The generating function ensures the coefficients align with B_{2k+2}/(2k+2)!, completing the step. By induction, the formula holds for all finite orders.^[20]^[19] For polynomials of degree d < 2k, the proof terminates exactly, as higher-order forward differences \Delta^{d+1} f(x) = 0 and corresponding derivatives f^{(d+1)} = 0, making the remainder zero and yielding a closed-form sum. This discrete inductive structure complements continuous derivations, such as integration by parts, by emphasizing finite difference properties for exact polynomial evaluations.^[19]

Applications

Integral approximations of sums

The Euler–Maclaurin formula provides a systematic way to approximate finite sums by integrals, particularly useful for smooth functions f over large intervals. For a sum \sum_{k=1}^n f(k), the basic approximation is \sum_{k=1}^n f(k) \approx \int_1^n f(x) \, dx + \frac{f(1) + f(n)}{2}, where the endpoint correction \frac{f(1) + f(n)}{2} accounts for the trapezoidal rule's adjustment over the crude rectangular or midpoint Riemann sums.^[1] Higher-order terms involving Bernoulli numbers and derivatives of f can then be added to refine this estimate, yielding \sum_{k=1}^n f(k) \approx \int_1^n f(x) \, dx + \frac{f(1) + f(n)}{2} + \sum_{k=1}^m \frac{B_{2k}}{(2k)!} \left( f^{(2k-1)}(n) - f^{(2k-1)}(1) \right) + R, with a remainder R that decreases as more terms are included.^[21] This approach is especially effective for large n, where the integral dominates and the corrections capture boundary effects efficiently.^[6] Incorporating additional terms from the formula significantly reduces the approximation error compared to basic Riemann sums. The zeroth-order integral alone often over- or underestimates the sum by an amount proportional to the function's variation, but the first correction term halves this error for linear functions and improves accuracy for higher smoothness.^[1] Each subsequent Bernoulli term eliminates contributions from higher even powers in the asymptotic error expansion, leading to exponential convergence in the number of terms for analytic functions, though practical limits arise from derivative computation costs.^[21] In numerical practice, truncating after a few terms balances accuracy and efficiency, outperforming equidistant quadrature rules for non-periodic integrands.^[6] A classic example is the approximation of the harmonic numbers H_n = \sum_{k=1}^n \frac{1}{k}, where applying the Euler–Maclaurin formula yields H_n \approx \ln n + \gamma + \frac{1}{2n} - \sum_{k=1}^m \frac{B_{2k}}{2k n^{2k}} + R, with \gamma \approx 0.57721 the Euler–Mascheroni constant emerging as the constant term in the expansion.^[1] For n=10, the first two terms give H_{10} \approx 2.92897, while the full approximation with up to B_4 achieves error below $10^{-5}, demonstrating rapid improvement over the simple \ln n + 1 estimate.^[22] This illustrates how the formula transforms a discrete sum into a computable integral plus asymptotic series, ideal for large n where direct summation is feasible but analytical insight is desired.^[1] The Euler–Maclaurin formula serves as the theoretical foundation for advanced numerical integration techniques that approximate sums via integrals. Romberg integration, for instance, applies Richardson extrapolation to trapezoidal rules derived from the formula's leading terms, achieving higher-order accuracy without explicit derivative evaluations.^[23] It also informs enhancements to Gaussian quadrature by providing error estimates for interpolatory rules on finite intervals, allowing adaptive node placement to minimize remainder contributions.^[24] These methods leverage the formula's structure to extend basic sum-to-integral conversions into robust algorithms for scientific computing.^[25]

Solution to the Basel problem

The Euler–Maclaurin formula can be applied to the Basel problem, which seeks the exact value of \zeta(2) = \sum_{k=1}^\infty \frac{1}{k^2}. Applying the formula to the partial sums \sum_{k=1}^n f(k) with f(x) = x^{-2} yields an asymptotic expansion that approximates the sum to high precision for large n, with terms involving Bernoulli numbers providing corrections to the integral \int_1^n x^{-2} \, dx = 1 - 1/n. As n \to \infty, the endpoint and derivative terms at the upper limit vanish, leaving contributions from the lower limit that relate to the Bernoulli numbers.^[1] This asymptotic series allows numerical evaluation of \zeta(2) by truncating at optimal terms, as the divergent nature of the full series is managed by stopping before divergence sets in. Euler employed an early version of the formula around 1735 to compute the sum numerically to 20 decimal places, verifying his exact result \zeta(2) = \frac{\pi^2}{6} obtained via the infinite product expansion of the sine function.^[1] The Euler–Maclaurin formula connects to the exact evaluation through its role in generating expressions involving Bernoulli numbers. Euler discovered the general formula for even positive integers: \zeta(2m) = (-1)^{m+1} \frac{B_{2m} (2\pi)^{2m}}{2 (2m)!}, where the Bernoulli numbers B_{2m} (with B_2 = \frac{1}{6}, B_4 = -\frac{1}{30}, etc.) provide the rational coefficients relating the sums to powers of \pi. This underscores the deep connection between discrete sums, integrals, and trigonometric functions.^[1]

Asymptotic expansions

The Euler–Maclaurin formula yields asymptotic expansions for sums \sum_{k=1}^n f(k) as n \to \infty, particularly when the higher derivatives of f decay sufficiently rapidly, allowing the series to provide corrections to the leading integral approximation.^[1] In such cases, the formula expresses the sum as an integral plus endpoint corrections and a series involving Bernoulli numbers applied to the derivatives at the upper limit, with contributions from the lower limit absorbed into constants or neglected for large n.^[1] The general asymptotic form is

\sum_{k=1}^n f(k) \sim \int_1^n f(x) \, dx + \frac{f(1) + f(n)}{2} + \sum_{k=1}^m \frac{B_{2k}}{(2k)!} f^{(2k-1)}(n) + R_m,

where B_{2k} are Bernoulli numbers and the remainder R_m is controlled by the next term in the expansion; terms at the lower limit 1 are fixed and contribute to the overall constant.^[1] This expansion is especially useful for functions f where the integral provides the dominant behavior, and the Bernoulli series refines it with inverse powers of n. A prominent application is the derivation of Stirling's approximation for the factorial n!. Applying the formula to f(x) = \log x, the sum \sum_{k=1}^n \log k = \log(n!) becomes

\log(n!) \approx \int_1^n \log x \, dx + \frac{\log n}{2} + \sum_{k=1}^m \frac{B_{2k}}{(2k)!} \frac{(2k-2)!}{n^{2k-1}} + C + R_m,

where \int_1^n \log x \, dx = n \log n - n + 1, the constant C = \frac{1}{2} \log(2\pi) arises from the remainder analysis as n \to \infty, and the series terms yield corrections like \frac{1}{12n} - \frac{1}{360n^3} + \cdots.^[1] Exponentiating gives the leading asymptotic n! \sim \sqrt{2\pi n} \, (n/e)^n, with higher-order terms from the divergent series improving accuracy for finite n.^[1] The formula extends naturally to the gamma function \Gamma(z+1), providing Stirling's series in the complex plane. For large |z| with |\arg z| < \pi - \delta (\delta > 0),

\log \Gamma(z+1) \sim \left(z + \frac{1}{2}\right) \log z - z + \frac{1}{2} \log(2\pi) + \sum_{k=1}^m \frac{B_{2k}}{(2k) \, z^{2k-1}} + R_m,

derived by applying Euler–Maclaurin to an integral representation or the digamma function's expansion, leading to \Gamma(z+1) \sim \sqrt{2\pi z} \, (z/e)^z.^[26] This holds in specified sectors away from the negative real axis, where the asymptotic validity is ensured by the decay of derivatives.^[26] Since the Bernoulli numbers grow factorially (|B_{2k}| \sim 4 \sqrt{\pi k} (k / (2\pi e))^{2k}), the series in the Euler–Maclaurin expansion is typically divergent, but it serves as an asymptotic series where partial sums up to the smallest term provide optimal approximation before the remainder grows.^[1] Truncation at the term minimizing the error—often around m \approx |z|/2 for the gamma case—yields exponentially small remainders relative to the leading term.^[26]

Sums of polynomials

When the function f(x) in the Euler–Maclaurin formula is a polynomial of degree p, the summation series terminates exactly after the (p+1)-th term because all higher-order derivatives of f vanish. This yields a closed-form expression for finite sums \sum_{k=1}^n f(k) without remainder, transforming the approximate relation into an equality.^[27] For the specific case f(x) = x^p where p is a positive integer, the formula simplifies to

\sum_{k=1}^n k^p = \frac{1}{p+1} \sum_{j=0}^p (-1)^j \binom{p+1}{j} B_j n^{p+1-j},

with Bernoulli numbers B_j (using the convention B_1 = -1/2) and B_j = 0 for odd j > 1. This expression, known as Faulhaber's formula in its Bernoulli number form, expresses the sum as a polynomial in n of degree p+1. Equivalently, it can be written as

\sum_{k=1}^n k^m = \frac{1}{m+1} \left[ n^{m+1} + \frac{1}{2} n^m + \sum_{k=1}^{\lfloor m/2 \rfloor} \frac{B_{2k}}{2k} \binom{m+1}{2k} n^{m+1-2k} \right],

highlighting the contribution from even-indexed Bernoulli numbers.^[28]^[29] A simple illustration is the first-order case p=1: truncating after the B_1 term gives

\sum_{k=1}^n k = \frac{1}{2} \left[ n^2 + (-1) \binom{2}{1} \left(-\frac{1}{2}\right) n \right] = \frac{n(n+1)}{2},

recovering the well-known arithmetic series sum. For higher powers, such as p=3, the full formula yields \sum_{k=1}^n k^3 = \left( \frac{n(n+1)}{2} \right)^2.^[27] The formula extends to sums over arbitrary integer intervals [a, b] using Bernoulli polynomials B_j(x), where

\sum_{k=a}^b k^p = \frac{B_{p+1}(b+1) - B_{p+1}(a)}{p+1},

with B_m(x) defined by the generating function \frac{t e^{xt}}{e^t - 1} = \sum_{m=0}^\infty B_m(x) \frac{t^m}{m!} and satisfying B_m(0) = B_m for the Bernoulli numbers. This generalization preserves the exactness for polynomials.^[1]

Generalizations

Extensions to higher dimensions

The Euler–Maclaurin formula extends to higher dimensions by generalizing the summation over lattice points in \mathbb{Z}^d to approximate integrals over \mathbb{R}^d with corrections involving boundary terms derived from multivariate Bernoulli polynomials or differential operators associated with the faces of the domain. For a smooth function f on \mathbb{R}^d, the multidimensional version takes the form

\sum_{k \in \mathbb{Z}^d} f(k) \approx \int_{\mathbb{R}^d} f(x) \, dx + \sum_{\text{boundary terms}} \text{corrections},

where the corrections incorporate higher-order contributions from the boundaries, often expressed using periodized multivariate Bernoulli polynomials \mathfrak{B}_{J,L}(x) indexed by multi-indices J and linear transformations L. These terms generalize the one-dimensional Bernoulli numbers to tensors or polynomial structures that account for the geometry of the lattice. This extension, developed for regular summands, has been further refined to handle functions with algebraic singularities via singular Euler–Maclaurin expansions, where operator coefficients are computed using derivatives of the Epstein zeta function.^[30]^[31] A specific case arises for sums over rectangular domains, such as double sums \sum_{i=1}^m \sum_{j=1}^n g(i,j), which approximate the double integral \iint g(x,y) \, dx \, dy plus edge integrals along the boundaries and adjustments at the corners. The formula includes line corrections analogous to one-dimensional boundary terms on each edge, combined with corner contributions that prevent overcounting, yielding

\sum_{i=1}^m \sum_{j=1}^n g(i,j) = \int_1^m \int_1^n g(x,y) \, dx \, dy + \sum_{\text{edges}} \int_{\text{edge}} g \, ds + \sum_{\text{corners}} g(c) + \text{higher-order terms}.

This structure arises from iterative application of the one-dimensional formula or direct multidimensional derivations, ensuring consistency for separable functions while capturing mixed partial derivative effects on non-rectangular boundaries.^[31] Applications of these extensions include counting lattice points in polyhedra, where the formula relates the number of integer points in a dilated polytope tP \cap \mathbb{Z}^d to its volume and lower-dimensional face contributions, directly informing the coefficients of Ehrhart polynomials. For a convex polytope P, the Ehrhart quasipolynomial is expressed as

\sum_{x \in tP \cap \mathbb{Z}^d} 1 = \sum_{f \in F(P)} \nu_0(P, f, t) \cdot \vol(f) \cdot t^{\dim f},

with periodic coefficients \nu_0 determined locally by the transverse geometry at each face f, providing an exact asymptotic expansion without remainder for polynomial indicators. However, challenges emerge in higher dimensions, as the remainder term involves increasingly complex infinite-order differential operators, and its estimation grows with d, necessitating decomposition algorithms like Barvinok's for computational feasibility in fixed dimensions.

Relation to other summation formulas

The Euler–Maclaurin formula shares a deep connection with the Poisson summation formula through the lens of Fourier analysis, where the former serves as a low-frequency approximation to the latter for sufficiently smooth functions. The Poisson summation formula states that for a suitable function f, \sum_{n=-\infty}^{\infty} f(n) = \sum_{k=-\infty}^{\infty} \hat{f}(k), with \hat{f} denoting the Fourier transform; expanding \hat{f}(k) in a Taylor series around k=0 yields the Euler–Maclaurin expansion as the leading terms, capturing the integral and endpoint corrections while higher k terms provide the remainder.^[32] This link is particularly evident in applications to spectral computations, such as on even-dimensional spheres, where Euler–Maclaurin approximates sums over polynomial eigenvalues in place of the full Poisson formula when Fourier multiplicities align with low-order polynomials.^[32] Closely related is the Abel–Plana formula, a complex-analytic extension that overlaps with Euler–Maclaurin in its correction terms but incorporates contour integrals for broader applicability to holomorphic functions. The Abel–Plana formula expresses

\sum_{n=0}^{\infty} f(n) = \int_{0}^{\infty} f(x) \, dx + \frac{1}{2} f(0) + i \int_{0}^{\infty} \frac{f(it) - f(-it)}{e^{2\pi t} - 1} \, dt + R,

where R is a remainder involving higher derivatives, providing an exact representation under suitable decay and analyticity conditions.^[13] Unlike Euler–Maclaurin, which relies on real-variable derivatives and Bernoulli polynomials for finite sums, Abel–Plana uses the cotangent or its variants in the complex plane, allowing evaluation of Euler–Maclaurin constants via residue calculus; the two formulas are interderivable through Cauchy's integral theorem and Weierstrass approximations for periodic extensions.^[33] Their corrections overlap in the endpoint and Bernoulli terms, making Abel–Plana a complementary tool for infinite sums where Euler–Maclaurin remainders diverge.^[13] Both the Euler–Maclaurin and Poisson summation formulas employ Bernoulli numbers in their generating functions, facilitating zeta regularization for divergent series. The generating function \frac{t}{e^t - 1} = \sum_{k=0}^{\infty} \frac{B_k t^k}{k!} underlies the Bernoulli corrections in Euler–Maclaurin, which, when applied to partial sums of \zeta(s), yield analytic continuations like \zeta(-m) = -\frac{B_{m+1}}{m+1} for positive integers m, regularizing expressions such as $1 + 2 + 3 + \cdots = -\frac{1}{12}.^[3] Poisson summation similarly leverages this via Fourier periodicity, connecting zeta values to functional equations, though Euler–Maclaurin provides the real-analytic bridge for low-order regularization.^[3] Conceptually, Euler–Maclaurin emphasizes real-analytic expansions along the line, suitable for asymptotic approximations of sums via derivatives, whereas Poisson summation operates in the spectral domain through Fourier duality, offering exact equalities for periodic or rapidly decaying functions.^[33] They are complementary for entire functions: Euler–Maclaurin excels in local Taylor-like corrections, while Poisson captures global harmonic structure; Abel–Plana bridges them by complexifying the real corrections of Euler–Maclaurin.^[13] Historically, Leonhard Euler's foundational work on series expansions in the 1730s directly influenced both the Poisson summation formula, developed by Siméon Denis Poisson in 1823 as a Fourier-based generalization of Euler's ideas, and the Abel–Plana formula, independently discovered by Niels Henrik Abel and Giovanni Plana around 1820–1823, which extended Euler's integral-sum relations into the complex plane.^[33] Euler's manipulations of infinite series, including early forms of Bernoulli number applications, laid the groundwork for these interconnections, as recognized in later analyses of their derivations from shared contour and periodicity principles.^[33]

References

[1]
[PDF] Euler-Maclaurin Formula 1 Introduction - People | MIT CSAIL
It is believed to be irrational and even transcendental by many mathematicians but there is no proof known up to date. Let consider a more general type of ...
[2]
[PDF] Euler-Maclaurin summation formula, PHYS 2400
Apr 25, 2024 · 4 The Euler-Maclaurin summation formula. Equation (24) is the Euler-Maclaurin summation formula. It can be rewritten for the case of a finite ...<|control11|><|separator|>
[3]
The Euler-Maclaurin formula, Bernoulli numbers, the zeta function ...
Apr 10, 2010 · Here is interesting result which one can use to come up with Euler -Maclaurin formula.I found out that 1/(1+X) is convergent for X>1 if we ...
[4]
[PDF] A Chronology of Interpolation: From Ancient Astronomy to Modern ...
Abstract—This paper presents a chronological overview of the developments in interpolation theory, from the earliest times to the present date.
[5]
https://encyclopediaofmath.org/wiki/Euler-MacLaurin_formula
[6]
[PDF] Euler-Maclaurin expansions without analytic derivatives
Aug 14, 2020 · The Euler-Maclaurin formulas relate sums and integrals. Discovered nearly 300 years ago, they have lost none of their importance over the ...
[7]
[PDF] Statistical Mechanics - Oberlin College and Conservatory
Aug 29, 2008 · The Euler-Maclaurin formula states that if f(x) → 0, f0(x) → 0, f00(x) → 0, etc. as x → ∞, then. ∞. X. `=0 f(`) ≈. Z ∞. 0 f(x) dx + 1. 2 f ...
[8]
[PDF] Accuracy of trace formulas
Abstract. Using quantum maps, we study the accuracy of semiclassical trace formulas. The role of chaos in improving the semiclassical accuracy in some ...
[9]
[PDF] 5. The Euler-Maclaurin Summation Formula
Mar 10, 2011 · With this function, we can state a first form of the Euler-Maclaurin summation formula. This formula shows how a sum can be approximated by an ...
[10]
[PDF] The Bernoulli Numbers: A Brief Primer - Whitman College
May 10, 2019 · Bernoulli numbers, and the Euler-Maclaurin Summation Formula are expressed. We hope that this will help the adventurous undergraduate reader ...
[11]
[PDF] Two Proofs of Euler-Maclaurin Formula, Its Generalizations and ...
Maclaurin applied it to the numerical computation of definite integrals, while Euler used it to calculate series. It has an extensive application in many ...
[12]
[PDF] resurgence of the euler-maclaurin summation formula
Abstract. The Euler-MacLaurin summation formula compares the sum of a function over the lattice points of an interval with its corresponding integral, ...<|separator|>
[13]
DLMF: §2.10 Sums and Sequences ‣ Areas ‣ Chapter 2 Asymptotic ...
The asymptotic behavior of entire functions defined by Maclaurin series can be approached by converting the sum into a contour integral by use of the residue ...
[14]
DLMF: §24.2 Definitions and Generating Functions ‣ Properties ‣ Chapter 24 Bernoulli and Euler Polynomials
### Summary of Asymptotic Expansion for Harmonic Numbers (Using Euler-Maclaurin)
[15]
Bernoulli Number -- from Wolfram MathWorld
Bernoulli numbers are a sequence of signed rational numbers important in number theory and analysis, arising in series expansions of trigonometric functions.Missing: factorial | Show results with:factorial
[16]
DLMF: §24.4 Basic Properties ‣ Properties ‣ Chapter 24 Bernoulli ...
For the relation of Bernoulli numbers to the Riemann zeta function see §25.6, and to the Eulerian numbers see (26.14.11). 24.3 Graphs24.5 Recurrence Relations.
[17]
DLMF: §24.9 Inequalities ‣ Properties ‣ Chapter 24 Bernoulli and ...
§24.9 Inequalities. ⓘ. Keywords: Bernoulli numbers, Bernoulli polynomials, Euler numbers, Euler polynomials, inequalities; Notes: For (24.9.1) and (24.9.2) ...<|separator|>
[18]
[PDF] Analysis of Series and Products. Part 1: The Euler–Maclaurin Formula
Aug 11, 2020 · This is the first order Euler–Maclaurin formula. Further formulae are obtained using integration by parts, differentiating f0 and integrating P1 ...
[19]
[PDF] understanding the euler-maclaurin formula
Thus he later derived the Euler-Maclaurin formula in 1732, which was subsequently derived by Maclaurin in 1742.Missing: primary sources 1738
[20]
[1705.09159] Approximating sums by integrals only - arXiv
May 24, 2017 · The Euler--Maclaurin (EM) summation formula is used in many theoretical studies and numerical calculations. It approximates the sum \sum_{k=0}^{ ...
[21]
[PDF] The Euler constant : γ
In fact, using the harmonic number notation ... A similar estimation is given in [6]. 2.2 Asymptotic expansion of the harmonic numbers. The Euler-Maclaurin ...
[22]
[PDF] Romberg Integration - Numerical Analysis
Oct 21, 2022 · Romberg Integration. 1. Adaptive Integration stopping the composite trapezoidal rule a Julia function. 2. The Euler-Maclaurin Summation Formula.
[23]
[PDF] quadrature.pdf - UMD MATH
... Euler-Maclaurin formula. This means that all corollaries that we have drawn ... • the Gaussian quadrature is exact for ˜p0, ˜p1, ..., ˜pn−1,. • P(xj) ...
[24]
[PDF] 9. Approximation of Integrals – Continued 9.1. Iterative Approaches ...
Richardson Extrapolation and Romberg Integration. ... Euler-Maclaurin summation formula (refr ... From the Euler-Maclaurin summation formula,. L(f) ...
[25]
None
### Summary of Stirling's Series Derivation for Gamma Function Using Euler-Maclaurin Formula
[26]
[PDF] Euler-Maclaurin summation formula, Physics 2400
Mar 29, 2016 · 4 The Euler-Maclaurin summation formula. Equation (24) is the Euler-Maclaurin summation formula. It can be rewritten for the case of a finite.
[27]
[PDF] An Extended Version of Faulhaber's Formula
Apr 7, 2016 · This paper presents an extended version of the well-known Faulhaber formula, which is used to compute the sum of the m-th powers of the first n ...<|control11|><|separator|>
[28]
[PDF] #A78 INTEGERS 22 (2022) FAULHABER'S FORMULA, ODD ...
Aug 3, 2022 · In order to write an expression for P km in n, we introduce the Bernoulli numbers. Set B0 = 1 and define Bn by n. X k=0 n + 1 k. Bk = 0, where ...
[29]
Singular Euler-Maclaurin expansion on multidimensional lattices
Feb 22, 2021 · First, the Euler-Maclaurin summation formula is generalised to lattices in higher dimensions, assuming a sufficiently regular summand function.
[30]
https://arxiv.org/abs/2102.10941
[31]
https://doi.org/10.1007/s00041-023-10011-9
[32]
(PDF) The summation formulae of Poisson, Plana, Euler-Maclaurin ...
Nov 9, 2015 · Besides the history of these formulas, several of their modifications and generalizations are considered. Connections between the Euler– ...Missing: ties | Show results with:ties