Behavior of the Numerical Integration Error

In this work, we consider different numerical methods for the approximation of definite integrals. The three basic methods used here are the Midpoint, the Trapezoidal, and Simpson’s rules. We trace the behavior of the error when we refine the mesh and show that Richardson’s extrapolation improves the rate of convergence of the basic methods when the integrands are sufficiently differentiable many times. However, Richardson’s extrapolation does not work when we approximate improper integrals or even proper integrals from functions without smooth derivatives. In order to save computational resources, we construct an adaptive recursive procedure. We also show that there is a lower limit to the error during computations with floating point arithmetic.


Introduction
(see for example [1] and [2]). In those cases the Fundamental Theorem of Calculus cannot be applied and the values of the integral have to be approximated numerically.
It is clear that when using numerical calculations, due to the round-off error, we usually cannot obtain the exact value I of the integral. So, we say that the numerical value of the integral I is obtained when we obtain I  approximate value of I , such that I I ε − <  , where 0 ε > is the acceptable error. Even in some cases when the Fundamental Theorem of Calculus is applicable, we still have to work with approximate value of the integral. For example 1 0 e e 1 1.7172 Since e is a transcendent number we cannot represent it as a finite decimal number. So, we have to use an approximation of the integral. Usually, the approximate value of the integral is calculated using the so-called formulas for numerical integration or quadrature formulas. The general idea is to approximate the integral as a sum A P x = = ∑ ∫ (6) In this work we consider so called Newton-Cotes formulas defined on a fixed set of nodes (see [3] for details). In this case, in order to obtain the coefficients i A , one has to solve the system of linear equations

The Error Term
The error term, ( ) If the function ( ) f x has a continuous m derivative, the error term can be represented in the form (see for details [3]

Midpoint Rule
The simplest formula, using only one node, is the Midpoint Rule. Suppose the node is the midpoint a b . By solving the Equation (7), we obtain the formula

Composite Midpoint Rule
In order to reduce the error term we partition the Figure 1). Suppose the length of each subinterval is the same, is given as where 24 b a c − = .
Suppose we have performed two calculations with different numbers of subintervals, N and 2N . Then and Therefore, if we double the number of subintervals, the error term is decreased four times. This is a practical way of confirming the power of h in the error term. In practice [5], the ratio is not exactly four, because f ′′ is usually not constant.

Numerical Experiments
We prepared a SAGE function MPloop for our numerical experiments with the Midpoint Rule. The code is given here: The variables a and b are the lower and upper limits of the integral, and the variable n is the number of subintervals.
The results obtained with the function MPloop for the integral ( ) is going to zero. One application of the Midpoint Rule can be found in [6].

Trapezoidal Rule
The formula based on two nodes, 0 x a = and 1 x b = , and exact for all second order polynomials, can be calculated via the system of two Equations (7). The resulting formula is called Trapezoidal Rule, because it produces exactly the area of the trapezoid with vertexes ( ) ,0 If the function ( ) f x has a continuous second derivative, the error term of the Trapezoidal Rule can be represented in the form (see [3])

Composite Trapezoidal Rule
The composite formula for the Trapezoidal Rule is   The nodes k x are defined in subsection 2.1 as given in Figure 1. The error term is (see [5] for details) As with the Midpoint Rule, we obtain

Numerical Experiments
We construct a SAGE function TPloop(a,b,n) for the composite Trapezoidal Rule: The variables a and b are the lower and upper limits of the integral, and the variable n is the number of subintervals.
The results obtained with function TRloop for the integral (18) are shown in Table 2. The last column shows the ratio of the errors Our numerical calculations confirm that if we double the number of subintervals, the error term is decreased four times. The error distribution for the same calculation is shown graphically in Figure 3.

Simpson's Rule
Suppose we want to use three nodes 0 x a = , 1 2 For this case the quadrature Formula (5) has three free (unknown) coefficients  and we can make the formula exact for second order polynomials. In other words, the formula must be exact for all polynomials of zero degree and particularly for ( ) 1 The formula must be exact for all first order polynomials and particularly for ( ) The formula must be exact for all second order polynomials and particularly for ( ) 2 P x x = . Therefore, 2 3 The last three equations form the system (7) for this case. The solution to the system is 0 , and the Simpson's Rule is given by The error term, under the condition that the fourth derivative , .

Composite Simpson's Formula
One way to construct the composite Simpson's formula is to divide the interval [ ] The error term in this case is , , As with the two previous rules, we obtain So, if the number of subintervals is multiplied by 2, the numerical error decreases by a factor of 1 16 .

Numerical Experiments
We prepared a SAGE function Simpsons for our numerical experiments with Simpson's Rule. The code is given here: The variables a and b are the lower and upper limits of the integral, and the variable n is the number of subintervals.
The results obtained with function Simpsons for the integral

Error Estimation from Numerical Results
How do we determine the number of subintervals such that the computed result meets some prescribed accuracy?
Next we will consider the Runge principle, Richardson's extrapolation, and the power of h estimation in the error term to the Midpoint, Trapezoidal, and Simpson's rules.

Midpoint Rule
A simple SAGE implementation of the Runge principle for estimating the error during the numerical calcula-tions is given here: The variables a and b in the definition of the function MPerror contain the lower and upper limits of the integral. The variable eps is the satisfactory error. The function returns: the first variable, Inew, is the approximation of the integral with the Midpoint Rule; the second variable, Iextr, is the value of Richardson's extrapolation for the integral using Iold and Inew; the third variable n gives the number of subintervals for the final approximation.
To verify the theoretical results for the Midpoint method, we calculate the integral

Trapezoidal Rule
When the number of subintervals N increases, the error in the composite Trapezoidal Rule reduces. On the other hand, when N increases, the number of the function values to be calculated also increases. We need to find a balance between the accuracy requirements and the cost (in operations or computing time) of the algorithm. This balance can be reached with the following algorithm: We calculate the approximation of the integral 0 I by applying the Trapezoidal Rule over the whole interval [ ] , a b . Then we calculate the approximation 1 I by dividing the interval [ ] , a b into two subintervals, and in each we apply the Trapezoidal Rule. In this partition the new node is only one in the middle of this interval. The other two nodes are involved in the calculation of 0 I and we do not need to calculate the value of the function on those nodes again. We calculate 2 I by introducing four subintervals. The new nodes are only two. So, we find 3 4 , , , ,  A simple SAGE implementation of the Runge principle for estimating the error during the numerical calculations with Trapezoidal Rule is given here: The variables a and b in the definition of the function TRerror contain the lower and upper limits of the integral. The variable eps is the satisfactory error. The function returns: the first variable, Inew, is the approximation of the integral with Trapezoidal Rule; the second variable, Iextr, is the value of Richardson's extrapolation for the integral using Iold and Inew; the third variable n gives the number of subintervals for the final approximation.
To verify the theoretical results for the Trapezoidal method, we calculate the integral (36) with a different number of subintervals and approximate the value of the integral with Richardson's extrapolation. The results are given in Table 5 The fact that Richardson's extrapolation applied to the Trapezoidal Rule gives fourth order of approximation has a simple explanation. Since,

Simpson's Rule
To verify the theoretical results for Simpson's method, we calculate the integral (36) with a different number of subintervals, and approximate the value of the integral with Richardson's extrapolation. The results are given in Table 6

Adaptive Procedures
According to [5]: "An adaptive quadrature routine is a numerical quadrature algorithm which uses one or two basic quadrature rules and which automatically determines the subinterval size so that the computed result meets some prescribed accuracy requirement.'' For our adaptive function we chose the Midpoint Rule, because it involves just one node to approximate the integral in the given subinterval. The SAGE function needs to be changed slightly if one wishes to use a different quadrature formula. The SAGE code for the function MPrecursion(a,b,eps) is given here: The parameters of the function MPrecursion are a and b-the lower and upper limits of the integral, and eps the accuracy. The function divides the interval into two subintervals, and, in each subinterval, approximates the integral twice using the Midpoint Rule for the full half interval and using the composite Midpoint Rule, which divides the half interval into two subintervals. If the error, according to Runge's principle, is less than eps, the value calculated with the composite rule is accepted as an approximation of the integral for this subinterval. If the error is greater than the prescribed accuracy, the function calls itself with a half-length interval and halved accuracy requirements. Therefore, we have one recursive algorithm which uses different mesh sizes for different parts of the interval. If the integrand is smooth and slowly varying, a relatively small number of subintervals is used. In the parts of the interval where the integrand is changing rapidly, the mesh is relatively fine.
An example of the behavior of the function MPrecursion is given in Figure 5.

Can We Always Trust the Runge Principle?
As we have shown in the previous section, the Runge principle for estimating the error during the numerical calculations works well under the condition that the respective derivative of the integrand exists and this derivative is a continuous function. If the derivative in the error term does not exist, the Runge principle fails, and Richardson's extrapolation cannot be applied. Consider the following examples.

Improper Integrals
The formal definition of improper integral is: Examples: Calculating improper integrals numerically is challenging. Usually, we need some a priori information about the integrand, because the techniques we apply depend on the properties of this function. Two main problems arise: 1) the integrand is not defined at one (or both) limits of the integral; 2) the integrand goes to infinity at one (or both) limits of the integral; Because of 1) Trapezoidal and Simpson's rules cannot be used. The Midpoint Rule can be used, since we do not need the value of the function at the endpoints a and b .
The results from calculations of the first improper integral from (40) with the Midpoint Rule are given in Table 7. The ratio of the errors in this case is about 1.414. Compared to the ratio of 4 for the "good'' functions, we see that the convergence is very slow. We will have a similar convergence for a "good'' integrand when we approximate the integral with a formula with error term ( )

Singularity in the Derivative
The slow convergence of the improper integral is not a surprise, since the integrand goes to infinity at one of the endpoints of the interval. Consider the integral  The function y x = is a continuous and bounded function, but the second derivative of this function goes to infinity when 0 x → . So, we can expect problems with practical convergence of the numerical calculation of the integral. The results from the numerical calculations of (41) with the Midpoint Rule are given in Table 8. The ratio of the errors in this case is about 2.8. Comparing this ratio to 4, for the "good'' functions, we see that the convergence is slow. We will have a similar convergence for a "good'' integrand when we approximate the integral with the formula with error term ( ) 3 3 , 1 R f N c h c N = = .

Behavior of the Error in Practice
We have to emphasize that in practical calculations the error in the final result depends on the accumulation of rounding errors. With the increasing number of subintervals the theoretical error decreases, but the number of arithmetic operations, each of which contains an error in the order of machine epsilon, increases. Therefore it is reasonable to use as many subintervals as we need to get the theoretical error term to exceed the machine epsilon. The total error starts to increase when the number of subintervals increases beyond a certain unknown number dependent on the integrand, the programming language, and the computer.
To illustrate this effect with our SAGE codes we use the class RealField with 32 bits. The effect when the Sage 64-bit floating-point system is used is similar but the computing time increases. The numerical results for three integrals, are shown in Figure 6. When the number of subintervals is relatively small, the error decreases as the number of subintervals increases. When the number of subintervals becomes greater than a certain number, the error starts to increase. The number after which the error starts to increase depends on the integrand.