Instead of finding the area under a curve defined by a set of data points, you may wish to integrate a function F(x). You could simply create a table of function values and use one of the methods described in earlier sections to calculate the area. But a more convenient solution would be to create a custom function that uses the Formula property of the cell to get the worksheet formula to be integrated, in the same way that was used in the preceding chapter, and uses the formula to find the area under the curve. This approach will be described in subsequent sections.
Integrating a Function
Defined by a Worksheet Formula by Means of a VBA Function Procedure
In this section, the trapezoidal and Simpson's rule methods are implemented as VBA custom functions, using an approach similar to that used in the differentiation functions of the previous chapter. The Formula property of the cell is used to get the worksheet formula to be differentiated into the VBA code as text. Then the SUBSTITUTE worksheet function is used to replace the variable of interest by an incremented value, and the Evaluate method used to get the new value of the formula. These values are used to calculate the area of each panel, and the areas of the panels are summed to obtain the area under the curve.
This function procedure can be used to integrate an expression F(x) defined by a worksheet formula, between specified lower and upper limits a and b respectively. A table of function values is not required.
The syntax of the function is In teg rate (express/on, variable, fromjower, to_upper). The argument expression is the integrand, the expression to be integrated. The argument variable is the variable of integration, and the arguments fromjower and to_upper are the lower and upper limits of integration, respectively. The VBA code is shown in Figure 7-8. Function procedures for both trapezoidal (IntegrateT) and Simpson's rule (Integrates) methods are shown.
The range of x values over which the integration is to be performed (to_upper - fromjower) is divided into N panels. The user can adjust the accuracy of the integration by changing the value of N in the procedure, with a concomitant increase in calculation time.
Option Explicit
Function lntegrateT(expression, variable, fromjower, to_upper) 'Simple trapezoidal area integration
Dim FormulaString As String, T As String, Xref As String
Dim H As Double, area As Double, X As Double
Dim N As Integer, K As Integer, J As Integer
Dim NRepI As Integer
Dim temp As String
Dim F1 As Double, F2 As Double
FormulaString = expression.Formula
T = Application.ConvertFormula(FormulaString, xlA1, xlA1, xIAbsolute) XRef = variable. Address
H = (to_upper - fromjower) / N area = 0 X = fromjower
NRepI = (Len(T) - Len(Application.Substitute(T, XRef,""))) / Len(XRef) For K = 1 To N
For J = NRepI To 1 Step -1 temp = Application.Substitute(T, XRef, X & "", J) If lsError(Evaluate(temp)) Then GoTo pt1 T = temp pt1: Next J F1 = Evaluate (T)
T = Application.ConvertFormula(FormulaString, xlA1, xlA1, xlAbsolute) For J = NRepI To 1 Step -1 temp = Application.Substitute(T, XRef, X + H & "", J) If lsError(Evaluate(temp)) Then GoTo pt2 T = temp pt2: Next J F2 = Evaluate(T)
IntegrateT = area
End Function_
Figure 7-8. VBA Function procedure to integrate a worksheet formula by the trapezoidal approximation method, (folder 'Chapter 07 Examples,' workbook 'Integration,' module 'Simplelntegration')
Function lntegrateS(expression, variable, fromjower, to_upper) 'Simpson's 1/3 rule area integration
Dim FormulaString As String, T As String, Xref As String Dim H As Double, area As Double, X As Double Dim N As Integer, K As Integer, J As Integer Dim NRepI As Integer Dim temp As String
Dim Y0 As Double, Y1 As Double, Y2 As Double
FormulaString = expression.Formula XRef = variable.Address
T = Application.ConvertFormula(FormulaString, xlA1, xlA1, xIAbsolute) NRepI = (Len(T) - Len(Application.Substitute(T, XRef,""))) / Len(XRef) For J = NRepI To 1 Step -1 temp = Application.Substitute(T, XRef, fromjower + X & "", J) If lsError(Evaluate(temp)) Then GoTo pt1 T = temp pt1: Next J
T = Application.ConvertFormula(FormulaString, xlA1, xlA1, xIAbsolute) For J = NRepI To 1 Step -1 temp = Application.Substitute^, XRef, fromjower + X + H & "", J) If lsError(Evaluate(temp)) Then GoTo pt2 T = temp pt2: Next J Y1 = Evaluate (T)
T = Application.ConvertFormula(FormulaString, xlA1, xlA1, xIAbsolute) For J = NRepI To 1 Step -1 temp = Application.Substitute(T, XRef, fromjower + X + 2 * H & "", J) If lsError(Evaluate(temp)) Then GoTo pt3 T = temp pt3: Next J Y2 = Evaluate (T)
Integrates = area
End Function_
Figure 7-9. VBA function procedure to integrate a worksheet formula by Simpson's method, (folder 'Chapter 07 Examples', workbook 'Integration', module 'Simplelntegration')
Some results returned by the IntegrateT and Integrates functions are shown in Figures 7-10 and 7-11, respectively. In general, results are more accurate when using the Simpson's method function.
A |
B |
C |
D |
E |
F G |
^ H I-';] | |
1 |
Integration |
function using simple trapezoidal appro*. | |||||
2 |
Function |
X |
Fix) |
from |
to |
Area exact |
% error |
3 |
X3 |
1 |
1 |
0 |
1 |
0.2500003 0.25 |
1.0E-04 |
4 |
4<0-r) |
1 |
0 |
0 |
1 |
3.14156 3.14159 |
1 2E-03 |
5 |
triangle' |
1 |
2 |
0 |
2 |
4.00000000 4 |
O.OE+OO |
6 |
Gaussian |
95 |
0.035 |
20 |
180 |
1.00000000 1 |
3.1E-12 |
7 | |||||||
8 |
* |
slope |
= 2, intercept= 4 | ||||
9 |
** ¡I |
= 100 |
,o= 10 |
A |
8 |
C |
D |
E |
F | G |
H | ||
1 |
Integration function using Simpson's method | |||||||
2 |
Function |
X |
FOO |
from |
to |
Area |
exact |
% error |
3 |
x3 |
1 |
1 |
0 |
1 |
0,2500000 |
0,25 |
4.4E-1 4 |
4 |
4-iCI-x2) |
1 |
0 |
0 |
1 |
3.141588 |
3.141593 |
1.6E-04 |
5 |
triangle |
1 |
2 |
0 |
2 |
4.0000000 |
4 |
2.2E-14 |
6 |
Gaussian |
95 |
0.035 |
20 |
180 |
1.0000000 |
1 |
1.3E-13 |
* ' | ||||||||
8 |
* |
slope |
= 2, intercept = 4 | |||||
9 |
** ^ |
= 100 |
, u= 10 |
Figure 7-11. Some results returned by the Integrates custom function, (folder 'Chapter 07 Examples', workbook 'Integration', sheet Simpson Integration Fn')
Figure 7-11. Some results returned by the Integrates custom function, (folder 'Chapter 07 Examples', workbook 'Integration', sheet Simpson Integration Fn')
Because some functions may require a large number of iterations, there may be a noticeable delay in calculation.
The preceding methods for numerical integration employ evenly spaced values of x at which the function is evaluated. Other formulas have been developed whereby the function is evaluated at specially selected values of x. These Gaussian quadrature formulas are significantly more efficient, in terms of the accuracy of the evaluation.
Gaussian quadrature formulas involve the evaluation of the function at a set of Xj values (nodes), with the use of a corresponding set of weights w„ in the following formula
The nodes and weights can be derived from certain kinds of polynomials. The Legendre polynomials will be used here to determine the values of x{ and Wj. The Legendre polynomials are a set of polynomials of degree N. Increasing N provides an increase in accuracy of evaluation but requires a concomitant increase in computation time. Values of Legendre polynomials for N up to 100 have been published.
The integration need not be limited solely to the interval -1 to 1. By employing a change of variable
the integral expression is b , 1 , A= JF(x)dx =j>
and equation 7-9 becomes b b-a^ J (b - a)zi +(b + a)
a which permits integration over any range.
The code shown in Figure 7-12 performs Gaussian quadrature using equation 7-12 and a tenth-order Legendre polynomial. Some results returned by the function are shown in Figure 7-13.
Option Explicit
Function lntegrate(expression, variable, fromjower, to_upper, Optional _ tolerance)
Dim FormulaString As String, XAddress As String Dim result As Double
FormulaString = expression.Formula XAddress = variable.Address 'Default is absolute FormulaString = Application.ConvertFormula(FormulaString, xlA1, xlA1, _ xIAbsolute)
Call GaussLeGendre10(FormulaString, XAddress, fromjower, to_upper, _
tolerance, result) Integrate = result End Function
Sub GaussLeGendre10(expression, XRef, fromjower, to_upper, tolerance, result)
'Uses ten-point Gauss-Legendre quadrature formula. 'Adapted from Shoup, p.203
Dim XJ As Variant, AJ As Variant
Dim TotalArea As Double, OldArea As Double, area As Double Dim T As String, temp As String
Dim I As Integer, J As Integer, K As Integer, JJ As Integer Dim N As Integer, NRepI As Integer
Dim A As Double, B As Double, C As Double, D As Double, F As Double Dim H As Double
XJ = Array(-0.973906528517172, -0.865063366688984, -0.679409568299024, -0.433395394129247, -0.148874338981631, 0.973906528517172,
0.865063366688984, 0.679409568299024, 0.433395394129247, 0.148874338981631) AJ = Array(0.066671344308688, 0.149451349915058, 0.219086362515982, 0.269266719309996, 0.295524224714753, 0.066671344308688, 0.149451349915058, 0.219086362515982, 0.269266719309996, 0.295524224714753)
If IsMissing(tolerance) Then tolerance = 0.0000000001 OldArea = 0 N = 1
For K = 1 To 10 'increments divided by 1,2,4,8,16,32,64,128,256,512 area = 0
A = from Jower + (I -1) • H B = A + H C = (B + A) / 2 D = (B - A) / 2
NRepI = (Len(T) - Len(Application.Substitute(T, XRef,""))) / Len(XRef) For J J = NRepI To 1 Step-1 temp = Application.Substitute^", XRef, C + D * XJ(J) & "", JJ) If lsError( Evaluate(temp)) Then GoTo pt1 T = temp pt1: Next JJ F = Evaluate(T) area = area + AJ(J) * F Next J Next I
If Abs((area - OldArea) / area) < tolerance Then GoTo AIIDone
OldArea = area
Next K
AIIDone:
result = area
End Sub
Figure 7-12. Integrate custom function, (folder 'Chapter 07 Examples', workbook 'Integration', module 'Legendreintegration')
A |
B I |
c |
D |
E |
F i G |
.H | |
1 |
Custom Function to Integrate Aiea Under a Curve | ||||||
2 |
-Integ rate {formu5a_in_ceS5, variable, from_ lower; fo_ upper, tolerance) | ||||||
3 |
Function |
X |
FOO |
ftom |
to |
area exact value |
% error |
4 |
4V(1 ->r)dx |
1 |
0 |
-1 |
1 |
3.1415928 3.1415927 |
4.0E-06 |
S |
x\jx |
1 |
1 |
0 |
1 |
0.2500000002 0,25 |
7.SE-08 |
6 |
(mx+b)dx |
2 |
0 |
0 |
2 |
4.000000003 4 |
7.6E-08 |
7 |
sin(x)dx |
o |
0 |
0 |
3.14 |
2.000000001 2 |
7.4E-08 |
9 |
(x e"')dx |
1 |
0.37 |
0 |
1000 |
1000000001 1 |
9.5E-08 |
9 |
(1A0 dx |
1 |
1 |
1 |
2 |
0.69314718 0.69314718 |
7.7E-08 |
Figure 7-13. Some results returned by the Integrate custom function, (folder 'Chapter 07 Examples', workbook 'Integration', sheet 'GaussLegendre Integration Fn')
Figure 7-13. Some results returned by the Integrate custom function, (folder 'Chapter 07 Examples', workbook 'Integration', sheet 'GaussLegendre Integration Fn')
Early versions of this program returned inaccurate results when the range h — a was large. The function Integrate illustrates one approach to overcoming this problem. First, the integral is evaluated over the total range b - a. Then the interval is divided into two halves and each "panel" is integrated separately. The sum of the two panels is compared to the previous value. If the difference is larger than a tolerance value, the interval is divided into quarters, the areas summed and so on. The process is continued for 10 cycles of iteration (512 panels) or until the area difference is less than a specified tolerance.
Because some functions may require a large number of iterations, there may be a noticeable delay in calculation. Increasing the value of tolerance should speed up calculation, but only at the expense of accuracy.
Integration with an Upper or Lower Limit of Infinity
Integrals such as
a can be evaluated by summing the areas of a number of panels covering the range from x = a to x = a suitably large value. It is to be expected that as x—> oo the area of panel(x) -» zero. Thus the integral can be evaluated by summing the integrals of a series of panels of increasing width (e.g., from 0-1, 1-10, 10-100, etc), ending the summation when the area of the last panel is suitably small. Manual adjustment of the panel widths is easily done by inspection of the results. Figure 7-14 shows a typical result.
1 A „ |
B |
C |
D |
E |
F i |
G |
H ! | ||
32 | |
Integrating to upper limit of infinity | ||||||||
33 |
Function |
X |
m |
from |
to |
integrated value |
I |
% error | |
34 |
y = 1 ¿£(1 +x)*sqtt(x) |
1 0.5 |
0 |
0.01 |
0.1989708 : |
0.1989708 |
93.67 | ||
35 |
1 |
0.5 |
0.01 |
01 |
0.4132174 |
0.6121882 |
80.51 | ||
36 |
1 |
0.5 |
IM: |
1 |
0.9582416 |
1.5704298 |
50.01 | ||
37 |
1 |
0.5 |
1 |
10 |
0.9582416 |
2.5286714 |
19 51 | ||
36 |
1 I 0.5 |
10 |
100 |
0.4132174 |
2.9418888 |
6.36 | |||
39 |
1 |
0.5 |
100 |
1000 |
0.1361128 |
3.0780016 |
2.02 | ||
40 |
0.5 |
1E+03 |
1E+04 |
0.0432252 |
3.1212268 |
0.65 | |||
41 |
1 |
0.5 |
1E+04 |
1E+05 |
0.0136748 |
3.1349016 |
0.21 | ||
42 |
1 |
0.5 |
1E+05 |
1E+06 |
0.0043245 |
3.1392261 |
0.08 | ||
43~ |
1 |
0.5 |
1E+06 |
1E+08 |
0.0018000 |
3.1410261 |
0.02 | ||
44 |
1 |
0.5 |
1E+08 |
1 E+10 |
0.0001800 |
3.1412061 |
0.01 | ||
45 |
1 |
0.5 |
1E+10 |
1E+14 |
0,0000198 |
3.1412259 |
0.01 | ||
46 |
1 |
0.5 |
1E+14 |
1E+18 |
0,0000002 |
3.1412261 |
0.01 | ||
47 |
exact value |
3.1415927 |
Figure 7-14. Integrating from a lower limit to an upper limit of infinity. Results returned by the Integrate custom function, (folder 'Chapter 07 Examples', workbook 'Integration', sheet 'Integrating to infinity by sum')
The length of a plane curve can be estimated by dividing the curve into segments, as in Figure 7-15, and approximating the length of the curve segment by the straight line AB. The length of AB = V(Ax)2 + (Ay)2 . The distance along
A |
B _ |
c | |
1 |
Distance Travelled Along a Curve | ||
2 |
1/4 circle with r = |
= 1 | |
3 |
(distance should be n/2) | ||
4 |
X |
y |
d |
5 ' |
0.000 i 1.000 | ||
6 |
0.050 |
0.999 |
0.050 |
7 ' |
0 075 |
0.997 |
0.G25 |
8 |
0 100 |
0.995 |
0.025 |
9 |
0.125 |
0.992 |
0.025 |
51 |
0.980 |
0.199 |
0.024 |
52 |
0.985 |
0.173 |
0.027 |
53 |
0.990 |
0.141 |
0.032 |
54 |
0.9925 |
0.122 |
0.019 |
55 |
0.9950 |
0.100 |
0.023 |
56 |
0.9975 |
0.071 |
0.029 |
57 |
0.999 |
0.045 |
0.026 |
58 |
1.000 |
0.000 |
0 045 |
59 |
I Sumx2 = |
3.14145 | |
60 |
% error = |
4 JE-03 |
Figure 7-16. Approximating the circumference of a circle of radius 1.
Note that the rows between 9 and 51 are hidden, (folder 'Chapter 07 Examples', workbook 'Curve Distance', sheet 'Curve Distance (Circle)')
Figure 7-16. Approximating the circumference of a circle of radius 1.
Note that the rows between 9 and 51 are hidden, (folder 'Chapter 07 Examples', workbook 'Curve Distance', sheet 'Curve Distance (Circle)')
The procedure is illustrated by estimating the length of one quarter of a circle of radius r= 1. The equation of the circle is x2 + y2 = 1, or y =-%\-x2 . As shown in Figure 7-16, the value of y and the distance d between successive points was calculated from x = 0 to * = 1, using an x increment of 0.025. Near the end of the range of x values, where y changes more rapidly, the x increment was decreased. The formula in cell C6 is
=SQRT((A8-A5)A2+(B8-B5)A2) The sum of the distances x 2, in cell C59 is a reasonable estimate of %.
Answers to the following problems are found in the folder "Ch. 07 (Integration)" in the "Problems & Solutions" folder on the CD.
00 x2
method.
2. Integrate the following expressions, using one of the custom functions for integration.
o lnx
3. Evaluate the elliptic integral nil j\/l - (1/2)sin2 xdx
4. An ellipse is a plane figure described by the locus of a point P(x, y) that moves such that the sum of its distances from two fixed points (foci) is a constant. If the ellipse has foci located at A (-c, 0) and B (c, 0) and the distance ACB is 2a, then by setting b = Va2 -c2 , the equation of the ellipse is simplified to
(a and b are termed the semiaxes of the ellipse).
Figure 7-17. Approximating the circumference of an ellipse.
For the ellipse shown in Figure 7-17, with foci at x = -0.5, y = 0 and x = 0.5, y = 0 and a= 1, determine the circumference of the ellipse.
5. Determine the area of the ellipse of problem 7-4.
6. Find the area between the curve y - 2x - x2 and the line y = -3.
7. Find the area between the curve y = 2x - xz and the line y = 2.5x -2.3.
8. Find the area enclosed between the two curves shown in Figure 7-7: y\ = x3 -
20x2 - IOOjc + 2000 and_y2 = 2x3 - 5x2 - 300x + 1000. The curves intersect in the region between x = -5 and x = 15.
(a and b are termed the semiaxes of the ellipse).
Figure 7-17. Approximating the circumference of an ellipse.
9. The area between the curve y = x2 and the horizontal line y = 4 is divided into two equal areas by the horizontal line y — c. Find c.
10. The area between the curve y = x2 + 3 and the line y = 12 is divided into two equal areas by the line y — c. Findc.
11. Integrate the following expression.
12. Integrate the following expressions, using the custom function for integration.
00 j
This Page Intentionally Left Blank
Was this article helpful?