Implementing the Trapezoidal Rule in VBA

Problem

You see how easy it is to apply the trapezoidal rule in a spreadsheet, but you'd like to implement the rule in VBA to save having to set up the columns of data as was done in the previous recipe. You'd also like to be able to quickly change the spacing between x-values to improve accuracy.

Solution

Add two custom VBA functions as discussed here. See Recipe 2.2 to learn the basics of writing VBA functions if you don't know them already.

Discussion

Let's reconsider the example from Recipe 10.1. This time, instead of setting up the table of data shown in Figure 10-1 and using the SUMPRODUCT function in cell C20, you need only write two simple VBA functions.

The first function we need is simply a function to compute a y-value for a given x-value using the known analytic function under consideration. Example 10-2 shows just such a function for computing the analytic function from the previous example.

Example 10-2. VBA function for computing y = f(x)

Public Function F(x As Double) As Double
 F = Exp(-(x ^ 2))
End Function

As you can see, it's a fairly simple function. It takes a given x-value as an argument and returns the corresponding y-value using the given analytic function. The other required VBA function is shown in Example 10-3.

Example 10-3. VBA function for trapezoidal rule

Public Function Trapezoidal(xMin As Double, xMax As Double, n As Integer) As Double
 Dim dx As Double
 Dim sum As Double
 Dim y As Double

 dx = (xMax - xMin) / n

 sum = (F(xMin) + F(xMax)) / 2#

 For i = 1 To (n - 1)
 y = F(xMin + (dx * i))
 sum = sum + y
 Next i

 sum = sum * dx

 Trapezoidal = sum
End Function

This function takes three arguments: the lower and upper bounds of the interval of integration, xMin and xMax, and the number of intervals over which to divide the range of x-values, n.

The first task upon entering this function is to declare a few local variables. dx will store the computed spacing of the x-values given the number of intervals, n. The computed y-value is stored in the variable y. sum is used to accumulate the products of the coefficients and the y values.

After declaring these local variables, the function goes on to compute dx, which is just the x-range divided by the number of intervals. sum is also initialized to the average of the results of the function evaluated at the upper and lower bounds of integration. This step takes care of the first and last terms in the trapezoidal rule formula involving the 1/2 coefficient. Next, the function enters a For loop on the index variable i from 1 to n-1 to accumulate the function evaluated at each computed x-point. This step takes care of all the middle sums in the trapezoidal rule formula.

The first calculation within the For loop evaluates y for a given x value by calling the function F as defined in Example 10-2. Here, x is computed by multiplying dx and i. Upon exiting the loop, sum is multiplied by dx to yield the final result.

To actually use these functions in a spreadsheet you need only type the formula =trapezoidal(0,1,10) into any given cell. Here, I used the same range and number of samples as used in the example from Recipe 10.1. As expected, the result in this case is exactly the same as before, that is, 0.746.

Using the VBA approach discussed here allows you to easily change the number of samples used in the computation to improve the accuracy of the result. Increasing the argument n results in more samples at a smaller x-spacing, resulting in greater accuracy. You can also change the upper and lower bounds of integration very easily using this approach. You can achieve both of these changes by taking the spreadsheet approach, but it requires you to reconstruct your columns of data.

You can integrate other analytic functions using the VBA functions shown here. All you need to do is change the F function to use your new analytic expression; the function trapezoidal remains unaffected.

I mentioned in Recipe 10.1 that some practitioners apply an iterative approach using the trapezoidal rule to adjust the number of sample intervals until the end result stops changing by an amount greater than some given tolerance. You can easily implement such a scheme using the VBA functions shown in this recipe. As an illustration, I added another custom VBA function as shown in Example 10-4.

Example 10-4. IterateTrapezoidal function

Public Function IterateTrapezoidal(xMin As Double, xMax As Double, n As Integer, tol As
Double) As Double
 Dim value As Double
 Dim previousValue As Double
 Dim difference As Double
 Dim counter As Integer
 Dim nNew As Integer

 previousValue = Trapezoidal(xMin, xMax, n)
 difference = 99999
 counter = 0
 nNew = n

 While ((counter < 100) And (difference > tol))
 nNew = nNew * 2
 previousValue = value
 value = Trapezoidal(xMin, xMax, nNew)
 difference = Abs(value - previousValue)
 counter = counter + 1
 Wend

 IterateTrapezoidal = value
End Function

This function, called IterateTrapezoidal, takes four arguments. The first three are the same as discussed earlier for the function trapezoidal. The fourth argument is new and represents the tolerance, tol, which tells the iterative algorithm when to stop.

The scheme goes like this: first, we compute a result using TRapezoidal just as we did earlier, using the supplied default number of sample intervals, n. Then we double n, recompute the result, and compare it with the previous result to see how much the result changed given the change in n. If the change is greater than tol, then we double n again and recompute. This repeats until the difference is less than tol, at which time the loop is exited and the last computed result is returned. This algorithm works pretty well and is one way to automate the process of choosing an appropriate n when applying the trapezoidal rule.

You do have to be careful with this approach. You'll notice that I put a counter in the While loop shown in Example 10-4. You could encounter a situation where the specified tolerance is too low and the solution never converges. The counter will eventually kick execution out of the loop, to prevent the function from getting stuck in an infinite loop. You do also have to keep in mind round-off error and precision. Basically, don't specify an unrealistically small tolerance, and keep in mind the magnitude of the result you're expecting. If you're expecting a result on the order of 10,000 ± 100, a tolerance of 1e-5 is probably way too small!

See Also

Refer back to Recipe 10.1 for more discussions on the trapezoidal rule.