How to do simple numerical integration in Excel using the Gaussian Quadrature method.
Microsoft Excel is a very powerful tool for a surprising range of applications. The built-in library of mathematical functions cover many of the needs in day-to-day design task. However, integrating an arbitrary continuous function in Excel is not one of the built-in function. In many Engineering projects, one often needs to approximate integrals of a continuous function (think calculating power or energy). Fortunately there are ways around this limitation in Excel. The simplest of the methods is the so called "trapezoidal rule" integration. As the name indicates, we evaluate the function f(x) at a number of points and calaculate the total area (integral) as the sum of the areas of small trapezoids between two points. So the integral in the interval [a,b] can be simply estimated by the trapezoid area given by:
S = (b-a) * [f(a) + f(b)]/2
Repeat this for many points, hopefully evaluated at very short spaces and you can get a reasonable estimate of the integral. Simple as the method is, it's also the leat accurate of them all. It may be OK for a first-order estimate, but the accuracy can be pretty bad as demonstrated in the example further below. This is especially true when trying to integrate exponential functions which as luck would have it, are extremely common in Engineering problems (and nature in general).
Fortunately, there are many alternative methods of integration which are much more precise and are not very computationally intensive either. Most people learn about the Simpson rule for example in an introductory numerical anlysis class. It's an improvement over the trapezoidal method but not that accurate either. Recently, I learned about a class of methods knows as "Gaussin Quadrature" that fit the bill perfectly: simple and quite accurate for the effort involved!
Methods such as the Simpson rule estimate an integral by evaluating a function f(x) at a point equidistant from a and b in each interval [a,b]. For example, with the Simpson rule, the function needs to be evaluated at x = (a+b)/2 which is exactly the middle of the interval. There is really no reason why it needs to be so. The beauty of Gaussian integration is that is does not fix the evaluation points at a fixed fraction of the interval, but rather it selects them at an optimum point (or points) in the interval. Figure 1 shows the two-point Gaussian quadrature for example. Note that the number of points can be much larger; the more points the more accurate the estimation. However, for the purpose of this article we will stick to two-point quadrature. As stated in the first equation, the integral is estimated by evaluating f(x) at two points: x = x1 and x = x2 and multiplying it by two constants c1 and c2. These points are chosen to yield minimum error and the associated formulas are shown below:
Figure 1 - The Two-point Gaussian Quadrature
To illustrate this, I made an example spreadsheet available for download. The function I chose to integrate (exponential) is a trivial case since the integral can be easily evaluated analytically:
However, this allows us to compare the integration results using both the trapezoidal method and the Gaussian quadrature against the know integral value. Once you understand the procedure, you will see that it is quite easy to adapt for any other function. The spreadsheet can be download through the link below:
|Gaussian Quadrature Example, Rev 1|
The worksheet shown in Figure 2 contains the exp(x) function evaluated at a number of points. Notice that the points are not linearly distributed but rather geometrically (xn = k * xn-1). The number of points per decade can be entered in the "# Points / Decade" field. For comparison, the integral is evaluated using both the trapezoidal and the two-point Gaussian quadrature methods. Notice the the points x = x1 and x = x2 are calculated in columns G and H.
Figure 2 - Gaussian Quadrature Spreadsheet
Since for this simple integral the exact value can also be easily calculated, we can compare the error obtained using these two approaches. Figure 3 shows the error percentage comparison as a function of x.
Figure 3 - Error Comparison
As evidenced by Figure 3, the Gaussian quadrature does a much better job of estimating the integral than the simple trapezoidal method. The trapezoidal error grows fast when it tries to approximate a convex fast-growing function such as exp(x) (i.e. in this case, it overestimates the area of the integral).
I hope this spreadsheet can be useful to others whenever simple and relatively accurate integral estimation is needed. For a more useful example application of this method, please have a look at my "Transistor Amplifier Calculator". In that spreadsheet, the Gaussian quadrature method is used to integrate the PSD function over the frequency range specified and thus calculate the total noise power.
As mentioned earlier, if even greater accuracy is needed, one can use more points in the Gaussian quadrature. In those cases it may be worth developing a VBA function as it can become tedious to add the points x1 ... xn in separate spreadsheet columns as I've done here. This may be the topic for a future article:)
Comments, questions, suggestions? You can reach me at: contact (at sign) paulorenato (dot) com