CH-05 (Pao - Engineering Analysis)

5Numerical Integrationand Program Volume5.1 INTRODUCTIONSometimes, one cannot help wonder why appears so often in a wide range ofmathematical problems and why it has a value of approximately equal to 3.1416.One may want to calculate this 16th letter in the Greek alphabet and would like toobtain its value as accurate as 3.14159265358979 achieved by the expert.1 Geometrically, represents the ratio of the circumference to the diameter of a circle.

It iscommonly known that if the radius of a circle is r, the diameter is equal to 2r, thecircumference is equal to 2r, and the area is equal to r2. Hence, to calculate thediameter we simply double the radius but to calculate the circumference and thearea of a circle is more involved. The transcendental number is the result ofcalculating the circumference or area of a circle by numerical integration.In this chapter, we discuss various methods that can be adopted for the need ofnumerical integration. Before we elaborate on determining the value of , let usdescribe the problem of numerical integration in general.Consider the common need of finding the area inside a closed contour C1 suchas the one shown in Figure 1A, or the area between the outside contour C2 and theinside contour C3 shown in Figure 1B.

The latter could be a practical problem ofdetermining the usable land size of a surveyed lot which include a pond. To evaluatethe area enclosed by the contour C1 approximately by application of digital computer,the contour can be treated as two separate curves divided by two points situated atits extreme left and right, denoted as PL and PR, respectively. A rectangular coordinatesystem can be chosen to adequately describe these two points with coordinates(XL,YL) and (XR,YR). Here, for convenience, we shall always place the entire contourC1 in the first quadrant of the X-Y plane. Such an arrangement makes possible tohave the coordinate (X,Y) values any point on C1 being greater than or equal to zero.The area enclosed in the contour C1 can be estimated by subtracting the area ABbetween the bottom-branch curve PLPBPR and the X-axis, from the area AT betweenthe top-branch curve PLPTPR and the X-axis.

Approximated evaluation of the areasAB and AT by application of digital computer proceeds first with selection of a finitenumber of points Pi from PL to PR. That is, to approximate a curve such as PLPTPRby a series of linear segments. Let N be the number of points selected on the curvePLPTPR, then the coordinates of a typical point are (Xi,Yi) for i ranges from 1 to Nand in particular (X1,Y1) = (XL,YL) and (XN,YN) = (XR,YR). The area between atypical linear segment PiPi + 1 and the X-axis is simply equal to:A i = (Yi + Yi +1 ) ( X i +1 − X i ) 2© 2001 by CRC Press LLC(1)FIGURE 1A.

The common need of finding the area inside a closed contour C1.FIGURE 1B. The common need of finding the area between the outside contour C2 and theinside contour C3.© 2001 by CRC Press LLCFIGURE 2. (Yi + Yi + 1)/2 is the average height and (Xi + 1–Xi) is the width of the shaded strip.Notice that (Yi + Yi + 1)/2 is the average height and (Xi + 1–Xi) is the width of theshaded strip shown in Figure 2.

Obviously, the total area AT between the top branchof contour C1, PLPTPR, and the X-axis is the sum of all strips under the N–1 linearsegments PiPi + 1 for i = 1,2,…,N. In other words, we may mathematically write:N −1AT =N −1∑ A = ∑ (Y + Yii =1ii +1i =1) (Xi+1 − Xi ) 2(2)To obtain the area AB between the bottom branch of contour C1 and the X-axis,we follow the same procedure as for the area AT except that the first point is to beassigned to PR and the last point to PL.

Suppose that there are M points selectedalong PRPBPL, then the coordinates of these points are (Xi,Yi) for i = 1,2,…,M andin particular (X1,Y1) = (XR,YR) and (XM,YM) = (XL,YL). Consequently, the area ABcan be calculated, similar to Equation 2 as:M −1AB =∑i =1© 2001 by CRC Press LLCM −1Ai =∑ (Y + Yii =1i +1) (Xi+1 − Xi ) 2(3)Since the points are numbered in increasing order from PR through PB and PL,it is then clear that Xi + 1 is always less than Xi.

AB thus carries a minus sign.Based on the above discussion, the area enclosed by contour C1 can thereforebe calculated by adding AT and AB if the numbering of the points selected on thecontour follows a clockwise direction. Let the total number of points selected aroundthe contour C1 be denoted as K, then K = N + (M–2) because PR and PL are re-usedin consideration of the bottom branch. Hence, the area enclosed in C1 is:KA = AT + AB =∑ (Y + Yii =1i +1) (Xi+1 − Xi ) 2(4)where the Nth point has coordinates (XN,YN) = (XR,YR) and the first and last pointshave coordinates (X1,Y1) = (XK + 1,YK + 1) = (XL,YL).

And it should be evident thatin case of a cut-out, such as the contour C2 shown in Figure 1(B), the subtractionof the area enclosed by the cut-out can be replaced by an addition of the value ofthe area when it is calculated by using Equation 4 but the numbering of the pointson contour C3 is ordered in counterclockwise sense.5.2 PROGRAM NUINTGRA — NUMERICAL INTEGRATION BYAPPLICATION OF THE TRAPEZOIDAL AND SIMPSON RULESProgram NuIntGra is designed for the need of performing numerical integrationby use of either trapezoidal rule or Simpson’s rule. These two rules will be explainedlater.

First, let us discuss why we need numerical integration.Figure 3 shows a number of commonly encountered cross-sectional shapes inengineering and scientific applications. The interactive computer program NuIntGrahas an option of allowing keyboard input of the coordinates of the vertices of thecross section and then carrying out the area computation of cross-sectional areabased on Equation 4.The following gives some detailed printout of the results for the cross sectionsshown in Figure 3. It is important to point out that the points on the contours describingthe cross-sectional shapes should be numbered as indicated in Figure 3, namely, clockwise around the outer boundary and counterclockwise around the inner boundary.© 2001 by CRC Press LLCFIGURE 3.

Commonly used cross sections in engineering and scientific applications.© 2001 by CRC Press LLC© 2001 by CRC Press LLCBy use of a Function subprogram F(X) which defines the upper branch of acircle of radius equal to 1 and having its center located at X = 1 as listed below,program NuIntGra also has been applied for calculating the value of . The screendisplay of this interactive run is also listed below after the Function F(X).QuickBASIC VersionFORTRAN VersionFUNCTION F(X)F = SQR(1(X-1)^2)END FUNCTIONFUNCTION F(X)F = SQRT(1.-(X-1)**2)RETURNEND© 2001 by CRC Press LLCMore discussion on the accuracy of will be given later after we have introducedboth the trapezoidal and Simpson’s rules which have already been incorporated inthe program NuIntgra.TRAPEZOIDAL RULEReturning to Figure 2, we notice that if in approximating the curve by a seriesof linear segments the points selected on the curve are equally spaced in X, Equation2 can be considerably simplified.

In that case, we have:∆X = X 2 − X1 = X 3 − X 2 = … = X N − X N −1(5)and Equation 2 can be written as:N −1A T = ∆X∑ (Y + Yii =1i +1)2(6)Or, in a different form for easy interpretation, it may also be written as:N −1A T = ∆X Y1 + 2Yi + YN  2i=2∑(7)All the in-between heights, Yi for i ranging from 2 to N–1 that is the next to thelast, are appearing twice because they are shared by two adjacent strips whereas thefirst and last heights, Y1 and YN, only appear once in Equation 7.

X in Equation 7is the common width of all strips used in summing the area.Equation 7 is the well known Trapezoidal Rule for numerical integration. In ageneral case, it can be applied for approximate evaluation of an integral by theformula:© 2001 by CRC Press LLC∫XRXLN −1f ( X)dX = ∆X f1 + 2fi + fN i=2∑(8)where the increment in X, X, is simply:∆X = ( X L − X R ) ( N − 1)(9)when N points are selected on the interval of integration from XL to XR.

It shouldbe understood that in Equation 8 fi is the value of the integrand function f(X)calculated at X = Xi. That is:fi ≡ f ( X i )(10)X i = (i − 1)∆X + X L(11)where for i = 1,2,…,NX1 = X LandXN = XR(12,13)Program NuIntGra allows the user to define the integrand function f(X) byspecifying a supporting Subprogram FUNCTION F(X) and to interactively inputthe integration limits, XL and XR, along with the total number of points, N, todetermine the value of an integral based on Equation 8. As an example, we illustratebelow the estimation of a semi-circular area specifying XL = 0, XR = 2, and N = 21and defining the integrand function f(X) in FUNCTION F(X) with a statementF = SQRT(1.-(X–1.)*(X–1.))This statement describes that the center of the circle is located at (1,0) and radiusis equal to 1. When N is increased from 21 to 101 with an increment of 20, thefollowing table shows that the accuracy of trapezoidal rule is steadily increased whenthe estimated value of the semi-circle is approaching the exact value of /2.NAreaError211.5521.21%411.5640.45%611.5670.26%811.5680.19%1011.5690.13%SIMPSON’S RULEThe Trapezoidal Rule approximates the integrand function f(X) in Equation 8by a series of connected straight-line segments as illustrated in Figure 2.