Euler-Maclaurin method to integrate sums

Can we precisly estimate discrete sums using continuous integrals? Can we precisely estimate contnuous integrals using discrete sums? In some cases, the answer is yes. The Scottish mathematician Colin Maclaurin apparently discovered this technique around the same time that the Swiss genius Leonhard Euler did.

This is meant to be a concise summary of Euler's original result [1]. It took me some time to follow his derivation so it felt worthwhile to make a record of it, in notation that makes sense to me, focussed on the desired result. My interest in this topic started with using it to derive asymptotic series. After the derivation, one example of an asymptotic approximation and one example of an infinite sum are included.

We want to know, or approximate, the sum of a series with a known form. One way to write this sum is

$S={\Sigma }_{i=1}^{\mathrm{n}}f\left(x_{\mathrm{i}}\right)$

If the spacing between values of x is 1, then we can write the sum without the last term as

$S-f(x_{\mathrm{n}})=S-\frac{dS}{dx}+\frac{1}{2}\frac{d^{2}S}{d{x}^2}-\frac{1}{6}\frac{d^{3}S}{d{x}^3}+\dots$

This is just estimating the value of S(x-1) from S(x) via a Taylor series. It allows us to estimate the last term in the series from the sum:

$f(x_{\mathrm{n}})=\frac{dS}{dx}-\frac{1}{2}\frac{d^{2}S}{d{x}^2}+\frac{1}{6}\frac{d^{3}S}{d{x}^3}-\dots$

What we actually want is an estimate of the sum. There is a natural analogy between a sum and an integral, so that the rate of change of the sum should be proportional to the value of the last term in the series. To this end, write

$\frac{dS}{dx}=\alpha f_{\mathrm{n}}+\beta f_{\mathrm{n}}^{\mathrm{\text{'}}}+\gamma f_{\mathrm{n}}^{\mathrm{\text{'}\text{'}}}+\delta f_{\mathrm{n}}^{\mathrm{(3)}}+\epsilon f_{\mathrm{n}}^{\mathrm{(4)}}+\zeta f_{\mathrm{n}}^{\mathrm{(5)}}+\eta f_{\mathrm{n}}^{\mathrm{(6)}}+\theta f_{\mathrm{n}}^{\mathrm{(7)}}+\iota f_{\mathrm{n}}^{\mathrm{(8)}}+\dots$,

where the subscript n denotes evaluation at x_n and the superscripts denote derivatives. Now we integrate this to find the sum:

$S=\alpha {\int }_{x_1}^{x_{\mathrm{n}}}\mathrm{dx}f\mathrm{(x)}+\beta [f\mathrm{(x)}{]}_{x_1}^{x_{\mathrm{n}}}+\gamma [f^{\mathrm{\text{'}}}\mathrm{(x)}{]}_{x_1}^{x_{\mathrm{n}}}+\delta [f^{\mathrm{\text{'}\text{'}}}\mathrm{(x)}{]}_{x_1}^{x_{\mathrm{n}}}+\epsilon [f^{\mathrm{(3)}}\mathrm{(x)}{]}_{x_1}^{x_{\mathrm{n}}}+\dots$

Later, when we apply the formula, we will need to consider an integration constant. Now we substitute the previous series for f(x) in terms of S into this series. For example, for the first term, we get

$S=\alpha [S-\frac{1}{2}\frac{dS}{dx}+\dots {]}_{x_1}^{x_{\mathrm{n}}}+\dots$

When we substitute the series in the following terms, none of them will have the sum without a derivative. This tells us that α must be 1. The expansion for S and its derivatives will be consistent if all of the derivatives are set to zero with appropriate values of the Greek letters. Only the α and β terms have the first derivative of S. When we require the first derivative to vanish, we get

$\beta =\alpha /2$

Pursuing the next derivative, we find

$0=\alpha /6-\beta /2+\gamma ,$

or

$\gamma =-\alpha /6+\beta /2$

There is a pattern here. The next derivative (third derivative of S) gives

$\delta =\alpha /\mathrm{24}-\beta /6+\gamma /2$

Now we write at once the last five that we will pursue here.

$\epsilon =-\alpha /\mathrm{120}+\beta /\mathrm{24}-\gamma /6+\delta /2$

$\zeta =\alpha /\mathrm{720}-\beta /\mathrm{120}+\gamma /\mathrm{24}-\delta /6+\epsilon /2$

$\eta =-\alpha /\mathrm{5040}+\beta /\mathrm{720}-\gamma /\mathrm{120}+\delta /\mathrm{24}-\epsilon /6+\zeta /2$

$\theta =\alpha /8!-\beta /7!+\gamma /6!-\delta /5!+\epsilon /4!-\zeta /6+\eta /2$

$\iota =-\alpha /9!+\beta /8!-\gamma /7!+\delta /6!-\epsilon /5!+\zeta /4!-\eta /6+\theta /2$

This should be more than enough terms to demonstrate the pattern of the series. Now here are numerical values for these coefficients just for completeness.

$\alpha =1$

$\beta =1/2$

$\gamma =1/\mathrm{12}$

$\delta =0$

$\epsilon =-1/\mathrm{720}$

$\zeta =0$

$\eta =1/\mathrm{30240}$

$\theta =0$

$\iota =-1/\mathrm{1209600}$

Euler does give explicit formulae for more terms in this series, but this is plenty of terms for many applications. The last one listed here, ι, corresponds to the seventh derivative. The coefficient corresponding to the ninth derivative is 1/47900160. Now we can use these coefficients to write an estimate for a sum in terms of a functional form.

$S={\int }^{x_{\mathrm{n}}}\mathrm{dx}f\mathrm{(x)}+\frac{f(x_{\mathrm{n}})}{2}+\frac{1}{\mathrm{12}}\frac{\mathrm{df}}{\mathrm{dx}}(x_{\mathrm{n}})-\frac{1}{\mathrm{720}}{f}^{\mathrm{(3)}}(x_{\mathrm{n}})+\frac{1}{\mathrm{30240}}{f}^{\mathrm{(5)}}(x_{\mathrm{n}})+\dots$   (higher order terms)

We have dropped the lower integration limit because we will determine an integration constant in any case.

Applications of the formula

Now we can apply this formula. First, consider the logarithm of n!, which is the sum of the logarithm of all the positive integers up to n. If we take just the first two terms, then the result is

$logn!\approx C+nlogn-n+\frac{1}{2}logn$

If we set the value of C to half of the logarithm of 2π, then these first two terms give Stirling's approximation. The accuracy of a higher order application can be seen in the last plot on the factorial estimates page.

Next we can use the formula for a famous sum, which is also included in the original paper. Consider

$S={\Sigma }_{\mathrm{n=1}}^{\mathrm{\infty }}\frac{1}{n^2}$

We can apply the formula to get

$S\approx C-\frac{1}{n}+\frac{1}{2n^2}-\frac{1}{6n^3}+\frac{1}{\mathrm{30}{n}^5}-\frac{1}{\mathrm{42}{n}^7}$

If we simply calculate the first ten terms of the sum, by hand or computer, the answer is 1.54976773. In the formula for S including powers through 1/n⁷, the right-hand side evaluated at 10 (without C) gives -0.09516634. If we set C to the sum minus this right-hand side evaluation, the answer is 1.64493407. This is the estimate of the infinite sum, since everything but C vanishes as n approaches infinity.

Euler presented this technique in October of 1735. In December of the same year, Euler proved that this sum is equal to π²/6 [2], which to this number of decimals is 1.64493407. Euler used the integral technique outlined here to calculate the sum beyond double floating point precision on a modern computer, before proving that its actual value is π²/6. (He gives 21 digits in his value!)

Just to show the power of this technique, if we only compute the terms through n=4, then the sum is 1.42361111. The right-hand side without C gives -0.22132307. The solution for C this time is 1.64493418. So calculating only through n=4 in the series, then including contributions up to the fifth derivative in the formula, gives the correct infinite sum to more than 7 decimal places.

[1] "Finding the sum of any series for a given general term," Leonhard Euler, originally Inventio summae cuiusque seriei ex dato termino generali presented to St Petersburg Academy in October, 1735
[2] "On the sums of series of reciprocals," Leonhard Euler, originally De summis serierum reciprocarum presented to St Petersburg Academy in December, 1735

• Mathematics and physics
• Ski Sickness
• Patient history - chronological index