Integrating a Function

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

Figure 7-10. Some results returned by the IntegrateT custom function, (folder 'Chapter 07 Examples', workbook 'Integration', sheet 'Trapezoidal Integration Fn')

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.

Gaussian Quadrature

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')

Distance Traveled Along a Curved Path

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

Length Segment Circle
Figure 7-15. Approximating the distance along a curve AB by the length of the straight line segment AB. (folder 'Chapter 07 Examples', workbook 'Curve Distance', sheet 'Curve Distance (Circle)')

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 %.

Problems

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

+1 0

Responses

Post a comment