Integration plays a central role in the mathematical models that have been developed to price options. Integration usually involves finding the antiderivative of a function and applying the Fundamental Theorem of Calculus, according to which the antiderivative is evaluated at the endpoints of the integral only. This well-known analytical procedure yields the area under the curve defined by the function, in between the two endpoints.
Unfortunately, while many option pricing models require evaluating the integral of a function to find the area under the curve, it is often impossible to find an antiderivative for the function. One example is the option pricing model of Heston and Nandi (2000) reviewed in Chapter 6, for which they find a semiclosed form solution for an option price, on an underlying asset that has volatility defined by a generalized autoregressive conditional heteroskedasticity (GARCH) process. The solution is in semiclosed form because finding the option price involves a complex integral that must be evaluated by numerical integration. We will need also to apply numerical integration to calculate model-free volatility, skewness, and kurtosis, in later chapters.
This chapter introduces numerical integration, presents algorithms that are used to perform numerical integration, and illustrates the VBA code to implement these algorithms. There are two main families of algorithms that have been developed for numerical integration: Newton-Cotes formulas, also called quadratures, and Gaussian quadratures. The main difference between these algorithms is the way in which the function being integrated is divided into subintervals along the horizontal axis. Newton-Cotes formulas divide the axis into subintervals at each integration point, or abscissa. Lagrange polynomials of varying complexity are used to approximate the areas under the function between successive pairs of abscissas. These areas are then aggregated to produce an approximation to the original integral. The subintervals are usually chosen to be of equal width.
Gaussian quadratures use a much more complex algorithm to approximate the integral. The chief difference is that the abscissas are determined not by the analyst, but by the algorithm itself, as the roots of approximating polynomials. Gaussian quadratures involve the sum of the product of the function evaluated at the optimal abscissas and a weighting function. The fundamental theorem of Gaussian quadratures defines the abscissas as the roots of the orthogonal polynomial for the same interval and weighting function. Although more complicated than Newton-Coates formulas, Gaussian quadratures are more accurate because they fit polynomials up to degree 2n − 1 exactly, where n is the number of abscissas. With today’s high computer speeds, however, a very large number of points can be chosen, and Newton-Cotes approximations can be computed very quickly. In light of this, for many financial analysts the increased accuracy of the Gaussian quadratures does not justify the extra programming effort required to implement them. It is therefore not surprising that the more straightforward Newton-Cotes formulas are the most widely used of numerical integration algorithms. Consequently, in this chapter we focus mostly on Newton-Cotes formulas. In later chapters, however, the superiority of Gaussian quadratures over Newton-Cotes formulas will become evident.
In this section, seven popular Newton-Coates numerical integration formulas are introduced: the left-, right-, and midpoint rules, the trapezoidal rule, Simpson’s rule and Simpson’s three-eighths rule, and Boole’s rule. We present the VBA code for implementing each rule, and compare the accuracy of each using a simple example. Finally, we assess how the accuracy of each approximation increases as the number of integration points increases.
Finding the integral of a function is equivalent to finding the area under the curve traced by that function. One way to approximate that area is by dividing the area into rectangles and adding the rectangles together. Suppose we wish to approximate the area under the curve traced by the function f(x) = x3 + 10. This is, of course, a simple function for which an antiderivate exists. It can be integrated exactly, but it is used throughout this chapter to illustrate the methodology of each rule. To approximate the area under the curve within the interval [a, b] on the horizontal axis (x-axis), n + 1 equally spaced integration points, or abscissas, are used to divide the interval into a set of n subintervals of equal width:
where x0 ≡ a and xn ≡ b are the endpoints. A rectangle is constructed in each subinterval and its area is obtained. Aggregating the areas of the rectangles produces the approximation. The left-, right-, and midpoint rules differ only by how the height of the each rectangle is chosen.
According to the right-point rule, the height of each rectangle is defined as the value of the function f(x) at the rightmost endpoint of each subinterval. Figure 2.1 illustrates the right-point rule using our function as an example, with a = −2 and b = 2. This rule can be expressed mathematically as the approximation
FIGURE 2.1 Right-Point Rule Numerical Integration
where
Δx = (b − a)/n = the width of each subinterval
xi = the rightmost point of each subinterval
a, b = endpoints of the interval over which the integral is evaluated
f (xi) = the value of the function evaluated at xi
n = the number of subintervals.
Note that each abscissas xi can be expressed in terms of the leftmost endpoint a and the subinterval width Δx as xi = a + iΔx, for i = 0, 1, 2, …, n. This convenient representation will be exploited later in this chapter, when we implement numerical integration with VBA.
The left-point rule sets the height of each rectangle to the leftmost value of f(x) in each subinterval, so that the approximation is
Figure 2.2 illustrates our example using the left-point rule approximation rule. Figures 2.1 and 2.2 demonstrate the shortcomings of the left- and right-point rules. Although these rules each produce a numerical approximation to the integral, a substantial amount of information is lost between the n + 1 points defined by the subintervals. It easy to see that when the function is monotone nondecreasing (as in our example), the right-point rule produces an approximation that is too large, while the left-point rule produces one that is too small. The reverse is true for a nonincreasing monotone function.
FIGURE 2.2 Left-Point Rule Numerical Integration
There are two ways to incorporate information ignored by the left- and right-point rules. The first method requires new function outputs, namely the y-values defined at the midpoint of each subinterval. This is aptly named the midpoint rule. The midpoint rule is identical to the left- and right-point rules, except that it uses the value of f(x) obtained at the midpoint of each subinterval for the height of the rectangles.
Mathematically, the midpoint rule can be expressed as
where
As illustrated by Figure 2.3, the increase in accuracy brought on by choosing the midpoint rule over either the right-point or left-point rule is substantial. Nevertheless, these three methods are relatively crude compared to the methods that follow.
FIGURE 2.3 Midpoint Rule Numerical Integration
The trapezoidal rule does not require additional information other than the value of f(x) at each abscissa. Implementing this rule is no more complicated than implementing the left-, right-, or midpoint rules. The interval [a, b] is again divided into a set of n subintervals of equal widthΔx. The goal of this approach, however, is to approximate the rate at which the function f(x) is increasing or decreasing within each subinterval. This is achieved by forming a line segment to join f(x) at the endpoints of each subinterval, which produces a trapezoid instead of a rectangle. The area of the trapezoid defined by the subinterval (xi−1, xi) is
Aggregating this quantity for all subintervals produces the trapezoidal rule, which can be expressed mathematically as
Because trapezoids are used instead of rectangles, the accuracy of the trapezoidal approximation is much higher than that achieved with the midpoint rule. Figure 2.4 illustrates our example with the trapezoidal rule approximation.
FIGURE 2.4 Trapezoidal Rule Numerical Integration
Simpson’s rule is a refinement over the trapezoidal rule because it uses parabolas instead of line segments to join points of f(x). Simpson’s rule, however, requires that n, the number of subintervals, be an even number. This is needed because the rule jointly approximates each pair of two consecutive areas.
Simpson’s rule is straightforward to derive. For simplicity, assume that we wish to approximate the area under the function f(x) in the interval [−a, a], using n = 2 subintervals. We know that a parabola can be fitted to any three points in the xy-plane, so we choose the points f(−a), f(0), and f(a). Suppose that the parabola joining these points is of the form
Because the antiderivative of f(x) can be obtained explicitly, the exact area under the parabola in the interval (−a, a) is
The goal is to express this known area in terms of the three points f(−a), f(0), and f(−a) since these points constitute the only known information for implementing the approximation. Since the parabola passes through the points (−a, f(−a)), (0, f(0)), and (a, f(a)), we have the following three equations:
Multiplying the second equation by four and adding the three equations, we obtain
Comparing this last equation to the integral derived above, we see that the area under the parabola between −a and a is equal to 
. This area will not change if the parabola is shifted along the x-axis, so this relation holds for intervals other than those centered about zero. Considering this fact, we can use this last equation for any three consecutive points on a curve. Hence, we obtain Simpson’s rule for numerical integration:
where, as before, x0 = a, x2 = b, and Δx = (b − a)/2 is the width of each subinterval.
If the interval [a, b] is very narrow, then Simpson’s rule will produce an accurate approximation of the integral. In most applications, however, the interval is wide. The obvious remedy is to break up the interval into a series of small subintervals, and apply Simpson’s rule to each subinterval. We can use the rule to approximate the area for any even number of subintervals, by overlapping the three sums for each pair of subintervals. This produces an algorithm sometimes known as the composite Simpson’s rule, which can be expressed as the approximation
where again x0 = a, x2 = b, Δx = (b − a)/n is the width of each subinterval, and n is the number of subintervals. In order for this method to work, n must be an even number because each application of the rule requires two subintervals.
Simpson’s rule can be refined to yield an approximation method known as Simpson’s three-eighths rule. This rule uses a cubic polynomial instead of a parabola to approximate f(x). Hence, to construct the approximation under this rule, four points instead of three are needed in the interval [a, b]. The approximation can be derived in a manner similar to that which produces Simpson’s rule. Simpson’s three-eighths rule can be written as the approximation
As before, we can we can break up the interval [a, b] into a set of n subintervals, where now we require n to be divisible by three, apply Simpson’s three-eighths rule in each subinterval, and aggregate the resulting sums. This is known as the composite Simpson’s three-eighths rule.
Now that we have Simpson’s rule for integrating a function with an even number of integration subintervals, and Simpson’s three-eighths rule for a number divisible by three, we can perform a numerical integration for any number of segments by combining the two methods. For example, suppose we have a function that we can evaluate at six evenly spaced points, thus creating five integration subintervals. We could use Simpson’s three-eighths rule for the first three segments and Simpson’s rule for the remaining two segments. The approximation would therefore be
The last Newton-Cotes method we introduce is Boole’s rule. This approximation requires four subintervals of integration and can be expressed where
As with Simpson’s rule and Simpson’s three-eighths rule, a composite version of this rule can be applied to a set of n subintervals within [a, b]. Since four points are needed for each subinterval, the number of subintervals used in the approximation must be divisible by four.
Finally, the trapezoidal rule, Simpson’s rule, Simpson’s three-eighths rule, and Boole’s rule all use special cases of Lagrange interpolating polynomials. In their most general form, each rule can be written as
where 
and Ln, i (x) is the Lagrange interpolating polynomial of degree n.
In the four approximation rules presented so far in this chapter, equally spaced abscissas are specified, or equivalently, subintervals of equal width. While this simplifies the algebraic formulas for the rules, there is no reason to impose equal width in any of these approximations. The VBA functions presented later in this chapter allow the abscissas to be arbitrarily spaced, and not evenly spaced only. This permits subintervals of any width to be used in the approximation. The left-, right-, and midpoint rules, and the trapezoidal rule, allow for intervals of arbitrary width. For example, for arbitrarily spaced points (x0, x1, …, xn) the right-point approximation rule (2.1) becomes
As before, the endpoints of the integral are x0 = a and xn = b.
Adaptive quadrature methods are those that specify small subinterval widths in regions of high functional variation, and larger widths in regions where the function is flatter. It is straightforward to implement adaptive quadrature in the approximation rules. For example, suppose the approximation is over the integration interval [a, b]. The interval is broken up into two subintervals [a, (a + b)/2] and [(a + b)/2, b], and two approximations are obtained. One approximation uses the entire interval, and the other uses the sum of the approximations from each subinterval. If the difference between these approximations is larger than some tolerance level, the smaller subintervals are used; otherwise, [a, b] is used.
When an integration rule uses all the abscissas x0, x1, …, xn in its formulation, the rule is said to be a closed rule. When the endpoints x0 = a and xn = b, are excluded, however, the rule is an open rule. Among the rules introduced in this chapter, the midpoint rule is an open rule, while the other rules are all closed. Closed rules can be fully open, so that both endpoints are excluded from the approximation, or semiopen, so that only one endpoint is excluded.
Open rules are needed to approximate improper integrals for which the function cannot be evaluated at the endpoints of the integration interval. For example, the improper integral
can be solved analytically, and has the value 
. Using Simpson’s rule on this integral, however, would yield an infinite value for the approximation because the function explodes at the left endpoint a = −1. Hence, in this example a semiopen rule would be required.
There are at least two simple ways to transform closed rules into semiopen or open rules. The easiest way is to replace the endpoints of the integration integral by values that are shifted by a small amount, ε, so that a and b are replaced by the points a + ε and b − ε, respectively. Any closed rule can now be applied to this new set of integration points. If ε is too small, however, a large number of abscissas may be required for the approximation to work properly since the function may not be well-behaved at points too close to the endpoints. In the above example, it easy to see that 
could take on very large values if ε is too small. The only way to mitigate this inflated value in the approximation is to define a very small subinterval width.
Open rules can also be created by choosing 
and 
as the endpoints. As the number of abscissas increases, these endpoints approach a and b, but the subinterval width becomes smaller. While f(x) becomes evaluated at points progressively closer to a and b, and may take on progressively larger values, f(x) gets multiplied by a subinterval width that becomes progressively smaller. Hence, this approach avoids some of the problems associated with choosing ε too small. Many software packages for numerical integration use this approach to evaluate improper integrals.
Another method to evaluate improper integrals is to write the integrand f(x) as the sum of two functions and create two integrals, at least one of which is not improper so that it can be approximated directly. The other integral may be an improper integral, but it must have an analytical solution so that it can be evaluated without a numerical approximation. In the option pricing models we cover, however, most integrands are too complicated for this method to work.
Finally, integrals with limits a = −∞ or b = ∞, or both, are often encountered in financial models. One way to deal with this problem is to transform the integral by change of variable, so that the new integral is finite. In most option pricing models, however, this is not a practical solution. The more popular method is to truncate the limits of integration so that a and b are replaced by very large negative and positive numbers, respectively.
While improper integrals pose a problem with Newton-Cotes rules, this is not always the case for Gaussian quadratures. Later in this chapter we introduce Gaussian quadratures that are designed to work with improper integrals that have limits of integration (0, ∞) and (−∞, ∞).
There are two different ways of implementing numerical integration in VBA, by passing vector parameters to VBA functions and by using inline functions. We first describe implementing Newton-Cotes methods by passing vector parameters, then by using inline functions. We also compare the accuracy of each rule using a simple functional example for which the exact integral is known.
This method of implementing numerical integration refers to how the inputs needed for processing the algorithm are handled. Inputs such as the vector of abscissas or the interval width must be created in a spreadsheet. These inputs are subsequently passed to the function via the function declaration statement.
The Excel file Chapter2Newton contains VBA code to implement the Newton-Cotes rules described in this chapter. The VBA functions RPRnumint(), LPRnumint(), MPRnumint(), and TRAPnumint() construct the right-point, left-point, midpoint, and trapezoidal rules, respectively. Each function requires two vectors as inputs. For all four functions, the vector x being passed contains the abscissas for the interval [a, b], while the vector y contains values of the function to be integrated, evaluated at each element of x for the left- and right-point rules and for the trapezoidal rule, and at the midpoints of the subintervals defined by x for the midpoint rule. These four functions allow the use of subintervals of different width.
The right-point rule uses a summation involving the terms (xi − xi−1)yi, where yi ≡ f(xi), as shown in the following code for the RPRnumint() function:
Function RPRnumint(x, y) As Double
n = Application.Count(x): RPRnumint = 0
For t = 2 To n
RPRnumint = RPRnumint + (x(t) - x(t-1))*y(t)
Next t
End Function
On the other hand, the function LPRnumint() for the left-point rule uses the terms (xi − xi−1) yi−1 in the summation.
Function LPRnumint(x, y) As Double
n = Application.Count(x): LPRnumint = 0
For t = 2 To n LPRnumint =
LPRnumint + (x(t) - x(t-1))*y(t-1)
Next t
End Function
In the midpoint rule, which we code as the function MPRnumint() that appears below, the vector x being passed contains the abscissas, just as in the previous two functions for the right- and left-point rules. The vector y, however, must contain the function evaluated at the midpoints of the subintervals defined by x. Hence, the function MPRnumint() requires a different y input than the two previous functions.
Function MPRnumint(x, y) As Double
n = Application.Count(x): MPRnumint = 0
For t = 2 To n
MPRnumint = MPRnumint + (x(t) - x(t-1))*y(t-1)
Next t
End Function
The function TRAPnumint() uses the terms 
in the summation to produce the trapezoidal approximation.
Function TRAPnumint(x, y) As Double
n = Application.Count(x): TRAPnumint = 0
For t = 2 To n
TRAPnumint = TRAPnumint + 0.5*(x(t)-x(t-1))*(y(t-1)+y(t))
Next t
End Function
For this first set of four Newton-Cotes methods the VBA functions are simply creating rectangles and trapezoids of different sizes and summing them. Note that each function allows subintervals of different width to be used, so that the integration points passed in the vector x need not be equally spaced. Using more points and subintervals of smaller width in areas where the function is known to be jagged or steep, as done by adaptive quadrature methods, would lead to more precise approximations. However, as mentioned earlier in this chapter, with the computing power available today we may use a high density of abscissas everywhere in the interval [a, b] and arrive at the same accuracy with subintervals of equal width. In the next set of rules we code with VBA, only subintervals of equal width are used.
The functions sSIMPnumint(), SIMPnumint(), sSIMP38numint(), sBOOLEnumint(), and BOOLEnumint() implement Simpson’s rule, Simpson’s three-eighths rule, and Boole’s rule, respectively. We present functions for both versions of each rule, the simple version and the composite version. The simple versions are meant to be applied on individual subintervals, and the composite versions on a set of subintervals. Both versions of each rule require as inputs the vector y of function values and the interval width. The VBA functions for the three simple versions follow.
Function sSIMPnumint(deltax, y) As Double
sSIMPnumint = (1/3)*(y(1) + 4 * y(2) + y(3))*deltax
End Function
Function sSIMP38numint(deltax, y) As Double
sSIMP38numint = (3/8)*(y(1) + 3*y(2) + 3*y(3) + y(4))*deltax
End Function
Function sBOOLEnumint(deltax, y) As Double
sBOOLEnumint = (1/45)*(14*y(1) + 64*y(2) + 24*y(3)
+ 64*y(4) + 14*y(5))*deltax
End Function
The next three VBA functions are for the composite versions of these rules.
Function SIMPnumint(deltax, y) As Double
n = Application.Count(y): SIMPnumint = 0
For t = 1 To (n-1)/2
ind = (t*2) - 1
SIMPnumint = SIMPnumint + (1/3)*(y(ind)
+ 4*y(ind + 1) + y(ind+2))*deltax
Next t
End Function
Function SIMP38numint(deltax, y) As Double
n = Application.Count(y): SIMP38numint = 0
For t = 1 To (n-1)/3
ind = (t*3) - 2
SIMP38numint = SIMP38numint + (3/8)*(y(ind) + 3*y(ind+1)
+ 3*y(ind+2) + y(ind+3))*deltax
Next t
End Function
Function BOOLEnumint(deltax, y) As Double
n = Application.Count(y): BOOLEnumint = 0
For t = 1 To (n - 1) / 4
ind = (t * 4) - 3
BOOLEnumint = (1/45)*(14*y(ind) + 64*y(ind+1) + 24*y(ind+2)
+ 64*y(ind+3) + 14*y(ind+4))*deltax
Next t
End Function
It is important to have simple and composite versions of each rule available, since in many cases we do not have a number of subintervals that allows any single rule to be used exclusively. For example, specifying n = 29 subintervals would preclude the use of Simpson’s three-eighths rule directly, since this number is not divisible by three. We might therefore apply the simple version of Simpson’s rule for the first two subintervals, and apply the composite Simpson’s three-eighths rule to the remaining n = 27 subintervals.
Now that we have all the Newton-Cotes integration rules coded in VBA, we can use Excel to compare their accuracy in approximating the integral of our example function f(x) = x3 + 10. This comparison is in the Newton-cotes worksheet of the Excel file Chapter2Newton, and is illustrated in Figure 2.5. The first step is to create a column of abscissas, which appear in column B of the figure. In the adjacent column (column C), the formula for the function evaluated at each point is entered. To compare across functions, we must ensure that the number of subintervals is divisible by two for Simpson’s rule, by three for Simpson’s three-eighths rule, and by four for Boole’s rule. We will therefore integrate from a = −3 to b = 10, with a subinterval width Δx = 1, corresponding to n = 12 subintervals. This allows each rule to be implemented. To implement the midpoint rule two additional columns must be created, one containing the midpoints and the other containing the function evaluated at each of the midpoints. These appear in columns D and E of Figure 2.5, respectively.
FIGURE 2.5 Comparison of Newton-Cotes Numerical Integration Methods
To implement each approximation rule, its VBA function is invoked using the “=” sign and the function name, and the values needed for the approximation must be specified as inputs. For example, in Figure 2.5, cell G11 contains the trapezoidal approximation, which requires the x and y vectors to be passed to the TRAPnumint() function. The x values appear in the range B7 to B19 while the y values are in the range C7 to C19. Hence, in cell G11 we type
Alternately, in cell G11 we insert the TRAPnumint() function from the Excel menu, and highlight the required ranges for the x- and y-values in the dialog box that appears. Invoking Simpson’s three-eighths rule requires the y-values and the interval width Δx in cell C4 to be passed to the function. Hence, in cell G13 of Figure 2.5 we type
or we insert the SIMP38numint() function from the Excel menu and enter the required ranges in the dialog box. The other approximation rules are invoked in a similar manner.
Figure 2.5 presents a comparison of the rules covered so far. The true value of the integral is 2,616. As mentioned earlier, since this is a nondecreasing function the left-point rule produces an underapproximation (2,136), while the right-point rule produces an overapproximation (3,144). The midpoint and trapezoidal approximations are much more accurate, 2,604 and 2,640, respectively, but the two Simpson’s rules and Boole’s rule each produce approximations that are identical to the true value. The accuracy gained by using the two Simpson’s rules and Boole’s rule justifies the increased complexity of programming these rules, and the increased computing resources required to implement them. However, as we shall see later in this chapter, increasing the number of subintervals dramatically improves accuracy of the simpler rules.
The term inline means the function to be integrated resides inside the numerical integration procedure. The chief advantage of inline functions is that they are less cumbersome to use because inputs such as the x and y vectors do not have to be created in the spreadsheet—these are created internally by the function. This is an important feature that saves a lot of spreadsheet space when the number of integration points is large. The previous example the function f(x) = x3 + 10 is used to demonstrate how inline functions work. First, a VBA function for f(x) must be defined.
Function Examplef(x) As Double
Examplef = x ^ 3 + 10
End Function
The function Examplef() can now be passed to the parent integration function NumericalInt() that appears below. The parent function uses trapezoidal rule, invoked by the TRAPnumint() function.
Function NumericalInt(startPoint, endPoint, n, fName As String) As Double
Dim pass_x() As Double, pass_y() As Double
ReDim pass_x(n + 1) As Double, pass_y(n + 1) As Double
delta_x = (endPoint - startPoint) / n
For cnt = 1 To n+1
pass_x(cnt) = startPoint + (cnt-1)*delta_x
pass_y(cnt) = Run(fName, pass_x(cnt))
Next
NumericalInt = TRAPnumint(pass_x, pass_y)
End Function
Note that in order for this function to work properly, “Option Base 1” must be included at the top of the module in the VBA Editor. This specifies that all arrays begin indexing at one and not at zero.
The function requires the inputs startPoint and endPoint, corresponding to the endpoints a and b, respectively; n, the number of subintervals; and fName, which is a character input for the VBA name of the function to integrate. The variable delta_x creates the subinterval width Δx = (b − a)/n while the statement involving pass_x() creates each integration point in terms of the first endpoint of the interval of integration (a, b) and the subinterval width. The statement involving pass_y() fits f(x) = x3 + 10 to the vector of integration points. Hence, pass_x() and pass_y() are the two vectors needed to invoke the trapezoidal approximation that appears in the second-to-last statement of the code. To run this function in Excel for f(x) between a = −2 and b = 10 with n = 5 subintervals, type the following in any cell in the spreadsheet:
In this example we implement trapezoidal numerical integration, but any of the approximation rules can be implemented instead. This requires replacing TRAPnumint() in the second-to-last statement of the preceding VBA code with the function corresponding to the rule of choice. With the VBA Select statement, however, it is possible to include each integration rule within the same parent function. In the function declaration statement, the choice of rule appears as an extra parameter that must be specified when invoking the function.
Function NumericalInt(startPoint, endPoint, n, fName As String, Method) As
Double
Dim pass_x() As Double, pass_y() As Double
ReDim pass_x(n+1) As Double, pass_y(n+1) As Double
delta_x = (endPoint-startPoint)/n
For cnt = 1 To n+1
pass_x(cnt) = startPoint + (cnt-1) * delta_x
pass_y(cnt) = Run(fName, pass_x(cnt))
Next
Select Case Method
Case 1: NumericalInt = LPRnumint(pass_x, pass_y)
Case 2: NumericalInt = RPRnumint(pass_x, pass_y)
Case 3: NumericalInt = TRAPnumint(pass_x, pass_y)
Case 4: NumericalInt = SIMPnumint(delta_x, pass_y)
Case 5: NumericalInt = SIMP38numint(delta_x, pass_y)
Case 6: NumericalInt = BOOLEnumint(delta_x, pass_y)
End Select
End Function
To implement the integration of our function using Boole’s rule, for example, in any cell we type
With this version of the NumericalInt() function we are able to select between the left-point rule, the right-point rule, the trapezoidal rule, Simpson’s rule, Simpson’s three-eighths rule, and Boole’s rule, simply by passing the correct integer to the Method parameter when invoking the function. The NumericalInt() function does not include the midpoint rule function MPRnumint(), since the vector of x values to be passed is different than for the other functions. This function appers in the Excel file Chapter 2 Newton.
Up to this point we have ignored the choice of the subinterval width. The choice involves a trade-off between speed and accuracy, since a smaller width yields greater accuracy but requires a higher density of integration points and, consequently, more processing time. This section illustrates the added accuracy in narrowing the width of each subinterval by increasing the number of abscissas.
Figure 2.6 shows that the accuracy of the numerical approximation is improved when n is increased. The disadvantage is that the function needs to be evaluated at more points of integration. If the function is complicated and includes exponents, logarithms, and trigonometric functions, for example, the increase in computing time required to evaluate all the points can be substantial.
FIGURE 2.6 Increasing the Number of Subintervals for the Trapezoidal Rule
This aspect of numerical integration is crucial to understanding the benefits and weaknesses of different integration rules. We compare the accuracy of the different rules in approximating the integral of our example function f(x) = x3 + 10 from a = −2 to b = 10, brought on by increasing the number of abscissas. This is easily done using the inline functions, without having to create lengthy columns in the spreadsheet for the x and y values.
Figure 2.7 presents the results of this comparison, which are included in the increasing n worksheet of the Excel file Chapter2Newton. Cells C4 and C5 contain the interval endpoints a = −2 and b = 10, respectively. For example, to invoke Simpson’s three-eighths rule with n = 1, 200 integration points on our function defined by Examplef(), in cell E13 we type
FIGURE 2.7 Increasing Accuracy of Newton-Cotes Numerical Approximations
with the Method parameter chosen as 5. Figure 2.7 shows how the six approximation rules that can be invoked by the NumericalInt() function compare when the number of integration points is increased. When 120 abscissas or less are used, Simpson’s rule, Simpson’s three-eighths rule, and Boole’s rule are superior to the left-point, right-point, and trapezoidal rules. With 1,200 abscissas, the trapezoidal rule is very accurate also. Even with 60,000 abscissas, however, the left-point and right-point rules remain inferior to the other rules. Nevertheless, both approximations still converge to the true value of the integral as the number points increases.
From this experiment it is clear that all of the approximation rules converge to the real integral value as we increase the number of abscissas. In practice, as long as enough points are used, any approximation method can be implemented. Thankfully, there is no restriction on the number of abscissas that can be specified. Of the Newton-Cotes methods covered in this chapter, trapezoidal integration is the most commonly used in practice, since its accuracy is quite high and it is easy to implement.
While Newton-Cotes methods approximate the integral on a set of equally spaced and predetermined abscissas, Gaussian quadratures select the optimal abscissas x1, x2, …, xn at which the function should be evaluated, along with a set of weights w1, w2, …, wn. Once these abscissas and weights are obtained, the approximation is straightforward because Gaussian quadratures take the form of a summation involving the weights and the function evaluated at the abscissas.
The most difficult part of Gaussian quadratures is finding the abscissas and weights. Hence, it is often simpler to tabulate these into the VBA code, evaluate the function at the abscissas, and take the sum. In this section we present VBA code that uses tabulated values, and VBA code that finds the abscissas and weights for an arbitrary number of points.
Gauss-Legendre quadratures are constructed for integrals over [−1, 1] but can be generalized to other intervals by transformation of variables. The abscissas, which lie in [−1, 1], and weighting coefficients are chosen so that the approximation given by
yields the smallest possible error. Since there are 2n parameters involved (abscissas plus weights), choosing the abscissas as the n zeroes of the Legendre polynomial of degree n, along with the correct coefficients, yields exact solutions for polynomials up to degree 2n − 1. In most applications, 10 points or less are sufficient to yield accurate approximations.
To approximate the integral on a more general interval [a, b], the transformation y = (2x − a − b)/(b − a) is applied, so that the transformed integral has limits [−1, 1] and the Gaussian quadrature is applied to the transformed integrand. Hence, the integral (2.5) becomes
The same weights are used, but the function is evaluated at the transformed abscissas instead of at the abscissas themselves.
In this section we present three algorithms for Gaussian quadratures. The first algorithm implements Gaussian-Legendre quadrature with 10 abscissas and 10 weights tabulated within the VBA code. The second is Gauss-Legendre quadrature that does not use tabulated abscissas or weights, but rather obtains the abscissas and weights through construction of Legendre polynomials within the VBA function used to code the quadrature. The third is Gauss-Laguerre quadrature, which employs Laguerre polynomials instead of Legendre polynomials, and is designed for integrals over (0, ∞) only. In the appendix to this chapter we provide tabulated values of optimal abscissas and weights for Gauss-Legendre, Gauss-Laguerre, and Gauss-Hermite quadratures. The VBA functions in this section are adapted from Press et al. (2002).
In this section we use tabulated abscissas and weights to obtain a 10-point Gauss-Legendre quadrature. The Excel file Chapter2Gauss contains the VBA function GLquad10(), which uses our example f(x) = x3 + 10. This function requires as inputs the name of the function to be integrated, and the endpoints of the integration interval.
Function examplef(x) As Double
examplef = x ^ 3 + 10
End Function
Function GLquad10(fname As String, a, b)
Dim x(5) As Double
Dim w(5) As Double
x(1) = 0.1488743389: x(2) = 0.4333953941
x(3) = 0.6794095682: x(4) = 0.8650633666
x(5) = 0.9739065285
w(1) = 0.2955242247: w(2) = 0.2692667193
w(3) = 0.2190863625: w(4) = 0.1494513491
w(5) = 0.0666713443
xm = 0.5 * (b + a): xr = 0.5 * (b - a)
s = 0
For j = 1 To 5
dx = xr * x(j)
s = s + w(j)*(Run(fname, xm+dx) + Run(fname, xm-dx))
Next j
GLquad10 = s * xr
End Function
This VBA code is simple because it exploits the fact that the abscissas and weights in the algorithm are each symmetric about zero. Only five terms are involved in the summation that produces the approximation. If wj denotes the weight associated with abscissas xj and −xj, then the two terms corresponding to these abscissas can be grouped as
This appears as
w(j)*(Run(fname, xm+dx) + Run(fname, xm-dx))
in the function GLquad10() above.
The Excel file Chapter2Weights contains tabulated abscissas and weights to implement Gauss-Legendre quadratures, using up to 32 points.
This algorithm does not use weights and abscissas tabulated within the code, but rather uses Legendre orthogonal polynomials generated by P0(x) = 1, P1(x) = x, and Pj(x) through the recursive formula
The recursion for Legendre polynomials takes on this simple form because the weighting function is w(x) = 1. The function GLquad() is used to implement Gauss-Legendre numerical integration in VBA. This function requires as inputs the name of the function to be integrated, the endpoints, and the number of points to use in the approximation.
Function GLquad(fname As String, x1, x2, n)
Dim x(50) As Double
Dim w(50) As Double
m = (n + 1) / 2
xm = 0.5 * (x2 + x1): xl = 0.5 * (x2 - x1)
Pi = 3.141592654: EPS = 3E-16
For i = 1 To m
z = Cos(Pi*(i - 0.25) / (n + 0.5))
Do
p1 = 1: p2 = 0
For j = 1 To n
p3 = p2: p2 = p1
p1 = ((2 * j - 1)*z*p2 - (j - 1)*p3) / j
Next j
pp = n * (z * p1 - p2) / (z * z - 1)
z1 = z: z = z1 - p1 / pp
Loop Until (Abs(z - z1) < EPS)
x(i) = xm - xl * z
x(n + 1 - i) = xm + xl * z
w(i) = 2 * xl / ((1 - z * z) * pp * pp)
w(n + 1 - i) = w(i)
Next i
GLquad = 0
For i = 1 To n
GLquad = GLquad + w(i) * Run(fname, x(i))
Next i
End Function
The GLquad() function uses the roots of Chebychev polynomials, each given by 
, as initial approximations to the abscissas. It then applies Newton’s methods to find more precise abscissas. This requires the first derivative of the Legendre polynomials, given as 
and 
, and found using the recursive relation
The GLquad() function stores the abscissas in array x, stores the weights in array w, evaluates the function at the abscissas, and constructs the sum which produces the approximation to the integral.
The third and last method of Gaussian quadrature presented in this section uses Laguerre polynomials as interpolating polynomials. The function GLAU() in the Excel file Chapter2Gauss contains the VBA code for implementing Gauss-Laguerre quadrature. Since Laguerre polynomials use a weighting function of the form w(x) = xαe−x, and since this approximation is designed for integrals over the interval (0, ∞), this function requires as inputs the function to be integrated, a value for α, and the number of points to use in the approximation.
Function GLAU(fname As String, n, alpha)
Dim x(100) As Double
Dim w(100) As Double
MAXITER = 10
EPS = 0.0000000000003
For i = 1 To n
If (i = 1) Then
z = (1 + alpha) * (3 + 0.92 * alpha)
/ (1 + 2.4 * n + 1.8 * alpha)
ElseIf (i = 2) Then
z = z + (15 + 6.25 * alpha) / (1 + 0.9 * alpha + 2.5 * n)
Else
ai = i - 2
z = z + ((1 + 2.55 * ai) / (1.9 * ai) + 1.26 * ai * alpha
/ (1 + 3.5 * ai)) * (z - x(i - 2)) / (1 + 0.3 * alpha)
End If
For its = 1 To MAXITER
p1 = 1
p2 = 0
For j = 1 To n
p3 = p2
p2 = p1
p1 = ((2 * j - 1 + alpha - z) * p2 - (j - 1 + alpha) * p3) / j
Next j
pp = (n * p1 - (n + alpha) * p2) / z
z1 = z
z = z1 - p1 / pp
If (Abs(z - z1) < EPS) Then
Exit For
End If
Next its
If (its > MAXITER) Then
MsgBox “Too many iterations”
End If
x(i) = z
w(i) = -Exp(Application.GammaLn(alpha + n)
- Application.GammaLn(n)) / (pp * n * p2)
w(i) = Exp(x(i)) * w(i)
Next i
GLAU = 0
For i = 1 To n
GLAU = GLAU + w(i) * Run(fname, x(i))
Next i
End Function
Similar to the Gauss-Legendre quadrature, the GLAU() function first uses initial guesses for the roots of the Laguerre polynomial, then uses Newton’s method to refine the guesses and obtain more precise estimates of the roots. It also uses the Excel function GammaLn() to calculate the log of the Gamma function.
The GLAU() function presented in this section involves finding the roots of the Laguerre polynomials, which are subsequently used as abscissas in the approximation for the integral. Of course, it is possible to simply use tabulated values of the abscissas and weights in the VBA code directly, as in the function GLquad10() described earlier. In that case the algorithm is very simple, since it involves the sum of the product of the weights and the function evaluated at the abscissas. The Excel file Chapter2Weights contains tabulated weights and abscissas to implement Gauss-Legendre and Gauss-Laguerre quadratures, using up to 32 points. The file also includes abscissas and weights to implement Gauss-Hermite quadratures, and is described in more detail in Appendix 2.1.
In this section we compare the accuracy of the three Gaussian quadrature methods. Figure 2.8 illustrates the Excel spreadsheet Chapter2Gauss. The 10-Point Gauss-Legendre VBA function requires only the endpoints a and b to be specified, along with the function to be integrated. Hence, in cell F10 we type
FIGURE 2.8 Comparison of Gaussian Quadratures
where Examplef() corresponds to the function f(x) = x3 + 10. The Variable-Point Gauss-Legendre function requires the endpoints and the number of points to be specified, so in cell F8 we type
Because f(x) is a polynomial of order three, only two points are needed to yield the exact result of 2,616.
We also evaluate the approximation of the integral of f(x) = 1/(1 + x)3 over the interval (0, ∞), which has exact value 0.5. The VBA function Examplef2() contains f(x), so in cell F18 we type
which invokes the Gauss-Laguerre method and yields an approximation that is very close to the true value even when n = 12 points are used.
This chapter presented some of the more popular algorithms for performing numerical integration, as well as VBA functions to implement them. Newton-Cotes formulas are easy to understand, and writing VBA code for these methods is straightforward. In most applications, either the trapezoidal rule or Simpson’s rule will produce sufficiently precise approximations, provided that a large number of abscissas are used. With the speed of today’s computers, using a large number of points will not lead to substantial increases in computation or execution time.
We have introduced two Gaussian quadratures, one that can be used over any interval of integration, and one that can be used only on (0, ∞). Gaussian quadratures are much more difficult to understand than Newton-Cotes formulas, but they produce accurate results with much fewer abscissas. Even though their algorithms are more complicated, since the roots of orthogonal polynomials must be found, Gaussian quadratures produce very accurate approximations with little execution time.
Finally, we note that other Gaussian quadratures can be constructed using different weighting functions to generate orthogonal polynomials. Examples of these methods include Gauss-Hermite, Gauss-Chebychev, and Gauss-Jacobi quadratures. Since the focus of this book is not numerical integration, we leave these methods out of this chapter and refer readers to Burden and Faires (2001) or other textbooks on numerical analysis.
This section provides reinforcement of the concepts of numerical integration presented in this chapter. Solutions to these problems are in the Excel file Chapter 2 Exercises provided on the CD-ROM.
2.1 Reproduce Figure 2.5 using the following functions over the interval (−3, 3), using Δx = 0.5. Recall that Figure 2.5 uses the method of passing parameters.
2.2 Repeat Question 1 using in-line functions, using n = 12, n = 120, and n = 1, 200 subintervals in the approximation.
2.3 Write a VBA function for implementing Simpson’s rule. If the number of subintervals is odd, use Simpson’s three-eighths rule over the first three subintervals, and Simpson’s rule over the remaining n − 3 subintervals. Apply this to the functions in Question 2.1 above.
2.4 Write VBA code for the sixteen-point Gauss-Legendre quadrature, using the tabulated values for abscissas and weights presented in the appendix to this chapter, which also appears on the CD-ROM. Apply this approximation to the functions in Question 2.1 above.
2.5 Gaussian quadratures produce exact integration for polynomials of degree 2n − 1 or less, where n is the number of subintervals in the Gaussian quadrature. Use the variable-point Gauss-Legendre quadrature included on the CD-ROM to show that the quadrature yields an exact solution for the function in Question 2.1(b).
2.1 The left- and right-point rules produce very poor approximations for the integrals involving functions f1 and f3. The other rules produce approximations that are accurate to two decimals. For function f2, only Boole’s rule produces an approximation that is close to the true value. Even that approximation, however, is poor. To save space, Figure 2.9 presents the results only for function f1.
FIGURE 2.9 Solution to Exercise 2.1
2.2 This exercise illustrates that dramatically increasing the number of subintervals will in some cases produce approximations that are not adequate. Even with n = 1, 200 approximations the left- and right-point rules cannot accurately approximate the integrals for f1 and f3. Boole’s rule produces an accurate approximation for function f2 even with only n = 120 subintervals. This exercise is illustrated in Figure 2.10.
FIGURE 2.10 Solution to Exercise 2.2
2.3 The VBA function SIMPGENnumint() implements Simpson’s rule for an even or odd number of subintervals. This is illustrated in Figure 2.11.
FIGURE 2.11 Solution to Exercise 2.3
2.4 The VBA function GLquad16() adapts the Ten-Point Gauss-Legendre function GLquad10() described in the section titled “Passing Vector Parameters,” to 16 abscissas. This is illustrated in Figure 2.12.
FIGURE 2.12 Solution to Exercise 2.4
2.5 The degree of the polynomial defined in f2 is 10, so at least n = 6 abscissas are needed provide an exact solution to the integral for f2. Since 2n − 1 = 2(6) − 1 = 11, six abscissas provides an exact solution to polynomials up to degree 11 also. Note that n = 5 abscissas are not sufficient in this example, since that would provide an exact solution for polynomials up to degree 2(5) − 1 = 9 only. This exercise is illustrated in Figure 2.13.
FIGURE 2.13 Solution to Exercise 2.5
This appendix contains abscissas and weights for implementing Gauss-Legendre, Gauss-Laguerre, and Gauss-Hermite quadratures, each tabulated for 8, 16, and 32 abscissas. The weights are available in the Excel file Chapter2Weights. Implementing these rules is straightforward, since they all involve a sum whose terms are the product of the tabulated weights, and the function evaluated at the abscissas (or for Gauss-Legendre quadratures, at a transformation of the abscissas).
For integrals over the interval [a, b] we use the Gauss-Legendre points and approximate the integral by Equation (2.6). The abscissas and weights are illustrated in Figure 2.14.
FIGURE 2.14 Abscissas and Weights for Gauss-Legendre Quadratures
For integrals over (0, ∞) we use Gauss-Laguerre points and approximate the integral by
where, as before, wi and xi denote the tabulated values of the weights and abscissas, respectively. These appear in Figure 2.15.
FIGURE 2.15 Abscissas and Weights for Gauss-Laguerre Quadratures
Finally, for integrals over (−∞, ∞) we use Gauss-Hermite points and approximate the integral by
The abscissas and weights for this method appear in Figure 2.16. Note that the Gauss-Legendre quadrature (2.6) is the only one that requires the function to be evaluated at the transformed values of the abscissas. In both other quadratures, the function is evaluated at the abscissas themselves.
FIGURE 2.16 Abscissas and Weights for Gauss-Hermite Quadratures
