François-Joseph Servois (1767–1847) was a French priest, artillery officer, professor of mathematics, and museum curator. Servois conducted research in many areas of mathematics, including geometry and the differential calculus. In 1817, after a three-year period with no publications, Servois published his "Memoir on Quadratures," where he tackled a new area of mathematics: numerical integration. This paper examines the "Quadratures" paper, in which Servois addressed a debate on numerical integration techniques among the mathematicians Christian Kramp (1760–1826), Joseph-Diez Gergonne (1771–1859), and Joseph-Balthazard Bérard (1763–1844?). Our paper provides a readers' guide to Servois' 1817 paper and concludes with several student activities that can be incorporated into a numerical analysis course. We also offer readers a complete English translation of Servois' "Memoir on Quadratures" (pdf).
François-Joseph Servois (1767–1847) was born on July 19, 1767, in the village of Mont-de-Laval, located in the French Department of Doubs, near the Swiss border. Servois pursued several careers during his life: he was a priest, artillery officer, professor of mathematics, and Curator of the Artillery Museum, located in the 7th Arrondissement of Paris. He died on April 17, 1847, in Mont-de-Laval, where "his legacy remains legendary" [Boyer 1895a, p. 313]. Readers interested in a more in-depth biography of Servois can refer to Bradley [2002], Petrilli [2010], and Petrilli [2017].
In a list of all mathematicians, Servois would rank as lesser-known. To the extent that his name is known at all, it is for introducing the words "distributive", "commutative", and "pole" into mathematics. However, Servois made significant contributions, in particular to the field of the Calculus of Derivations (or Operator Calculus). We feel his three major works were:
Solutions peu connues de différents problèmes de géométrie-pratique^{1}
"Essai sur un nouveau mode d’exposition des principes du calcul différentiel"^{2}
"Réflexions sur les divers systèmes d’exposition des principes du calcul différentiel, et, en particulier, sur la doctrine des infiniment petits^{"}^{3}
The first is Servois’ only book, which was published in 1804. This geometry text, which Servois compiled during his career in the military, was intended to be a reference on applied geometry for military officers, presenting constructions in ruler-geometry that could be used on the battlefield at any time [Bradley 2002]. This book brought him fame among specialists—indeed, the great mathematician Jean-Victor Poncelet (1788–1867) consulted Servois for his expertise on geometry several times during the writing of the 1822 edition of the Traité des propriétés projectives [Taton 1972].
Readers interested in Servois’ geometry textbook can refer to the English translation located in [Servois 1804]. As this stands as Servois’ first major work in mathematics, it exposes us to the type of mathematician that Servois would become. It is clear that military personnel would not be interested in the geometric justifications for these constructions, but rather would prefer the simplest possible steps for completing these tasks quickly. Nevertheless, Servois provided proofs of all of his applied constructions throughout the book. Later in his mathematical career, we see that Servois became a strong advocate of rigor^{4} in mathematics.
In 1814, Servois, following in the footsteps of Joseph-Louis Lagrange (1736–1813), proposed a foundation for calculus based on power series through the publication of his "Essay" [Servois 1814a]. Servois developed therein the known rules of the differential calculus using his definition of the differential:
\[dz = \Delta z - \frac{1}{2} \Delta^2 z + \frac{1}{3} \Delta^3 z - \ldots,\]
and the linear operator properties of the differential. He wrote a follow-up paper, the "Reflections" [Servois 1814b], as a philosophical justification of his work in the "Essa"' and an explanation of how basing the calculus on power series brought rigor to the subject. Readers interested in Servois' contributions to the foundations of calculus can refer to [Bradley and Petrilli 2010a] and [Bradley and Petrilli 2010b].
After a three-year break, Servois returned to research in calculus, but this time with a focus on numerical integration [Servois 1817]. From a psychological point of view, this “Memoir on Quadratures” illustrates Servois' return to the frame of mind that he was in when he wrote his geometry textbook. In this paper, Servois was not concerned with rigor and foundations, but rather with the practical applications of numerical integration. All the same, he provided derivations and proofs appropriate to this stage of the development of calculus. He also engaged in debates over the utility of divergent series, a topic that was of relevance to the subject.
In the next section, we summarize the dispute that motivated Servois to write his “Quadratures'” paper.
1. Little-known Solutions to Various Problems in Practical Geometry. Petrilli [2017] includes an English translation of this monograph.
2. "Essay on a New Method of Exposition of the Principles of Differential Calculus." Bradley and Petrilli [2010a] includes an English translation of this paper.
3. "Reflections on the Various Systems of Exposition of the Principles of the Differential Calculus and, in particular, on the Doctrine of the Infinitely Small." Bradley and Petrilli [2010b] includes an English translation of this paper.
4. In mathematics, the term rigor has many definitions. The most appropriate definition for describing Servois’ notion of rigor is provided by Grabiner [1981, p. 5]: First, every concept of the subject had to be explicitly defined in terms of concepts whose nature was held to be already known. ... Second, theorems had to be proved, with every step in the proof justified by a previously proved theorem, by a definition, or by an explicitly stated axiom. Third, the definitions chosen, and the theorems proved, had to be sufficiently broad to support the entire structure of valid results belonging to the subject.
Joseph-Diez Gergonne's (1771–1859) Annales de mathématiques pures et appliquées had its fair share of debates and disputes within its pages. Servois was involved in what was probably the most famous of these, regarding the use of geometry to represent complex numbers, along with Jean-Robert Argand (1768–1822) and Jacques Français (1775–1833).
Servois' “Quadratures” paper was inspired by another debate, which occurred from 1815 to 1817. Servois [1817] stated:
For my part, I confess that I read with great satisfaction, in the Annales de mathématiques, the detailed expositions of three new methods of approximation that came from good sources, because they belong to professors Dobenheim^{1}, Kramp, and Bérard. I have participated, with more or less the aptitude of an interested party, in the debates that they have initiated [p. 73].
In 1815 Christian Kramp (1760–1826) published a paper in Gergonne's Annales, in which he gave an approximation for the area under a curve using the Trapezoidal Rule and an ad hoc extrapolation method that was first proposed by Alexandre-Magnier (Magnus) d'Obenheim (1753–1840), who was a professor at various French artillery schools from 1801 to 1835.
In particular, Kramp provided approximations to the area using the Composite Trapezoidal Rule for \(n =\) 1, 2, 3, 4, 6, and 12 subdivisions. His extrapolation procedure used these cruder approximations to generate a much better approximation for the area under the curve. He showed that his method of approximation is very accurate when used to approximate the values of \(\ln 2\) and \(\frac{\pi}{4}\).
Shortly after Kramp's article was published, Gergonne wrote a follow-up piece stating that Kramp's idea was truly remarkable, but the approximation it produced was not very useful in a practical sense without an error term. In the examples of \(\ln 2\) and \(\frac{\pi}{4}\), these values are well-known, but there is no way to assess the accuracy of an approximation without an error bound when the true value is unknown. Gergonne argued that we should not adopt methods of numerical integration that sacrifice precision for the sake of speed.
Kramp wrote a follow-up paper to Gergonne's rebuttal. In this paper, Kramp abandoned his original idea, which was based on the work of d'Obenheim. Instead, he took on the problem of expressing a curve as a Taylor Series and integrating it term-by-term. He derived what are now known as the Newton-Cotes formulas for several cases; however, there was an error in his formula for calculating the approximation in the case of \(n=12\) subdivisions.
Joseph-Balthazard Bérard (1763–1844?) entered the debate in 1816. In his paper, Bérard discussed issues with the d'Obenheim-Kramp method of finding the area under a curve and the improved method that Kramp introduced in his second article. Bérard asserted that the method in Kramp's first article was arbitrary and that he only succeeded in making accurate approximations by chance. He noted that Gergonne's proposed improvement was just as arbitrary and gave less accurate approximations. In fact, he claimed that adapting Gergonne's method to \(n= 16\) subdivisions gave rise to a poorer approximate of \(\ln 2\) than the case of \(n=12\). Instead, he proposed to follow the methods of Kramp's second paper and gave much more efficient ways of calculating the coefficient formulas needed for the approximation. He was even able to calculate the formulas for \(n = 24\) and \(n = 120\) subdivisions. Along the way, he turned up an error in Kramp's formula for the case of \(n = 12\) subdivisions.
Finally, later in 1816, Kramp wrote a third paper that defended the method he introduced in his second paper and critiqued Bérard's paper. Kramp argued that his formula for \(n = 12\) subdivisions was the correct one, even though Bérard's formula gave a slightly more accurate approximate value of \(\ln 2\). There is little if any new mathematics in this paper, which is more of a polemical review of Bérard's work.
In his “Memoir on Quadratures,” Servois attempted to give an overview of numerical integration that encompassed the points of view expressed in these works of Kramp, Gergonne, and Bérard. Let us now explore Servois' insights into this debate. We begin by presenting an overview (in today's notation) of the numerical integration formulas discussed by Servois: the Trapezoidal Rule and the Newton-Cotes Formulas. We then provide a section-by-section thematic breakdown of Servois' “Quadratures” paper. We next provide introductions to the finite difference notation that Servois himself employed and to his general approach to questions in analysis, sometimes called the operational calculus. Finally, we examine two specific topics from Servois' memoir in more detail: his use of Bernoulli Numbers and his derivation of the Euler-MacLaurin formula, which served as his primary tool for resolving the debate on numerical integration.
1. In [Servois 1817], Servois used this alternate spelling of d'Obenheim.
Most students learn in calculus that a simple approximation to the definite integral is
\[\int^b_a f(x) \, dx \approx \frac{b-a}{2}[f(a) + f(b)].\]
This is called the Trapezoidal Rule, because the right-hand side is the area of the trapezoid bounded by the interval \([a,b]\) on the \(x\)-axis, the lines \(x=a\) and \(x=b\) and the straight line segment joining \((a,f(a))\) to \((b,f(b))\). Although this approximation generally isn't very accurate, it can be greatly improved by breaking the interval \([a,b]\) into subintervals.
In particular, we let \(n\) be a natural number, \(\omega = \frac{b-a}{n}\), and \(x_k = a +k \omega\), for \(k=0, 1, 2, \ldots, n\). Then if we apply the Trapezoid Rule to each subinterval \([x_{k-1},x_k]\), for \(k= 1, 2, \ldots, n\), and sum the results, we have^{1}
\[\int^b_a f(x) \, dx \approx \frac{\omega}{2}\left[f(x_0) +2 \sum^{n-1}_{k=1}f(x_k) + f(x_n)\right]. \,\,\,\,\,\,\,\,\,\,\,\,\,\,(I)\]
This is called the Composite Trapezoidal Rule, or just the Trapezoidal Rule in some books. Intuitively, the approximation should get better and better with larger values of \(n\). This is indeed the case for nice functions \(f\). For example, under the hypothesis that \(f''(x)\) is bounded on \([a,b]\), the error in the approximation is roughly proportional to \(n^{-2}\); see [Burden 2016, p. 205].
In [Kramp 1815a], the Composite Trapezoidal Rule is applied repeatedly for \(n=1\), 2, 3, 4, 6, and 12, and these six results are substituted into an extrapolating formula that gives a new approximation that is generally much better even than the approximation corresponding to \(n=12\). Kramp showed that this gave surprisingly good estimates (0.69314806 and 0.78539271 respectively) of
\[\ln 2 = \int^2_1 \frac{dx}{x} \quad \mbox{and}\quad \frac{\pi}{4} = \int^1_0 \frac{dx}{1+x^2}.\]
Gergonne criticized Kramp's result both because of the arbitrariness of the extrapolation rule and the fact that, without an error bound, we have no way of knowing how good the approximation is when the true value of the integral is unknown–and if the value is known, then there's no need for an approximation procedure!
1. Equation numbers in this article are given in Roman numerals to distinguish them from the original equation numbers used by Servois in [Servois 1817]. There is no correlation between these Roman numerals and his Hindu Arabic numerals.
Many modern calculus books also include Simpson's Rule. With the notation of the previous section and \(n = 2,\) Simpson's Rule says: \[\int^b_a f(x) \, dx \approx \frac{\omega}{3}[f(x_0) + 4f(x_1) + f(x_2)].\] This rule can be derived by fitting a parabola to the three points \((x_0,f(x_0))\), \((x_1,f(x_1))\) and \((x_2,f(x_2))\), and using the area under that parabola between \(x=a\) and \(x=b\) to approximate the area under \(y=f(x)\). The Trapezoidal Rule and Simpson's Rule are just the first two cases (\(n=1\) and 2, respectively) of a family of formulas known as the Closed Newton-Cotes Formulas. In [Kramp 1815b] and [Bérard 1816], the authors derive the Newton-Cotes formulas for many values of \(n\), and give different answers in the case of \(n=12\).
There is a unique polynomial of degree \(\leq n\) that passes through the \(n+1\) distinct points \[ (x_0,f(x_0)), \; (x_1,f(x_1)), \; \ldots, \; (x_n,f(x_n)).\] The Newton-Cotes formulas are derived by finding the polynomial in a useful form and integrating it to approximate the integral of \(f(x)\). The existence of such a polynomial was proven by Lagrange in 1795 and the form in which he gave it is called the Lagrange Interpolating Polynomial [Lagrange 1795, p. 276ff]. For the \(n+1\) given points, let \[L_{n,k} = \frac{\prod^n_{i=0, i\ne k}(x-x_i)}{\prod^n_{i=0, i\ne k}(x_k-x_i)},\] for \(k=0, 1,\ldots n\). Notice that \(L_{n,k}(x_j) = \delta_{j,k}\) in terms of Kronecker's delta function. Therefore \[P(x) = \sum_{k=0}^n f(x_k)L_{n,k}(x) \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\mbox{(II)}\] satisfies \(P(x_k) = f(x_k)\) for \(k=0, 1,\ldots n\).
The achievements of Kramp and Bérard in deriving Newton-Cotes formulas for large \(n\) may be impressive, but as a matter of practice, such formulas are almost never used. Not only are they inconvenient to use, but the polynomial (II) has an oscillatory nature that limits the accuracy of these formulas for large \(n\). Instead, modern numerical integration is often based on the composite version of Simpson's Rule inside an adaptive algorithm, such as Romberg Integration; see e.g. [Burden 2016, p. 207].
It's not clear if Kramp and Bérard knew of the work of Roger Cotes (1682–1716). He was an English mathematician, a fellow of Trinity College, Cambridge and the editor of the second edition of Newton's Principia. He published only one paper in his lifetime, but after his untimely death at the age of 33, his papers were published posthumously [Cotes 1722]. In the last part of this book (Opera Miscellanea), Cotes described how to derive the Newton-Cotes formulas, and on page 33 gave the specific formulas for \(n=3,4, \ldots 11\).
Servois' memoir is a long paper, some sections of which can be omitted without affecting understanding of the material that follows. Here we provide a thematic breakdown of the sections in Servois' paper, where page numbers refer to the pagination of the original article, which are given in square brackets in the English translation of Servois' "Memoir on Quadratures" (pdf).
When tackling questions regarding numerical integration in his "Memoir on Quadratures,'' Servois used slightly different notation than the current notation that we introduced earlier in this article. He considered a function \(F\), fixed a constant difference \(\Delta x=\omega\), and then defined a sequence of \(n+1\) points
\[(a,F(a)), (a+\omega,F(a+\omega)), (a + 2\omega,F(a + 2 \omega)), \ldots (x,F(x))\]
with equally-spaced \(x\)-coordinates (which Servois usually called "abscissas").
Although subscripts were not in general use during Servois' time, it will make our lives simpler to use them in discussing his memoir. Servois considered a sequence on the \(x\)-axis
\[x_k = a +k \omega; \; k=0, 1, 2, \ldots, n\]
and the corresponding \(y\)-values \(y_k = F(x_k)\). He denoted the left endpoint \((x_0,y_0)\) by \((a,v)\) and the right endpoint \((x_n,y_n)\) by \((x,y)\), which can be a little confusing, given that \(x\) and \(y\) are also the general names of the variables. Although he didn't use subscripts, Servois had access to a useful notation, due to the French mathematician Louis François Antoine Arbogast (1759–1803). In his influential textbook Du calcul des dérivations [Arbogast 1800], Arbogast defined an operator \(E\) called the varied state by the equation
\[E y = F(x + \Delta x),\]
which is \(F(x + \omega)\) in Servois' setting. The varied state can be iterated, so that^{1}
\[E^ny=F(x + n\omega),\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\mbox{(III)}\]
which is the definition given by Servois in his equation (1).
The varied state also has an inverse \(E^{-1}\) which can also be iterated
\[E^{-n}y = F(x-n\omega).\]
In fact, Servois even considered expressions like
\[E^{\frac{1}{2}}y = F\left(x + \frac{1}{2} \omega\right);\]
thus, in equation (III) we could even think of \(n\) as representing a continuous parameter.
The varied state is useful for the differential calculus, because it can be used in the definition of the finite difference \(\Delta y\):
\[\Delta y = F(x + \Delta x) - F(x)= Ey - y.\]
Because \(Ey - y\) may be thought of as \(Ey - 1y\), Servois wrote
\[\Delta y = (E-1)y\]
in his equation (2), and even abstracted this to the equation \(\Delta = E -1\), defining the operator \(\Delta\) in terms of the varied state and the identity operator. This rather forward-looking kind of abstraction was at the core of what was called the “Calculus of Derivations” or the “Operational Calculus,'” a field pioneered by Arbogast to which Servois made important contributions in [Servois 1814a].
1. Recall that equation numbers in this article are given in Roman numerals to distinguish them from the original equation numbers used by Servois in [Servois 1817]. There is no correlation between these Roman numerals and his Hindu Arabic numerals.
In Continental Europe before the 19th century, calculus was done in terms of differentials and not derivatives. For example, taking the differential of the equation \(y=x^2\) would give rise to \(dy = 2x \, dx\) by the rules of Gottfried Wilhelm Leibniz's (1646–1716) differential calculus. Today, we might formally divide this equation through by \(dx\) to get \(\frac{dy}{dx} = 2x\), the expression of the derivative in Leibniz' notation, but to mathematicians in Servois' time and earlier, the equation \(dy = 2x \, dx\) was meaningful in itself, giving the relationship between an "infinitely small" increment in the \(x\)-direction and the corresponding increment in the \(y\) direction—also infinitely small.
At the time Servois was writing his "Memoir on Quadratures,'' this point of view was gradually being replaced by the scheme we use today: finding the quotient of the finite differences \(\Delta y\) and \(\Delta x\), then taking the limit as \(\Delta x \rightarrow 0\). For more on the foundation of calculus in Servois' time, and his contributions to the search for a foundation, see [Bradley and Petrilli 2010b].
Instead of taking limits, Servois defined the differential \(dy\) to be the following alternating series:
\[dy = \Delta y - \frac{1}{2} \Delta^2 y + \frac{1}{3} \Delta^3 y - \ldots .\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \mbox{(IV)}\]
He gave a derivation of this series in [Servois 1814a]. We can get an idea of how it works by considering the following simple example.
Set \(y=x^2\). Then \(\Delta y = Ey - y = (x + \Delta x)^2 - x^2 = 2 x \Delta x + (\Delta x)^2\), and
\[\begin{align*}\Delta^2 y &= \Delta (\Delta y) \\&= \Delta \left(2 x \Delta x + (\Delta x)^2\right)\\&= \left[2(x + \Delta x)\Delta x + (\Delta x)^2\right] - \left[2 x \Delta x + (\Delta x)^2\right]\\&= 2(\Delta x)^2.\end{align*}\]
Furthermore, \(\Delta^3 y = 0\) and so \(\Delta^k y = 0\) for \(k \geq 3\). Substituting these values into (IV), we thus have
\[dy = 2 x \Delta x + (\Delta x)^2 -\frac{1}{2}\left(2(\Delta x)^2\right) = 2x \Delta x,\]
or \(dy = 2x \,dx\) upon interpreting \(dx\) for \(\Delta x\). This nominal equivalence between \(dx\) and \(\Delta x\) explains why Servois wrote \(\omega = \Delta x = dx\) just above his equation (1).
Servois, and a number of mathematicians before him, noticed that the series (IV) has the same form as the MacLaurin series for \(\ln(1 + x)\):
\[\ln(1+x) = x - \frac{x^2}{2} + \frac{x^3}{3} + \ldots.\]
Servois used Log to denote the natural logarithm and wrote \(dy = \mbox{Log}(1+ \Delta)y\) in his equation (3). Cancellation of the \(y\)’s in this equation then gives \(d = \mbox{Log}(1+ \Delta)\), so that Servois' equation (3) can be formally expressed as \(e^d y = (1 + \Delta)y = Ey\).
We note that Servois' equation (1) is incorrect as written, whether the error is typographical or otherwise. It should probably have been given as
\[E^n(y) = F(x+n\omega) \quad \mbox{and} \quad Ey = e^d y.\]
In his equation (4), Servois derived the inverse \(\sum\) of the finite difference operator \(\Delta\).
To find an expression for \(\sum y = \Delta^{-1} y = (E -1)^{-1}y\), we observe that by the definition of inverse, we must have \((E-1)(E-1)^{-1}y = y\). We then check that
\[\begin{align*}(E-1) & [E^{-1}y + E^{-2}y + E^{-3}y + \cdots] \\& = (y - E^{-1}y) + (E^{-1}y - E^{-2}y) + (E^{-2}y - E^{-3}y) + \cdots,\end{align*}\]
which is a telescoping sum, formally equal to \(y\); this gives us a candidate for \(\sum y\).
There is also the observation that if \(K\) is any quantity at all satisfying \(\Delta K=0\), then \((E-1)K=0\).
Such a quantity \(K\) plays the role for \(\sum\) that is analogous to the arbitrary constant \(C\) in an indefinite integral. Therefore,
\[\sum y = E^{-1} y + E^{-2} y + E^{-3} y + \cdots + K,\]
a formal sum that we might write as
\[\begin{align*}\sum^{\infty}_{k=1} f(x - k\omega) + K&= \sum^{\infty}_{k=1} f(x_0 + n\omega - k\omega) + K \\&= \sum^{\infty}_{k=n+1} f(x_0 + (n-k)\omega) + \sum^{n}_{k=1} f(x_0 + (n-k)\omega) + K \\ &= \sum^{\infty}_{k=1} f(x_0 - k\omega) + \sum^{n-1}_{k=0} f(x_k) + K.\end{align*}\]
When a modern reader encounters a series like this, the question of convergence naturally arises; to us the series only makes sense if the terms \(f(x_0 - k\omega)\) go to zero fast enough for convergence. To Arbogast, Servois and their contemporaries, this series was generally viewed as a formal expression, which can be manipulated by the rules of algebra without concern over convergence. However, Servois was quite concerned with questions of convergence and divergence when they arose in relation to numerical aspects of series. In [Servois 1817, pp.103–106], he even criticized the practice employed by some mathematicians of this period of using divergent series for numerical approximation.
Where did Servois actually use an expression like \(\sum y\)? It turns up in the second part of his equation (9), where we have the expression \(\sum y - \sum v\). But because \(v=f(x_0)\), we have
\[\sum y - \sum v = \left[\sum^{\infty}_{k=1} f(x_0 - k\omega) + \sum^{n-1}_{k=0} f(x_k) + K\right] - \left[\sum^{\infty}_{k=1} f(x_0 - k\omega) + K\right].\]
Formally subtracting the infinite sum from itself, we are left with the finite sum
\[\sum^{n-1}_{k=0} f(x_k).\]
In an analogous fashion, we must understand Servois' use of the notation \(\int y \,dx\) and \(\int v \, da\) to mean the improper integrals
\[\int_{-\infty}^x F(t)\,dt \quad \mbox{and} \quad \int_{-\infty}^a F(t)\,dt,\]
again ignoring questions of convergence. Formal subtraction in the expression for \(Z\)
in Servois' equation (9) gives
\[Z = \int_{-\infty}^x F(t)\,dt - \int_{-\infty}^a F(t)\,dt = \int^x_a F(t) \,dt,\]
an ordinary definite integral. We note that Servois did not use the modern notation for definite integrals. In a different place in the paper, on page [103], he used the notation
\[\int Fx dx \left\{\frac{x}{a}\right\}\]
to stand for the expression that we would write today as
\[\int^x_a F(t)\,dt.\]
The Bernoulli Numbers are named in honor of Jakob Bernoulli (1655–1705), the older brother of Johann Bernoulli (1667–1748) and uncle of Daniel Bernoulli (1700–1782). He studied the formulas for the sums of powers and presented some results about them in his posthumous book Ars conjectandi [Bernoulli 1713]. Bernoulli was not the first person to study these formulas, but he did collect together all the formulas for the powers \(n=1, 2, 3, \ldots, 10\) and noticed a remarkable pattern in the coefficients of these formulas.
Let's use modern notation and define
\[S_m(n) = \sum_{k=1}^n k^m = 1^m + 2^m + 3^m + \cdots + n^m.\]
Students usually learn the first few of these formulas in a calculus class and maybe prove them using mathematical induction:
\[\begin{array}{rll}S_1(n) &= {\displaystyle \frac{n(n+1)}{2}} &= \frac{1}{2} n^2 + \frac{1}{2} n,\\ S_2(n) &= {\displaystyle \frac{n(n+1)(2n+1)}{6}} &= \frac{1}{3} n^3 + \frac{1}{2} n^2 + \frac{1}{6} n,\\ S_3(n) &= {\displaystyle \left(\frac{n(n+1)}{2}\right)^2} & = \frac{1}{4} n^4 + \frac{1}{2} n^3 + \frac{1}{3} n^2.\end{array}\]
There are some patterns here. For example, the first term is always \(\frac{1}{m+1}n^{m+1}\) and the second is always \(\frac{1}{2}n^m\). With the help of the binomial coefficients
\[\binom{n}{k} = \frac{n\cdot(n-1) \cdots (n-k+1)}{k!}\]
Bernoulli eventually found that that there is a sequence of constants\(b_0, b_1, b_2, \ldots\), such that
\[ S_m(n) = \frac{1}{m+1} \sum_{k=0}^m b_k \binom{m+1}{k} n^{m+1-k}.\]
These numbers are called the Bernoulli Numbers. They are usually denoted \(B_k\) instead of \(b_k\), but because Servois used the notation \(B_n\) for something slightly different, we are using \(b_k\) to denote the Bernoulli numbers as they are usually defined today.
It's easy to see from our knowledge of \(S_1(n)\), \(S_2(n)\) and \(S_3(n)\) that
\[b_0 = 1, \quad b_1 = \frac{1}{2}, \quad \mbox{and } b_2 = \frac{1}{6}.\]
What is much less obvious is that \(b_k\) is zero for odd values of \(k\) bigger than one. The even values are
\[b_2 = \frac{1}{6}, \quad b_4 = -\frac{1}{30}, \quad b_6 = \frac{1}{42}, \quad b_8 = -\frac{1}{30} \ldots.\]
There is an alternating sign, but otherwise the pattern followed by the Bernoulli numbers of even index is not a simple one. The Bernoulli numbers can be given by a simple recursion and they can also be found in a variety of infinite series. For example,
\[ \cot x = \frac{1}{x} \sum^{\infty}_{n=0}\frac{(-1)^n b_{2n} (2x)^{2n}}{(2n)!}.\,\,\,\,\,\,\,\,\,\,\mbox{(V)}\]
Because all the Bernoulli numbers of odd index are zero, except the first, some authors use \(B_n\) to denote \(b_{2n}\). When Servois introduced the Bernoulli numbers in formula (6) of his “Memoir on Quadratures” he did something similar; he made the definition \(B_n = |b_{2n}|\). With this notation, the MacLaurin series for cotangent is
\[\cot x = \frac{1}{x} \left[ 1 - \sum^{\infty}_{n=1}\frac{B_n (2x)^{2n}}{(2n)!} \right].\,\,\,\,\,\,\,\,\,\,\,\mbox{(VI)}\]
The Bernoulli Numbers appear in Servois' memoir in a curious way. Because
\[\sum y = (E-1)^{-1} y \quad \mbox{and} \quad Ey = e^d y\]
we must have
\[ (e^d - 1)\sum y = y.\,\,\,\,\,\,\,\,\,\,\,\mbox{(VII)}\]
Servois assumed that the operator \(\sum\) has the form \(Ad^{-1} + Bd^0 + C d^1 + D d^2 + \cdots\) and, using the series expansion of \(e^d\), solved the operator version of equation (VII):
\[\left(d + \frac{d^2}{2!} + \frac{d^3}{3!} + \cdots\right) \left(Ad^{-1} + B + C d + Dd^2 + \cdots \right) = 1.\]
Expanding the left hand side, one term of the first series at a time, we get
\[\begin{array}{rllllllllll} A &+& Bd &+& Cd^2 &+& Dd^3 &+& E d^4 &+& \cdots \\ \\ &+& \frac{A}{2} d &+&\frac{B}{2} d^2 &+& \frac{C}{2} d^3 &+& \frac{D}{2} d^4 &+& \cdots \\ \\ && &+& \frac{A}{6} d^2 &+& \frac{B}{6} d^3 &+& \frac{C}{6} d^4 &+& \cdots \\ \\ && && &+& \frac{A}{24} d^3 &+& \frac{B}{24} d^4 &+& \cdots \\ \\ && && && &+& \cdots &+& \cdots. \end{array}\]
Now sum the columns: the first column sums to \(1\) and the others to \(0\). This gives
\[A=1, \quad B=-\frac{1}{2}, \quad C=\frac{\frac{1}{6}}{2!}, \quad D=0, \quad E=\frac{-\frac{1}{30}}{4!}, \ldots .\,\,\,\,\,\,\,\,\,\,\,\mbox{(VIII)}\]
That is, the numbers \(A, B, C, D, E, \ldots\) are equal to
\[\frac{b_k}{k!}\]
for \(k=0, 1, 2, \ldots\), with the exception that \(B=-b_1\). This is why the Bernoulli numbers appear in Servois' equation (6).
When we apply the integral test to a series, we consider a function \(g\) where \(a_k = g(k)\) and draw conclusions about the convergence of
\[\sum_{k=1}^{\infty} a_k\]
based on the convergence of
\[ \int_1^{\infty} g(x)\, dx.\]
Put another way, the finite integral
\[ \int_1^n g(x)\, dx\]
gives us information about the partial sum
\[ \sum_{k=1}^{n-1} g(k).\]
Of course, the finite integral and the partial sum are not equal, but they are close enough so that each one converges if and only if the other one does. Additionally, the reason they are close together is that the partial sum of the sequence is also a Riemann sum of the integral, with \(\Delta x = 1\) and the left-hand endpoint of each subinterval used in forming the sum.
In the 18th century, both Leonhard Euler (1707–1783) and Colin MacLaurin (1698–1746) worked on the problem of determining the relationship between finite sums and the corresponding integrals. The Euler-MacLaurin Formula says that if \(g\) is \(m\) times continuously differentiable on the interval \([i,j]\), then
\[ \sum^{j-1}_{k=i} g(k) = \int_i^j g(t) \, dt - \frac{1}{2} \left[ g(j) - g(i) \right] + \sum_{k=2}^{m} \frac{b_k}{k!} \left[ g^{(k-1)}(j) - g^{(k-1)}(i) \right] + R, \,\,\,\,\,\,\,\,\,\,\mbox{(IX)}\]
where \(R\) is a remainder term that depends on the function \(g\), the interval \([i,j]\) and the degree \(m\) of the approximation [Kac 2002, p. 96]. We can think of the derivatives \(g^{(k)}(t)\) as providing information concerning how well the partial sum approximates the integral.
In the case of \(i=0\) and \(j=n\), solving equation (IX) for the integral, we have
\[\int_0^n g(t) \, dt = \sum^{n-1}_{k=0} g(k) + \frac{1}{2} g(n) - \frac{1}{2} g(0) - \sum_{k=2}^{m} \frac{b_k}{k!} \left[ g^{(k-1)}(n) - g^{(k-1)}(0) \right] - R. \,\,\,\,\,\,\,\,\,\,\mbox{(X)}\]
Although he never called it by that name, Servois derived the Euler-MacLaurin Formula in his “Memoir on Quadratures” and then used it to resolve questions about the error of approximation in numerical integration formulas. The Euler-MacLaurin Formula first appears in a preliminary form in his equation (6), and then in a form equivalent to formula (X) in his equation (10).
We note that Servois gave the formula as an infinite series involving the derivatives of the function \(g\), rather than with an error term. It was Siméon Denis Poisson (1781–1840) who first derived the error term associated with this formula.
In this section and the next, we provide suggestions for student exercises related to Servois' memoir on numerical integration. We begin with some additional detail about Kramp's notation and algorithm which forms the basis of the exercises at the end of this section.
The Composite Trapezoidal Rule (or just the Trapezoidal Rule) gives an approximate value of a definite integral
\[ \int^b_a f(x) \, dx. \]
To apply the rule, we let \(\Delta x = \frac{b-a}{n}\) and define \(x_i = a + k\Delta x\) for \(0 \le i \le n\). Then
\[ \frac{\Delta x}{2}\left[f(x_0) +2 f(x_1) + 2 f(x_2) + \cdots + 2 f(x_{n-1}) + f(x_n)\right] \,\,\,\,\,\,\,\mbox{(XI)}\]
provides an approximate value of the definite integral, with better and better approximations as the value of \(n\) increases.
In his article [Kramp 1815a], Christian Kramp presented an algorithm that takes a few such approximate values of a definite integral, for different small values of \(n\), and uses an extrapolation procedure to produce a new estimate that is much more accurate than those given by the Trapezoidal Rule.
Let's illustrate Kramp's algorithm by applying it to the definite integral
\[ \int_1^2 \frac{1}{x} dx. \]
Kramp chose this example because we can compare the approximations to the exact value of the integral, which is \(\ln 2 = 0.69314718\cdots\).
The first approximation comes from taking \(n=1\). This gives
\[ S_4 = \frac{1}{2} \left[f(1)+f(2)\right] = \frac{3}{4} = 0.75. \]
The second approximation is given by \(n=2\):
\[ S_2 = \frac{\frac{1}{2}}{2} \left[f(1)+ 2f\left(\frac{3}{2}\right) + f(2)\right] = \frac{17}{24} = 0.708\overline{3}. \]
We take one more approximation, this time with \(n=4\):
\[ S_1 = \frac{\frac{1}{4}}{2} \left[f(1)+ 2f\left(\frac{5}{4}\right)+ 2f\left(\frac{3}{2}\right)+ 2f\left(\frac{7}{4}\right) + f(2)\right] = \frac{1171}{1680} \approx 0.697023810. \]
Let's pause to make some comments about the subscripts, which seem to be wrong. The first observation is that Kramp actually used this notation. In the 1810s some mathematicians used subscripts, which were then a recent invention, but many others did not. Nowhere in his “Memoir on Quadratures” [Servois 1817], nor in any of his other works, did Servois make use of subscripts. None of the great mathematicians of the 17th and 18th centuries had access to this useful notational device.
As for the case in point, Kramp took the smallest value of \(\Delta x = \frac{1}{4}\) as a sort of unit for all of the approximations. Let's call it \(\omega\). Then the first approximation above, with \(n=1\), comes from taking \(\Delta x = 4 \omega\), and so he denoted it \(S_4\). The second approximation has \(\Delta x = 2 \omega\), so it is called \(S_2\). Finally, the last and most accurate estimate uses \(\Delta x = 1 \omega\), so it is called \(S_1\). The larger \(n\) is, the smaller \(\Delta x\) is and the more accurate the approximation is. So intuitively, the value of \(S_0\), which may be thought of as a limit as \(\Delta x\) goes to zero, should be the value of the integral. The problem then is to use \(S_4\), \(S_2\) and \(S_1\) to estimate \(S_0\).
We let \(S(x)\) be the function that gives the values of \(S_n\) (Kramp used \(F\) for this function). Kramp assumed that the function \(S\) was an even function, although he didn't explain why, other than to say that he wanted \(S'(0)=0\). Therefore the MacLaurin series for \(S\) is
\[ S(x) = A + Bx^2 + Cx^4 + D x^6 + \cdots .\]
We then expect that \(S(0)=A\) should be the value of the definite integral, or at least a good estimate to it.
In the current example, we only have three estimates, so we can only get the approximate values of the first three coefficients of the the MacLaurin series. In other words, we'll fit our three data points to the equation
\[ S_n = A + Bn^2 +Cn^4. \]
This gives a \(3\times 3\) linear system:
\[\begin{array}{llrlrll} A & + & B & + & C & = & \frac{1171}{1680} \\ A & + & 4B & + & 16 C & = & \frac{17}{24} \\ A & + & 16B & + & 256C & = & \frac{3}{4} \end{array} \]
This system is small enough to solve by hand. The solution gives the value
\[ A = \frac{4367}{6300} \approx 0.693174603, \]
which is very close to \(\ln 2\). The absolute error of the approximation is \(|A - \ln 2| \approx 0.0000274\), which is less than .004% of the actual value of the integral. This compares with an absolute error of about \(0.00388\) in the estimate \(S_4\), more than 100 times bigger.
Let \(N\) be a natural number, preferably a composite number with many factors. Let \(m = \sigma(N)\), the number of factors of \(N\). We estimate the integral
\[ I = \int^b_a f(x) \, dx \]
by \(m\) different applications of the Trapezoidal Rule. Let \(\omega = \frac{b-a}{N}\).
Now for each factor \(n\) of \(N\), let \(k=\frac{N}{n}\) and let \(S_k\) be the Trapezoidal Rule approximation of the integral \(I\) with \(\Delta x = k \omega\) given by formula (XI).
Suppose that
\[ S(x) = A_0 + A_1 x^2 + \cdots A_m x^{2(m-1)} \]
and fit the values of \(S_k\) to this equation. This gives an \(m \times m\) linear system
\[\begin{array}{llrllrlcl} A_0 & + & A_1 & + & \cdots & + & A_m & = & S_1\\ A_0 & + & 2A_1 & + & \cdots & + & 2^{2(m-1)} A_m & = & S_2\\ &\vdots&&\vdots&&\vdots&&\vdots \\ A_0 & + & kA_1 & + & \cdots & + & k^{2(m-1)} A_m & = & S_k\\ &\vdots&&\vdots&&\vdots&&\vdots\\ A_0 & + & NA_1 & + & \cdots & + & N^{2(m-1)} A_m & = & S_N.\\ \end{array} \]
(Here we implicitly assumed that \(N\) is even. More generally, the second equation corresponds to the smallest proper factor \(k > 1\) of \(N\).)
This linear system is solved and the value of \(A_0\) is used to estimate the definite integral \(I\).
In his paper [Kramp 1815a], Kramp used \(N=12\), so that \(S(x)\) was a polynomial of degree 10
and the estimates \(S_1\), \(S_2\), \(S_3\), \(S_4\), \(S_6\) and \(S_{12}\) were fitted to the curve; see figure 7. This gave him an estimate
\[ \int_1^2 \frac{1}{x} dx \approx 0.69314806, \]
which has an absolute error of about \(8.8 \times 10^{-7}\), a little more than one part in a million of the true value. This is an impressively accurate estimate!
Servois' “Quadratures” paper also contains several gems that could be incorporated as historically-based projects within a numerical analysis course.
The complete English translation of Servois' paper itself can be found here (pdf).
We recall that the equation numbers denoted by Roman numerals below refer back to equations that appeared in earlier sections of this article, in order to distinguish them from the original equation numbers used by Servois in his memoir. There is no correlation between these Roman numerals and his Hindu Arabic numerals.
Download the authors' English translation of Servois' “Memoir on Quadratures” (pdf).
Rob Bradley ([email protected]) is a professor in the Department of Mathematics and Computer Science at Adelphi University. He’s the president of the Euler Society, past Chair of HOM SIGMAA (the History of Mathematics Special Interest Group of the MAA), and past president of the CSHPM (Canadian Society for History and Philosophy of Mathematics). With Ed Sandifer, he wrote Cauchy's Cours d'analyse: An Annotated Translation and edited Leonhard Euler: Life, Work and Legacy. He and Ed collaborated with Sal Petrilli on the definitive English translation of l’Hôpital’s calculus textbook.
Salvatore J. Petrilli, Jr. ([email protected]) is an Associate Professor of Mathematics and Department Chair at Adelphi University. His general research interests include history of mathematics, mathematics education, and applied statistics. However, the majority of his research has been devoted to the life and mathematical contributions of François-Joseph Servois. He is the author of Servois' Solutions peu connues de différens problèmes de géométrie-pratique and, along with Robert E. Bradley and C. Edward Sandifer, published L'Hôpital's Analyse des infiniments petits: An Annotated Translation with Source Material by Johann Bernoulli.