A Simple Derivation of Newton-Cotes Formulas with Realistic Errors

In order to approximate the integral $I(f)=\int_a^b f(x) dx$, where $f$ is a sufficiently smooth function, models for quadrature rules are developed using a given {\it panel} of $n (n\geq 2)$ equally spaced points. These models arise from the undetermined coefficients method, using a Newton's basis for polynomials. Although part of the final product is algebraically equivalent to the well known closed Newton-Cotes rules, the algorithms obtained are not the classical ones. In the basic model the most simple quadrature rule $Q_n$ is adopted (the so-called left rectangle rule) and a correction $\tilde E_n$ is constructed, so that the final rule $S_n=Q_n+\tilde E_n$ is interpolatory. The correction $\tilde E_n$, depending on the divided differences of the data, might be considered a {\em realistic correction} for $Q_n$, in the sense that $\tilde E_n$ should be close to the magnitude of the true error of $Q_n$, having also the correct sign. The analysis of the theoretical error of the rule $S_n$ as well as some classical properties for divided differences suggest the inclusion of one or two new points in the given panel. When $n$ is even it is included one point and two points otherwise. In both cases this approach enables the computation of a {\em realistic error} $\bar E_{S_n}$ for the {\it extended or corrected} rule $S_n$. The respective output $(Q_n,\tilde E_n, S_n, \bar E_{S_n})$ contains reliable information on the quality of the approximations $Q_n$ and $S_n$, provided certain conditions involving ratios for the derivatives of the function $f$ are fulfilled. These simple rules are easily converted into {\it composite} ones. Numerical examples are presented showing that these quadrature rules are useful as a computational alternative to the classical Newton-Cotes formulas.

sufficiently smooth function, models for quadrature rules are developed using a given panel of n (n ≥ 2) equally spaced points. These models arise from the undetermined coefficients method, using a Newton's basis for polynomials. Although part of the final product is algebraically equivalent to the well known closed Newton-Cotes rules, the algorithms obtained are not the classical ones. In the basic model the most simple quadrature rule Qn is adopted (the so-called left rectangle rule) and a correctionẼn is constructed, so that the final rule Sn = Qn +Ẽn is interpolatory. The correctionẼn, depending on the divided differences of the data, might be considered a realistic correction for Qn, in the sense thatẼn should be close to the magnitude of the true error of Qn, having also the correct sign. The analysis of the theoretical error of the rule Sn as well as some classical properties for divided differences suggest the inclusion of one or two new points in the given panel. When n is even it is included one point and two points otherwise. In both cases this approach enables the computation of a realistic errorĒS n for the extended or corrected rule Sn. The respective output (Qn,Ẽn, Sn,ĒS n ) contains reliable information on the quality of the approximations Qn and Sn, provided certain conditions involving ratios for the derivatives of the function f are fulfilled. These simple rules are easily converted into composite ones. Numerical examples are presented showing that these quadrature rules are useful as a computational alternative to the classical Newton-Cotes formulas.

Introduction
The first two quadrature rules taught in any numerical analysis course belong to a group known as closed Newton-Cotes rules. They are used to approximate the integral I(f ) = b a f (x)dx of a sufficiently smooth function f in the finite interval [a, b]. The basic rules are known as trapezoidal rule and the Simpson's rule. The trapezoidal rule is Q(f ) = h/2 (f (a) + f (b)), for which h = b − a, and has the theoretical error while the Simpson's rule is Q(f ) = h/3 (f (a) + 4 f ((a + b)/2) + f (b)), with h = (b − a)/2, and its error is The error formulas (1) and (2) are of existential type. In fact, assuming that f (2) and f (4) are (respectively) continuous, the expressions (1) and (2) say that there exist a point ξ, somewhere in the interval (a, b), for which the respective error has the displayed form. From a computational point of view the utility of these error expressions is rather limited since in general is quite difficult or even impossible to obtain expressions for the derivatives f (2) or f (4) , and consequently bounds for |I(f ) − Q(f )| ∞ . Even in the case one obtains such bounds they generally overestimate the true error of Q(f ). Under mild assumption on the smoothness of the integrand function f , our aim is to determine certain quadrature rules, say R(f ), as well as approximations for its errorẼ(f ), using only the information contained in the table of values or panel arising from the discretization of the problem. The algorithm to be constructed will produce the numerical value R(f ), the correction or estimated errorẼ(f ) as well as the value of the interpolatory rule The true error of S(f ) should be much less than the estimated error of R(f ), that is, for a sufficiently small step h. In such case we say thatẼ(f ) is a realistic correction for R(f ). Unlike the usual approach where one builds a quadrature formula Q(f ) (like the trapezoidal or Simpson's rule) which is supposed to be a reasonable approximation to the exact value of the integral, here we do not care wether the approximation R(f ) is eventually bad, provided that the correctionẼ(f ) has been well modeled. In this case the value S(f ) will be a good approximation to the exact value of the integral I(f ). Besides the values R(f ), E(f ) and S(f ) we are also interested in computing a good estimationĒ S (f ) for the true error of S(f ), in the following sense. If the true error the approximationĒ(f ) is said to be realistic if its decimal form has the same sign as E(f ) and its first digit in the mantissa differs at most one unit, that is, E(f ) = ±0.(d 1 ± 1) · · · × 10 −k (the dots represent any decimal digit). Finally, the algorithm to be used will produce the values (R(f ),Ẽ(f ), S(f ),Ē(f )).
In section 1.

Two models
In this work we consider to be given a panel of n (n ≥ 2) points {(x 1 , f 1 ), (x 2 , f 2 ), . . . , (x n , f n )}, in the interval [a, b], having the nodes x i equally spaced with step where f a sufficiently smooth function in the interval. We consider the following two models: Model A Using only the first node of the panel we construct a quadrature rule Q n (f ) adding a correctionẼ n (f ), so that the corrected or extended rule S n (f ) = Q n (f ) +Ẽ n (f ) is interpolatory for the whole panel, where f [x 1 , x 2 , . . . , x n ] denotes the (n − 1)-th divided difference and a 1 , a 2 , . . ., a n are weights to be determined. Note that Q n (f ) is simply the so-called left rectangle rule, thusẼ n (f ) = n j=2 a j f [x 1 , . . . , x j ] can be seen as a correction to such a rule.

Model B
The rule Q n (f ) uses the first n − 1 points of the panel (therefore is not interpolatory in the whole panel), and it is added a correction termẼ n (f ), so that the corrected or extended rule S n (f ) is interpolatory, Since the interpolating polynomial of the panel is unique, the value computed for S n (f ) using either model is the same and equal to the value one finds if the simple closed Newton-Cotes rule for n equally spaced points has been applied to the data. This means that the extended rules (5) and (6) are both algebraically equivalent to the referred simple Newton-Cotes rules. However, the algorithms associated to each of the models (5) and (6) are not the classical ones for the referred rules. In particular, we can show that that for n odd, the rules Q n (f ) in model B are open Newton-Cotes formulas [3]. Therefore, the extended rule S n in model B can be seen as a bridge between open and closed Newton-Cotes rules.
The method of undetermined coefficients applied to a Newton's basis of polynomials is used in order to obtain S n (f ). The associated system of equations is diagonal, The same method can also be applied applied to get any hybrid model obtained from the models A and B. For instance, an hybrid extended rule using n = 3 points could be written as In this work our study is mainly focused in model A.

Notation and background
Definition 2.1. (Canonical and Newton's basis) Let P k be the vector space of real polynomials of degree less or equal to k (k a nonnegative integer). The set is the canonical basis for P n−1 .
is known as the Newton's basis for P n−1 .
A polynomial interpolatory quadrature rule R n (f ) obtained from a given panel where the coefficients (or weights) c j (1 ≤ j ≤ n) can be computed assuming the quadrature rule is exact for any polynomial q of degree less or equal to n− 1, that is, deg(R n (f )) = n − 1, according to the following definition The degree of the quadrature rule is denoted by deg(R). When deg(R) = n − 1 the rule is called interpolatory.
In particular for a n-point panel the interpolating polynomial p satisfies where p can be written in Newton's form (see for instance [6], p. 23, or any standard text in numerical analysis) and the remainder-term is where . . , n. Therefore from (10) we obtain where E Rn (f ) denotes the true error of the rule R n (f ). Thus, The expression (14) suggests that the application of the undetermined coefficient method using the Newton's basis for polynomials should be rewarding since the successive divided differences are trivial for such a basis. In particular, the weights for the extended rule in model A are trivially computed.
Proposition 2.1. The weights for the rule S n (f ) in model A are Proof. The divided differences do not depend on a particular node but on the distance between nodes. Thus for any given n-point panel of constant step h, we can assume without loss of generality that The undetermined coefficients method applied to the Newton's basis < w 0 (t), w 1 (t), . . . , w n−1 (t) > leads to the n conditions a i = I(w i−1 ) or, equivalently, to a diagonal system of linear equations whose matrix is the identity. The equalities in (15) can also be obtained directly from (14).
Theoretical expressions for the error E Rn (f ) in (13) can be obtained either via the mean value theorem for integrals or by considering the so-called Peano kernel ( [1], p. 285, [2], p. 176). However, we will use the method of undetermined coefficients whenever theoretical expressions for the errors E Qn and E Sn are needed.
For sufficiently smooth functions f , the fundamental relationship between divided differences on a given panel and the derivatives of f is given by the following well known result, Applying (16) to the canonical or Newton's basis, we get (17) By construction, the rules S n in models A and B are at least of degree n − 1 of precision according to Definition 2.2.
In this work the undetermined coefficients method enables us to obtain both the weights and theoretical error formulas. This apparently contradicts the following assertion due to Walter Gautschi ([2], p. 176): "The method of undetermined coefficients, in contrast, generates only the coefficients in the approximation and gives no clue as to the approximation error".
Note that by Definition 2.2 the theoretical error (1) says that deg(Q) = 1 for the trapezoidal rule, and from (2) one concludes that deg(Q) = 3 for the Simpson's rule. This suggests the following assumption.
Assumption 2.1. Let be given a n-point (n ≥ 1) panel with constant step h > 0, a sufficiently smooth function f defined on the interval [a, b], and a quadrature rule R(f ) (interpolatory or not) of degree m, there exists a constant K h = 0 (depending on a certain power of h) and a point ξ, such that , and m is the least integer for which (18) holds.
The expression (18) is crucial in order to deduce formulas for the theoretical error of the rules in model A or B.
where ω m+1 (x) is either the element φ m+1 (x) of the canonical basis, or the element w m+1 (x) of the Newton's basis, or any polynomial of degree m + 1 taken from any basis for polynomials used to apply the undetermined coefficients method.

Proof. For any nonnegative integer
As m is the least integer for which the righthand side of (18) is non zero, and deg(R) = m, one has from which it follows (19). As for a fixed basis the interpolating polynomial is unique, it follows that the error for the corresponding interpolatory rule is unique as well.
In Proposition 2.4 we show that deg(S n (f )) depends on the parity of n. Thus, we recover a well known result about the precision of the Newton-Cotes rules, since S n (f ) is algebraically equivalent to a closed N ewton − Cotes formula with n nodes. Let us first prove the following lemma.
Lemma 2.1. Consider the Newton's polynomials Let n ≥ 1 be an integer and I [−1] , I [0] and I [−2] the following integrals: As the integrand is an odd function in I, we have I [−1] (w n ) = 0. For n even, we get where the integrand is an even function, thus The proof is analogous so it is ommited.
The degree of precision for the rules in models A and B, and the respective true errors can be easily obtained using the undetermined coefficients method, the Lemma 2.1 and Proposition 2.3. The next propositions ( 2.4, 2.5 and 2.6) establish the theoretical errors and degree of precision of these rules. In particular, in Proposition 2.4 we recover a classical result on the theoretical error for the rule S n (f ) -see for instance [4], p. 313.
(ii) If n is even and f ∈ C n (J), then Proof. We can assume without loss of generality that the panel is {(0, f 1 ), (h, f 2 ), . . ., ((n − 1)h, f n )} (just translate the point x 1 ) . For n even or odd, by construction of the interpolatory rule S n (f ), we have deg(S n ) ≥ n − 1 in model A or B. Taking the Newton's polynomial w n (t) = t(t − h) · · · (t − (n − 1)h)), whose zeros are 0, h, . . . , (n − 1) h, we get for the divided differences in (5), where f has been substituted by w n in (5). Thus, S n (w n ) = 0 and S j (w n ) = 0 for any j ≥ n + 1.
To the point θ it corresponds a point ξ in the interval J, and so (23) holds. w j−1 (t)dt, j = 1, 2, . . . , n. The entries in column 2 (heading a 1 ) should be multiplied by h; the entries in column 3 (heading a 2 ) multiplied by h 2 , and so on.

Realistic errors for model A
Note that, by construction, the weights in model A are positive for any n ≥ 2. Therefore, the respective extended rule S n (f ) does not suffers from the inconvenient observed in the traditional form for Newton-Cotes rules where, for n large (n ≥ 9) the weights are of mixed sign leading eventually to losses of significance by cancellation. The next Proposition 3.1 shows that a reliable computation of realistic errors for the rule S n (f ), for n ≥ 3, depends on the behavior of a certain function g(x, h) involving certain quotients between derivatives of higher order of f and its first derivative. Fortunately, in the applications, only a crude information on the function g(x, h) is needed, and in practice it will be sufficient to plot g(x, h) for some different values of the step h, as it is illustrated in the numerical examples given in this Section (for some simple rules) and in Section 4 (for some composite rules) .
where a i = I(w i−1 ), i = 1, 2, . . . , n. Denote by g(x, h) (or g(x) when h is fixed) the function Assuming that and then, for a sufficiently small step h > 0, the correctionẼ n (f ) is realistic for Q n (f ). Furthermore, the true error of S n (f ) can be estimated by the following realistic errors: (a) For n odd:Ē Proof. (a) By Proposition 2.4 (i), we have and from Proposition 2.2 the correctionẼ n (f ) can be written as x 3 ), . . ., θ n−2 ∈ (x 1 , x 2 , . . . , x n ). Therefore, Thus, using the hypothesis in (25) we obtain , that is, where c is a constant not depending on h. Thus, by (26) we get Therefore, for h sufficiently small, Finally, by Proposition 2.2 and the continuity of the function f (n+1) (x), we know that So, from (29) we obtain (27). (b) The proof is analogous so it is omitted.
The next proposition shows that Proposition 3.1 for the case n = 2 leads to the rule S 2 (f ) which is algebraically equivalent to the trapezoidal rule, and when n = 3 the rule S 3 (f ) is algebraically equivalent to the Simpson's rule.
thenẼ 2 (f ) is a realistic error for Q 2 (f ), for h = (b − a) sufficiently small. A realistic approximation for the true error of S 2 (f ) is Proof. By Proposition 2.2 there exists a point θ ∈ (x 1 , , soẼ(f ) = h 2 /2 f ′ (θ). Using Proposition 2.4, we know that As by hypothesis f ′ (θ) = 0, we get and so Therefore, for a sufficiently small h the correctionẼ 2 (f ) is realistic for Q 2 (f ). Since we obtain (32) from (33). Let x is not differentiable at x = 0. However f ′ (x) = 0 for x > 0, and the result (32) still holds. The numerical results (for 6 digits of precision) are: The true error for S 2 (f ) is
where x 1 = a, x 2 = (a + b)/2 and and then, for a sufficiently small h = (b − a)/2,Ẽ 3 (f ) is a realistic correction for Q 3 (f ). A realistic approximation to the true error of S 3 (f ) is Proof. As I(w 4 ) =  Let 0, the condition (38) holds with x 1 = 0 and x 3 = 2h. Consider As lim x→0 g(x, h) = 1 and for any 0 < x ≤ 1 we have g(x, h) > h, for 0 < h ≤ 1 (see Figure 1), thus the condition (39) is satisfied. Notice that g(x, h) gets closer to the value 1 as h decreases. Therefore one can assure that realistic estimates (40) can be computed to approximate the true error of S 3 (f ) for a step h ≤ 1.
In Table 3 is displayed the estimated errorsĒ S3 (f ) and the true error for S 3 (f ), respectively for h ∈ {1/2, 1/4, 1/8, 1/16}. As expected, the computed values forĒ S3 (f ) have the correct sign and closely agree with the true error. For h = 1/2, we have   From Table 1 we obtain the following expression for the rule S 5 (f ), where the coefficients a i are computed using (15). It can be observed in the plot in the Figure 2 that for h ∈ {1/8, 1/16, 1/32, 1/64} the condition g(x, h) > h is satisfied. Therefore, since n is odd, one concludes from Proposition 3.1 (a) that the following realistic estimation for the true error of S 5 (f ) is, Table 4 are displayed the computed realistic errorsĒ S5 (f ) for the steps h referred above. For instance for h = 1/8, we have hĒ

Composite rules
The rules whose weights have been given in Table 1 are here applied in order to obtain the so-called composite rules. The algorithm described hereafter for composite rules is illustrated by several examples presented in paragraph 4.1.
Since the best rules S n are the ones for which n is even (when deg(S n ) = n holds) the examples refer to S 3 , S 5 , S 7 and S 9 . Whenever the conditions of Proposition 3.1 for obtaining realistic errors are satisfied, these rules enable the computation of high precision approximations to the integral I(f ), as well as good approximations to the true error. This justifies the name realistic error adopted in this work. . , x N } into subsets of n points each, and for an offset of n − 1 points, we get i panels each one containing n successive nodes. To each panel we apply in succession the rule Q n and compute the respective realistic correctionẼ n as well as the estimated realistic errorĒ Sn for the rule S n (f ). For the output we compute the sum of the partial results obtained for each panel as described in (43): According to Proposition 2.4, for composite rules with n points by panel and step h > 0, in the favorable cases (those satisfying the hypotheses behind the theory) one can expect to be able to compute approximations S having a realistic error. Analytic proofs for realistic errors in composite rules for both models A and B will be treated in a forthcoming work [3].

Numerical examples for composite rules
Once computed realistic errors within each panel for a composite rule, we can expect the errorĒ in (43) to be also realistic. This happens in all the numerical examples worked below. Let In the interval J = [10 5 , 2 × 10 5 ], the function f (x) = 1/ ln(x) belongs to the class C ∞ (J). Since f ′ (x) = −(x ln 2 (x)) −1 = 0, for all x ∈ J, and for k ≥ 2 the quotients f (k) (x)/f ′ (x), with x ∈ J, are close to γ(x) = 0 and tends to the zero function γ(x) as k increases. Therefore, the lefthand side in the inequality (26) is very close to 1. That is, for n ≥ 2 the function is such that g(x, h) ≃ 1, for h sufficiently small. So Proposition 3.1 holds and one obtains realistic errors for the rules S n (f ). The behavior of the function g(x, h) is illustrated for the case n = 3 in Figure  3, for h ∈ {5, 5/2, 5/3}. Note that in this example g(x) ≃ 1 while the chosen steps h are greater than 1. However, realistic errors are still obtained. In Table 5 the realistic error is compared with the true error, respectively for the rules with n odd from 3 to 9 points (the step h is as tabulated).  For n = 7, the computed approximation for the integral I is S = 8406.24312084620270862164604369467068, where all the digits are correct. The simple rule S 7 (f ) is defined (see Table 1 In Appendix 4.1 a Mathematica code for the composite rule S, for n = 7, is given. The respective procedure is called Q7A and the code includes comments explaining the respective algorithm. Of course we could have adopted a more efficient programming style, but our goal here is simply to illustrate the algorithm described above for the composite rules. (* realistic error for q *) s = q + E7; (* s is a realistic approximation to the integral *)