CH-05 (Pao - Engineering Analysis), страница 2

Thesestraight-line segments can be expressed as linear functions of X. A straight line© 2001 by CRC Press LLCwhich passes through two points (Xi,Yi) and (Xi + 1,Yi + 1) may be described by theequation:Y = m i X + bi(14)where mi is the slope and bi is the intercept at the Y-axis of the ith straight line.Upon substituting the coordinates (Xi,Yi) and (Xi + 1,Yi + 1) into Equation 14, weobtain two equationsYi = m i X i + b iandYi +1 = m i X i +1 + b iThe slope mi can be easily solved by subtracting the first equation from thesecond equation to bem i = (Yi +1 − Yi ) ( X i +1 − X i )(15)Substituting mi into the Yi equation, we obtain the intercept to be:b i = Yi − (Yi +1 − Yi ) X i ( X i +1 − X i )(16)For the convenience of further development as well as for ease in computerprogramming, it is better to write the equation describing the straight line as:1Y = a 0 + a1X =∑a Xi(17)ii=0That is to replace mi and bi in Equation 14 by a1 and a0, respectively.

FromEquations 15 and 16, we therefore can have:a 0 = (Yi +1 − Yi ) ( X i +1 − X i )(18)a1 = Yi − (Yi +1 − Yi ) X i ( X i +1 − X i )(19)andA logical extension of Equation 17 is to express Y as a second-order polynomialof X, namely:2Y = a 0 + a1X + a 2 X 2 =∑a Xii=0© 2001 by CRC Press LLCi(20)FIGURE 4. Three points are required on the f(X) in order to determine the three coefficientsa0, a1, and a2.It is a quadratic equation describing a parabola.

If we select Equation 20 toapproximate a segment of f(X) for numerical evaluation of the integral in Equation8, three points are required on the f(X) curve as illustrated in Figure 4 in order todetermine the three coefficients a0, a1, and a2. For simplicity in derivation, let thethree points be (0,Y1), (X,Y2), and (2X,Y3). Upon substitution into Equation 20,we obtain:= Y1(21)a 0 + ( ∆X)a1 + ( ∆X) a 2= Y2(22)a 0 + (2 ∆X)a1 + (2 ∆X) a 2= Y3(23)a02and2As a0 is already determined in Equation 21, it can be eliminated fromEquations 22 and 23 to yield:a1 + ( ∆X)a 2 = (Y2 − Y1 ) ∆X(24)a1 + (2 ∆X)a 2 = (Y3 − Y1 ) 2 ∆X(25)and© 2001 by CRC Press LLCThe solutions of the above equations can be obtained by application of Cramer’sRule as:a1 = 2 (Y2 − Y1 ) ∆X − (Y3 − Y1 ) 2 ∆X = ( −3Y1 + 4Y2 − Y3 ) 2 ∆X(26)and[]a 2 = (Y2 − Y1 ) ∆X − a1 ∆X = (Y1 − 2Y2 + Y3 ) 2( ∆X)2(27)Having derived a0, a1, and a2 in terms of the ordinates Y1, Y2, and Y3, and theincrement in X, X, we are ready to substitute Equations 20, 21, 26, and 27 intoEquation 8 to compute the area A1ø3 in Figure 4 under the parabola.

That is:A1→3 =2 ∆X2∫ ∑02a i X i dX =i=0∑i=02 a i i +1  i + 1 X  ∆X0After simplification, it can be shown that the area1ø3 is related to the ordinatesY1, Y2, and Y3, and the increment X by the equation:A1→3 =∆X(Y1 + 4Y2 + Y3 )3(28)When the above-described procedure is applied for numerical integration byapproximating the curve of the integrand Y(X) as a series of connected parabolas,we can expand Equation 28 to cover the limits of integration to obtain:M −2A=∑i =1,3( odd )M −2A i→ i + 2 =∆X(Y + 4Yi+1 + Yi+2 )3 i =1,3 i∑(29)( odd )Notice that the limits of integration are treated by having M stations which mustbe an odd integer, and the stepsize X = (XU-XL)/(M–1). These stations are dividedinto groups of three stations.

Equation 28 has been successively employed forevaluating the adjacent areas A1ø3, A3ø5, …, AM–2øM in order to arrive at Equation 29.Equation 29 can also be written in the form of:M −1M −2∆X A=Y +4Yi + 2Yi + YM 3  1i =3,5i = 2 ,4( even )( odd )∑© 2001 by CRC Press LLC∑(30)which is the well-known Simpson’s rule. Program NuIntGra has the option of usingSimpson’s rule for numerical integration. It can be shown that when this program isapplied for the integration of the semi-circular area under the curve Y(X) = (1–X2).5,the Simpson’s rule using different M stations will result inMAreaError211.5640.45%411.5680.19%611.5690.13%811.5700.05%1011.5700.05%Presented below are the program listings of NuIntGra in both QuickBASICand FORTRAN versions.QUICKBASIC VERSION© 2001 by CRC Press LLCFORTRAN VERSION© 2001 by CRC Press LLCMATLAB APPLICATIONMATLAB has a file quad.m which can perform Simpson’s Rule.

To evaluatethe area of a semi-circle by application of Simpson’s Rule using quad.m, we firstprepare the integrand function as a m file as follows:If this file integrnd.m is stored on a disk which has been inserted in the diskdrive A, quad.m is to be applied as follows:>> Area = quad(‘A:integrnd’,0,2)Notice that quad has three arguments.

The first argument is the m file in whichthe integrand function is defined whereas the second and third arguments specifythe limits of integration. Since the center of the semi-circle is located at x = 1, thelimits of integration are x = 0 and x = 2. The display resulted from the execution ofthe above MATLAB statement is:Notice that warning messages have been printed but the numerical result is notaffected.© 2001 by CRC Press LLCMATHEMATICA APPLICATIONMathematica numerically integrate a function f(x) over the interval x = a to x =b by use of the function NIntegrate.

The following example demonstrates the computation of one quarter of a circle having a radius equal to 2:In[1]: = NIntegrate[Sqrt[4x^2], {x, 0, 2}]Out[1] = 3.141595.3 PROGRAM VOLUME — NUMERICAL APPROXIMATIONOF DOUBLE INTEGRATIONProgram Volume is designed for numerical calculation of double integration involving an integrand function of two variables. For convenience of graphical interpretation, the two variables x and y are usually chosen and the integrand function isdenoted as z(x,y).

If the double integration is to be carried for the region xL≤x≤xUand yL≤y≤yU, the value to be calculated is the volume bounded by the z surface, z =0 plane, and the four bounding planes x = xL, x = xU, y = yL, and y = yU where thesub-scripts L and U are used to indicate the lower and upper bounds, respectively.The rectangular region xL≤x≤xU and yL≤y≤yU on the z = 0 plane is called the basearea. The volume is there-fore a column which rises above the base area and boundedby the z(x,y) surface, assuming that z is always positive.

Mathematically, the volumecan be expressed as:V=∫yUyL∫xUxLz ( x , y )d(1)If we are interested in finding the volume of sphere of radius equal to R, thebounds can be selected as xL = yL = 0 and xU = yU = R, and let z = (R2x2y2).5. Equation1 can then be applied to find the one-eighth of spherical volume. In fact, the resultcan be obtained analytically for this z(x,y) function.

We are here, however, interestedin a computational method for the case when the integrand function z(x,y) is toocomplex to allow analytical solution.The trapezoidal rule for single integration discussed in the program NuIntgracan be extended to double integration by observing from Figure 1 that the totalvolume V can be estimated as a sum of all columns erected within the space boundedby the z surface and the base area. In Figure 5, the integrand functions used are:(z(x, y) = R 2 − x 2 − y 2)0.5, for x 2 + y 2 ≤ R 2(2)for x 2 + y 2 > R 2(3)andz(x, y) = 0,© 2001 by CRC Press LLCFIGURE 5. In this figure , the integrand functions used are Equations 2 and 3.Notice that Equation 3 is an added extension of Equation 2 because if we useEquation 1 and the upper limits are xU = yU = R, a point outside of the boundary x2+ y2 = R2 on the base area 0≤x≤R and 0≤y≤R is selected for evaluating z(x,y), theright-hand side of Equation 2 is an imaginary number.

Adding Equation 3 willremedy this situation.If we partition the base area into a gridwork by using uniform increments xandy along the x- and y-directions, respectively. If there are M and N equallyspaced stations along the x- and y-direction, respectively, then the increments canbe calculated by the equations:∆x = (x U − x L ) (M − 1) = R (M − 1)(4)∆y = (y U − y L ) ( N − 1)(5)and= R ( N − 1)At a typical grid-point on the base area, (xi,yj), there are three neighboring points(xi,yj + 1), (xi + 1, yj), and (xi + 1,yj + 1). The z values at these four points can be averagedfor calculation of the volume, Vij, of this column by the equation:© 2001 by CRC Press LLCVi, j =()1z +z +z +z∆x∆y4 i, j i, j+1 i +1, j i +1, j+1(6)where:(z i, j ≡ z x i , y j)(7)The total volume is then the sum of all Vi,j for i ranging from 1 to M and jranging from 1 to N.

Or,V=R∫∫0=R0∆x∆y4(x2).5+ y 2 dxdy∑∑(iz i, j + z i, j+1 + z i +1, j + z i +1, j+1j)(8)The two summations in Equation 8 are loosely stated. Actually, the heightscalculated at all MxN grid-points on the base area used in Equation 7 can be separatedinto three groups: (1) those heights at the corners whose coordinates are (0,0), (0,R),(R,0), and (R,R), are needed only once, (2) those heights on the edges of the basearea, excluding those at the corners, are needed twice because they are shared bytwo adjacent columns, and (3) all heights at interior grid-points are needed fourtimes in Equation 8 because they are shared by four adjacent columns.

That is tosay, in terms of the subscripts I and j the weighting coefficients, wi,j, for zi,j can besummarized as follows:wi,j = 1= 4= 2forforfor(i,j) = (1,1),(1,N),(M,1),(M,N),i = 2,3,…,M-1 and j = 2,3,…,N-1other i and j combinationsSubsequently, Equation 8 can be written as:V=∆x∆y4MNi =1j=1∑∑wzi, j i, j(9)A more precise way to express V in terms zi,j is to introduce a weightingcoefficient vector for Trapezoidal rule, {wt}. Since we have averaged the four heightsof each contributing column, that is linearly connecting the four heights.