Просмотреть файл в отдельном окне: 56e0506bcd9d653b019236b878f303be.pdf

5

Numerical Integration

and Program Volume

5.1 INTRODUCTION

Sometimes, one cannot help wonder why appears so often in a wide range of

mathematical 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 to

obtain 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 is

commonly known that if the radius of a circle is r, the diameter is equal to 2r, the

circumference is equal to 2r, and the area is equal to r2. Hence, to calculate the

diameter we simply double the radius but to calculate the circumference and the

area of a circle is more involved. The transcendental number is the result of

calculating the circumference or area of a circle by numerical integration.

In this chapter, we discuss various methods that can be adopted for the need of

numerical integration. Before we elaborate on determining the value of , let us

describe the problem of numerical integration in general.

Consider the common need of finding the area inside a closed contour C1 such

as the one shown in Figure 1A, or the area between the outside contour C2 and the

inside contour C3 shown in Figure 1B. The latter could be a practical problem of

determining the usable land size of a surveyed lot which include a pond. To evaluate

the 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 at

its extreme left and right, denoted as PL and PR, respectively. A rectangular coordinate

system 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 contour

C1 in the first quadrant of the X-Y plane. Such an arrangement makes possible to

have 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 AB

between the bottom-branch curve PLPBPR and the X-axis, from the area AT between

the top-branch curve PLPTPR and the X-axis. Approximated evaluation of the areas

AB and AT by application of digital computer proceeds first with selection of a finite

number of points Pi from PL to PR. That is, to approximate a curve such as PLPTPR

by a series of linear segments. Let N be the number of points selected on the curve

PLPTPR, then the coordinates of a typical point are (Xi,Yi) for i ranges from 1 to N

and in particular (X1,Y1) = (XL,YL) and (XN,YN) = (XR,YR). The area between a

typical 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 the

inside contour C3.

© 2001 by CRC Press LLC

FIGURE 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 the

shaded strip shown in Figure 2. Obviously, the total area AT between the top branch

of contour C1, PLPTPR, and the X-axis is the sum of all strips under the N–1 linear

segments PiPi + 1 for i = 1,2,…,N. In other words, we may mathematically write:

N ?1

AT =

N ?1

? A = ? (Y + Y

i

i =1

i

i +1

i =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 be

assigned to PR and the last point to PL. Suppose that there are M points selected

along PRPBPL, then the coordinates of these points are (Xi,Yi) for i = 1,2,…,M and

in particular (X1,Y1) = (XR,YR) and (XM,YM) = (XL,YL). Consequently, the area AB

can be calculated, similar to Equation 2 as:

M ?1

AB =

?

i =1

© 2001 by CRC Press LLC

M ?1

Ai =

? (Y + Y

i

i =1

i +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 therefore

be calculated by adding AT and AB if the numbering of the points selected on the

contour follows a clockwise direction. Let the total number of points selected around

the contour C1 be denoted as K, then K = N + (M–2) because PR and PL are re-used

in consideration of the bottom branch. Hence, the area enclosed in C1 is:

K

A = AT + AB =

? (Y + Y

i

i =1

i +1

) (Xi+1 ? Xi ) 2

(4)

where the Nth point has coordinates (XN,YN) = (XR,YR) and the first and last points

have coordinates (X1,Y1) = (XK + 1,YK + 1) = (XL,YL). And it should be evident that

in case of a cut-out, such as the contour C2 shown in Figure 1(B), the subtraction

of the area enclosed by the cut-out can be replaced by an addition of the value of

the area when it is calculated by using Equation 4 but the numbering of the points

on contour C3 is ordered in counterclockwise sense.

5.2 PROGRAM NUINTGRA — NUMERICAL INTEGRATION BY

APPLICATION OF THE TRAPEZOIDAL AND SIMPSON RULES

Program NuIntGra is designed for the need of performing numerical integration

by use of either trapezoidal rule or Simpson’s rule. These two rules will be explained

later. First, let us discuss why we need numerical integration.

Figure 3 shows a number of commonly encountered cross-sectional shapes in

engineering and scientific applications. The interactive computer program NuIntGra

has an option of allowing keyboard input of the coordinates of the vertices of the

cross section and then carrying out the area computation of cross-sectional area

based on Equation 4.

The following gives some detailed printout of the results for the cross sections

shown in Figure 3. It is important to point out that the points on the contours describing

the 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 LLC

FIGURE 3. Commonly used cross sections in engineering and scientific applications.

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

By use of a Function subprogram F(X) which defines the upper branch of a

circle 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 screen

display of this interactive run is also listed below after the Function F(X).

QuickBASIC Version

FORTRAN Version

FUNCTION F(X)

F = SQR(1(X-1)^2)

END FUNCTION

FUNCTION F(X)

F = SQRT(1.-(X-1)**2)

RETURN

END

© 2001 by CRC Press LLC

More discussion on the accuracy of will be given later after we have introduced

both the trapezoidal and Simpson’s rules which have already been incorporated in

the program NuIntgra.

TRAPEZOIDAL RULE

Returning to Figure 2, we notice that if in approximating the curve by a series

of linear segments the points selected on the curve are equally spaced in X, Equation

2 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 ?1

A T = ?X

? (Y + Y

i

i =1

i +1

)2

(6)

Or, in a different form for easy interpretation, it may also be written as:

N ?1

?

?

A T = ?X? Y1 + 2

Yi + YN ? 2

?

?

i=2

?

(7)

All the in-between heights, Yi for i ranging from 2 to N–1 that is the next to the

last, are appearing twice because they are shared by two adjacent strips whereas the

first and last heights, Y1 and YN, only appear once in Equation 7. X in Equation 7

is the common width of all strips used in summing the area.

Equation 7 is the well known Trapezoidal Rule for numerical integration. In a

general case, it can be applied for approximate evaluation of an integral by the

formula:

© 2001 by CRC Press LLC

?

XR

XL

N ?1

?

?

f ( X)dX = ?X? f1 + 2

fi + 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 should

be 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,…,N

X1 = X L

and

XN = XR

(12,13)

Program NuIntGra allows the user to define the integrand function f(X) by

specifying a supporting Subprogram FUNCTION F(X) and to interactively input

the integration limits, XL and XR, along with the total number of points, N, to

determine the value of an integral based on Equation 8. As an example, we illustrate

below the estimation of a semi-circular area specifying XL = 0, XR = 2, and N = 21

and defining the integrand function f(X) in FUNCTION F(X) with a statement

F = SQRT(1.-(X–1.)*(X–1.))

This statement describes that the center of the circle is located at (1,0) and radius

is equal to 1. When N is increased from 21 to 101 with an increment of 20, the

following table shows that the accuracy of trapezoidal rule is steadily increased when

the estimated value of the semi-circle is approaching the exact value of /2.

N

Area

Error

21

1.552

1.21%

41

1.564

0.45%

61

1.567

0.26%

81

1.568

0.19%

101

1.569

0.13%

SIMPSON’S RULE

The Trapezoidal Rule approximates the integrand function f(X) in Equation 8

by a series of connected straight-line segments as illustrated in Figure 2. These

straight-line segments can be expressed as linear functions of X. A straight line

© 2001 by CRC Press LLC

which passes through two points (Xi,Yi) and (Xi + 1,Yi + 1) may be described by the

equation:

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, we

obtain two equations

Yi = m i X i + b i

and

Yi +1 = m i X i +1 + b i

The slope mi can be easily solved by subtracting the first equation from the

second equation to be

m 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 computer

programming, it is better to write the equation describing the straight line as:

1

Y = a 0 + a1X =

?a X

i

(17)

i

i=0

That is to replace mi and bi in Equation 14 by a1 and a0, respectively. From

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

and

A logical extension of Equation 17 is to express Y as a second-order polynomial

of X, namely:

2

Y = a 0 + a1X + a 2 X 2 =

?a X

i

i=0

© 2001 by CRC Press LLC

i

(20)

FIGURE 4. Three points are required on the f(X) in order to determine the three coefficients

a0, a1, and a2.

It is a quadratic equation describing a parabola. If we select Equation 20 to

approximate a segment of f(X) for numerical evaluation of the integral in Equation

8, three points are required on the f(X) curve as illustrated in Figure 4 in order to

determine the three coefficients a0, a1, and a2. For simplicity in derivation, let the

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

a0

2

and

2

As a0 is already determined in Equation 21, it can be eliminated from

Equations 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 LLC

The solutions of the above equations can be obtained by application of Cramer’s

Rule 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 the

increment in X, X, we are ready to substitute Equations 20, 21, 26, and 27 into

Equation 8 to compute the area A1?3 in Figure 4 under the parabola. That is:

A1?3 =

2 ?X

2

? ?

0

2

a i X i dX =

i=0

?

i=0

2

? a i i +1 ?

?? i + 1 X ?? ?X

0

After simplification, it can be shown that the area1?3 is related to the ordinates

Y1, 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 by

approximating 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 ?2

A=

?

i =1,3

( odd )

M ?2

A 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 must

be an odd integer, and the stepsize X = (XU-XL)/(M–1). These stations are divided

into groups of three stations. Equation 28 has been successively employed for

evaluating 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 ?1

M ?2

?X ?

?

A=

Y +4

Yi + 2

Yi + YM ?

3 ? 1

i =3,5

i = 2 ,4

?

?

?

?

( even )

( odd )

?

© 2001 by CRC Press LLC

?

(30)

which is the well-known Simpson’s rule. Program NuIntGra has the option of using

Simpson’s rule for numerical integration. It can be shown that when this program is

applied 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 in

M

Area

Error

21

1.564

0.45%

41

1.568

0.19%

61

1.569

0.13%

81

1.570

0.05%

101

1.570

0.05%

Presented below are the program listings of NuIntGra in both QuickBASIC

and FORTRAN versions.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

FORTRAN VERSION

© 2001 by CRC Press LLC

MATLAB APPLICATION

MATLAB has a file quad.m which can perform Simpson’s Rule. To evaluate

the area of a semi-circle by application of Simpson’s Rule using quad.m, we first

prepare 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 disk

drive 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 which

the integrand function is defined whereas the second and third arguments specify

the limits of integration. Since the center of the semi-circle is located at x = 1, the

limits of integration are x = 0 and x = 2. The display resulted from the execution of

the above MATLAB statement is:

Notice that warning messages have been printed but the numerical result is not

affected.

© 2001 by CRC Press LLC

MATHEMATICA APPLICATION

Mathematica 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.14159

5.3 PROGRAM VOLUME — NUMERICAL APPROXIMATION

OF DOUBLE INTEGRATION

Program 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 is

denoted as z(x,y). If the double integration is to be carried for the region xL?x?xU

and 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 the

sub-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 base

area. The volume is there-fore a column which rises above the base area and bounded

by the z(x,y) surface, assuming that z is always positive. Mathematically, the volume

can be expressed as:

V=?

yU

yL

?

xU

xL

z ( x , y )d

(1)

If we are interested in finding the volume of sphere of radius equal to R, the

bounds can be selected as xL = yL = 0 and xU = yU = R, and let z = (R2x2y2).5. Equation

1 can then be applied to find the one-eighth of spherical volume. In fact, the result

can be obtained analytically for this z(x,y) function. We are here, however, interested

in a computational method for the case when the integrand function z(x,y) is too

complex to allow analytical solution.

The trapezoidal rule for single integration discussed in the program NuIntgra

can be extended to double integration by observing from Figure 1 that the total

volume V can be estimated as a sum of all columns erected within the space bounded

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

and

z(x, y) = 0,

© 2001 by CRC Press LLC

FIGURE 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 use

Equation 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), the

right-hand side of Equation 2 is an imaginary number. Adding Equation 3 will

remedy this situation.

If we partition the base area into a gridwork by using uniform increments x

andy along the x- and y-directions, respectively. If there are M and N equally

spaced stations along the x- and y-direction, respectively, then the increments can

be 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 averaged

for calculation of the volume, Vij, of this column by the equation:

© 2001 by CRC Press LLC

Vi, j =

(

)

1

z +z +z +z

?x?y

4 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 j

ranging from 1 to N. Or,

V=

R

??

0

=

R

0

?x?y

4

(x

2

)

.5

+ y 2 dxdy

??(

i

z i, j + z i, j+1 + z i +1, j + z i +1, j+1

j

)

(8)

The two summations in Equation 8 are loosely stated. Actually, the heights

calculated at all MxN grid-points on the base area used in Equation 7 can be separated

into 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 base

area, excluding those at the corners, are needed twice because they are shared by

two adjacent columns, and (3) all heights at interior grid-points are needed four

times in Equation 8 because they are shared by four adjacent columns. That is to

say, in terms of the subscripts I and j the weighting coefficients, wi,j, for zi,j can be

summarized as follows:

wi,j = 1

= 4

= 2

for

for

for

(i,j) = (1,1),(1,N),(M,1),(M,N),

i = 2,3,…,M-1 and j = 2,3,…,N-1

other i and j combinations

Subsequently, Equation 8 can be written as:

V=

?x?y

4

M

N

i =1

j=1

??w

z

i, j i, j

(9)

A more precise way to express V in terms zi,j is to introduce a weighting

coefficient vector for Trapezoidal rule, {wt}. Since we have averaged the four heights

of each contributing column, that is linearly connecting the four heights. That is,

the trapezoidal rule has been applied twice, one in x-direction and another in ydirection. When M and N stations are employed in x- and y-directions, respectively.

we may therefore define two weighting coefficient vectors

© 2001 by CRC Press LLC

{w } = [1

2 2 … 2 2 2]

T

(10)

{w } = [1

2 2 … 2 2 1]

T

(11)

t x

and

t y

It should be noted that the subscripts x and y are added to indicate their association with the x- and y-axes, respectively, and that the orders of these two vectors

are M and N, respectively, and that the beginning and ending components in both

vectors are equal to one and the other components are equal to two. Having defined

{wt}x and {wt}y, it is now easy to show that Equation 9 can be simply written as:

V=

T

?x

{wt }x [Z]{wt }y ?2y

2

(12)

where [Z] is a matrix of order M by N having zi,j as its elements. Since {wt}x is a

vector of order M by one, its transpose is of order one by M and {wt}y is of order

N by one, the matrix multiplication of the three matrices can be carried out and the

result does agree with the requirement on wi,j spelled out in Equation 9.

Use of weighting coefficient vectors has the advantage of extending the numerical evaluation of double integrals from trapezoidal rule to Simpson’s rule where

three adjacent heights are parabolically fitted (referring to program NuIntgra for

more details). To illustrate this point, let us first introduce a weighting coefficient

vector for Simpson’s rule as:

{w } = [1

s

4 2 … repeat of 4 and 2 … 4 1]

T

(13)

If we wish to integrate by application of Simpson’s rule in both x- and ydirections and using M and N (both must be odd) stations, the formula for the volume

is simply:

V=

T

?x

{ws}x [Z]{ws}y ?3y

3

(14)

If for some reason one wants to integrate using Simpson’s rule along x-direction

by adopting M (odd) stations, and using trapezoidal rule along y-direction by adopting N (no restriction whether odd or even) stations, then:

V=

T

?x

{ws}x [Z]{wt }y ?2y

3

(15)

Let us present a numerical example to further clarify the above concept of

numerical volume integration. Consider the problem of estimating the volume

© 2001 by CRC Press LLC

between the surface z(x,y) = 2x + 3y2 + 4 and the plane z = 0 for 0?x?2 and 1?y?2

by application of trapezoidal rule along the x direction using an increment of 0.5,

and Simpson’s rule along y direction using also an increment of 0.5. The increments

of o.5 in both x- and y-directions make M = 5 and N = 3. First, we calculate the

elements of [Z] which is of order 5 by 3:

Next, the volume is calculated to be:

Program Volume has been developed for interactive specification of the integrand function z(x,y), the integration limits xL, xU, yL, and yU, the method(s) of

integration (i.e., , trapezoidal or Simpson’s rule) and number of stations in both xand y-direction. The integrand function z(x,y) needs to be individually compiled.

Both QuickBASIC and FORTRAN versions are made available. Listings are provided below along with sample applications.

FORTRAN VERSION

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Sample Application

The FUNCTION Z(X,Y) listed above is for finding the volume under the surface

z(x,y) = (x2 + y2–4).5 over the region 0?x?2 and 0?y?2. The exact solution is

volume = 4.1889. For a sample run of the program Volume using trapezoidal rule

and 21 stations along both x- and y-directions, the screen display of interactive

communication through keyboard input and the calculated result is:

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Applications

The same calculation of one-eighth of a sphere of radius equal to 2 as in the

FORTRAN version is run but here using Simpson’s rule. The screen display is:

We have presented earlier the manual calculation of the double integration for

z(x,y) = 2x + 3y2 + 4, program Volume can be applied to have the results displayed

on the monitor screen as below. The answer is exactly the same as from manual

computation.

MATLAB APPLICATION

A Volume.m file can be created and added to MATLAB m files to calculate a

double integral when the integrand function is specified by another m file. For

integrating a function Z(X,Y) over the region X1?X?X2 and Y1?Y?Y2 by either

Trapezoidal or Simpson’s rules (designated as rule 1 or 2) and with NX and NY

stations along the X and Y directions, respectively, this file may be written as

followed:

© 2001 by CRC Press LLC

For each problem, the integrand function Z(X,Y) needs to be prepared as a m

file. In case that a hemisphere of radius 2 and centered at X = 0 and Y = 0, we may

write:

In case of Z(X,Y) = 2X + 3Y2 + 4, we may write a new file as:

Once the files Volume.m, FuncZ.m, and FuncZnew.m, the following MATLAB executions can be achieved:

© 2001 by CRC Press LLC

Notice the first and second integrations of the hemisphere use Trapezoidal and

Simpson’s rule in both X and Y, respectively. Both use 21 stations in X and Y

directions. The third integration of Z = 2X + 3Y2 + 4 over the region 0?X?2 and

1?Y?2 is carried out using Trapezoidal rule along X direction with 5 stations and

Simpson’s rule along Y direction with 3 stations.

MATLAB has a mesh plot capability of generating three-dimensional hiddenline surface. For example, when the function FuncZ is used to generate a hemispherical surface of radius equal to 2 described by a square matrix [Z], a plot shown in

Figure 6 can be obtained by entering MATLAB commands as follows:

We observe from Figure 6 that the hidden-line feature is apparent but the hemisphere appears like a semiellipsoid. This is due to the aspect ratios of the display

monitor and/or of the printer. mesh is the option of specifying different scale factors

for the X-, Y-, and Z-axes. To make the three-dimensional surface to appear as a

perfect hemisphere, user has to experiment with different scale factors for the three

axes. This is left as a homework problem. Also, mesh has option for displaying the

surface by viewing it from different angles, user is again urged to try generation of

different 3D hidden-line views.

MATHEMATICA APPLICATIONS

Mathematica has a three-dimensional plot function called Plot3d which can be

applied for drawing the hemispherical surface. Figure 7 is the result of entering the

statement:

© 2001 by CRC Press LLC

FIGURE 6.

FIGURE 7.

© 2001 by CRC Press LLC

Input]: = Sphere = Plot3D[If[4X^2Y^2>0, Sqrt[4X^2Y^2],0,

{X,–2,2},{Y,–2,2},PlotPoints->{60,60}]

The If command tests the first expression inside the brackets, it the condition

is true then the statement which follows is implemented and other the last statement

inside the bracket is implemented. In this case, the surface only rises over the base

circle of radius equal to 2. The PlotPoints command specifies how many gird points

along X- and Y-directions should be taken to plot the surface. The default number

of point is 15 in both directions. The greater the number of grid points, the smoother

the surface looks.

The same result can be obtained by first defining a surface function, say sf, and

then apply Plot3d for drawing the surface using sf as follows:

Input]: = sf[X_,Y_] = If[4X^2Y^2>0, Sqrt[4X^2Y^2], 0]

Input[2]: = Plot3D[sf[X,Y],{X,–2,2},{Y,–2,2},PlotPoints->{60,60}]

5.4 PROBLEMS

NUINTGRA

1. Having learned how to apply Trapezoidal Rule for numerical integration,

how would you find the area under the line y(x) = 1 + 2x and between

x = 1 and x = 2? Do it not by direct integration, but numerically. What

should be the stepsize for x in order to ensure an accurate result?

2. Having learned how to apply Simpson’s Rule for numerical integration,

how would you find the area under the parabolic curve y(x) = 1 + 2x +

3x2 and between x = 1 and x = 2? Do it not by direct integration but

numerically! What should be the stepsize for x in order to ensure an

accurate result?

3. If Trapezoidal Rule, instead of Simpson’s Rule, is applied for Problem 2,

find out how small should be the stepsize for x in order to achieve the

same result accurate to the fifth significant digit.

4. Could Simpson’s Rule be applied for Problem 1? Would the result be

different? If the result is the same, explain why.

5. Given five points (1,1), (2,3), (3,2), (4,5), and (5,4), use a stepsize of x =

1 to compute ydx by application of Simpson’s and Trapezoidal rules.

6. Use the trapezoidal and Simpson’s rules to find the area within the ellipse

described by the equation (x/a)2 + (y/b)2 = 1. Compare the numerical

results with the exact solution of ab.

7. Implement the integration of the function f(x) = 3e–2xsinx over the interval

from x = 0 to x = 1 (in radian) by applying both the Trapezoidal and

Simpson’s rules and using an increment of x = 0.25.

8. Find the exact solution of Problem 7 by referring to an integration formula

for f(x) from any calculus book. Decrease the increment of x (i.e., ,

© 2001 by CRC Press LLC

9.

10.

11.

12.

13.

14.

15.

16.

17.

18.

increase the number of points at which the integrand function is computed)

to try to achieve this analytical result using both Trapezoidal and Simpson’s Rules.

Apply the function Quad.m of MATLAB to solve Problem 1.

Apply the function Quad.m of MATLAB to solve Problem 2.

Apply the function Quad.m of MATLAB to solve Problem 6.

Apply the function Quad.m of MATLAB to solve Problem 7.

Apply MATLAB to spline curve-fit the five points given in Problem 5

and then integrate.

Apply the function NIntegrate of Mathematica to solve Problem 1.

Apply the function NIntegrate of Mathematica to solve Problem 2.

Apply the function NIntegrate of Mathematica to solve Problem 6.

Apply the function NIntegrate of Mathematica to solve Problem 7.

Problem 13 but apply Mathematica instead.

VOLUME

1. Apply trapezoidal rule for integration along the x direction and Simpson’s

rule along the y direction to calculate the volume under the surface z(x,y) =

3x + 2y2 + 1 over the rectangular region 0?x?2 and 0?y?4 using increments x = y = 1.

2. Rework Problem 1 except trapezoidal rule is applied for both x and y

directions.

3. Find by numerical integration of the ellipsoidal volume based on the

double integral 3[1–(x/5)2–(y/4)2]1/2dxdy and for x values ranging from

2 to 4 and y values ranging from 1 to 2. Three stations (for using Simpson’s

rule) for the x integration and two stations (for using trapezoidal rule) for

the y integration are to be adopted.

4. Find the volume between the z = 0 plane and the spherical surface z(x,y) =

[4x2 – y2]1/2 for x and y both ranging from 0 to 2 by applying the Simpson’s

rule for both x and y integrations. Three stations are to be taken along the

x direction and five stations along the y direction for the specified numerical integration.

5. Specify a FUNCTION Z(x,y) for program Volume so that the volume

enclosed by the ellipsoid (x/a)2 + (y/b)2 + (z/c)2 = 1 can be estimated by

numerical integration and compare to the exact solution of 4abc/3.2

6. In Figure 8, the shape and dimensions of a pyramid are described by the

coordinates of the five points (Xi,Yi) for I = 1,2,…,5. For application of

numerical integration to determine its volume by either trapezoidal or

Simpson’s rule, we have to partition the projected plane P2P3P4P5 into a

gridwork. At each interception point of the gridwork, (X,Y), the height

Z(X,Y) needs to be calculated which requires knowing the equations

describing the planes P1P2P3, P1P3P4, P1P4P5and P1P5P2. The equation of

a plane can be written in the form of 2(X–a) + m(Y–b) + n(Z–c) = 1

where (a,b,c) is a point on the plane and (2,m,n) are the directional cosines

© 2001 by CRC Press LLC

FIGURE 8. Problem 6.

7.

8.

9.

10.

of the unit normal vector of the plane.3 Apply the equation of plane and

assign proper values for the coordinates (Xi,Yi) describing the pyramid,

and then proceed to write a FUNCTION Z(X,Y) to determine its volume

by using program Volume.

Find the volume under the surface z = 3x2–4y + 15 over the base area of

0?x?2 and 1?y?2 by applying Simpson’s Rule along the x-direction using

an increment of x = 1, and Trapezoidal Rule along the y-direction using

an increment of y = 0.25.

How do you find the volume under the plane z = 2x–0.5y and above the

rectangular area bounded by x = 0, x = 1, y = 0, and y = 2 numerically

and not by actually integrating the z function? Explain which method and

stepsizes in x and y directions you will use, give the numerical result and

discuss how accurate it is.

Use the function FuncZnew which defines the equation Z = 2X + 3Y2 +

4 and plot the Z surface for 0?X?2 and 1?Y?2 by applying mesh of

MATLAB. Experiment with different increments of X and Y.

Modify the use of mesh by defining a vector {S} = [SX SY SZ} containing

the values of scaling factors for the three coordinate axes and then enter

mesh(Z,S) to try to improve the appearance of a hemisphere, better than

the one shown in Figure 2. Referring to Figure 2, the lowest point is the

original and the X-axis is directed to the right (width), Y-axis is directed

to the left (depth), and Z-axis is pointing upward (height). Since the

hemisphere has a radius equal to 2 and by actually measuring the width,

depth, and height to be in the approximate ratios of 2 7/8”: 2 7/16”: 2

3/4”. Based on these values, slowly adjust the values for SX, SY and SZ.

© 2001 by CRC Press LLC

FIGURE 9. Problem 11.

11. Figure 9 is obtained by using mesh to plot the surface Z = 1.5Re–2R and

R = (X2 + Y2)H for –15?X,Y?15 with increment of 1 in both X and Y

directions. Try to generate this surface by interactively entering MATLAB

commands. Apply the m file volume and modify the function FuncZnew

to accommodate this new integrand function to calculate the volume of

this surface above the 30x30 base area.

12. Apply Mathematica to solve Problem 6.

13. Apply Mathematica to solve Problem 7.

14. Apply Mathematica to solve Problem 9.

15. Apply Mathematica to solve Problem 11.

5.5 REFERENCES

1. M. Abramowitz and I. A. Stegum, editors, Handbook of Mathematical Functions with

Formulas, Graphs and Mathematical Tables, National Bureau of Standards Applied

Mathematics Series 55, Washington, DC, 1964.

2. R. C. Weast, Standard Mathematical Tables, the Chemical Rubber Co. (now CRC

Press LLC), Cleveland, OH, 13th edition, 1964.)

3. H. Flanders, R. R. Korfhage, and J. J. Price, A First Course in Calculus with Analytic

Geometry, Academic Press, New York, 1973.

© 2001 by CRC Press LLC

Numerical Integration

and Program Volume

5.1 INTRODUCTION

Sometimes, one cannot help wonder why appears so often in a wide range of

mathematical 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 to

obtain 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 is

commonly known that if the radius of a circle is r, the diameter is equal to 2r, the

circumference is equal to 2r, and the area is equal to r2. Hence, to calculate the

diameter we simply double the radius but to calculate the circumference and the

area of a circle is more involved. The transcendental number is the result of

calculating the circumference or area of a circle by numerical integration.

In this chapter, we discuss various methods that can be adopted for the need of

numerical integration. Before we elaborate on determining the value of , let us

describe the problem of numerical integration in general.

Consider the common need of finding the area inside a closed contour C1 such

as the one shown in Figure 1A, or the area between the outside contour C2 and the

inside contour C3 shown in Figure 1B. The latter could be a practical problem of

determining the usable land size of a surveyed lot which include a pond. To evaluate

the 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 at

its extreme left and right, denoted as PL and PR, respectively. A rectangular coordinate

system 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 contour

C1 in the first quadrant of the X-Y plane. Such an arrangement makes possible to

have 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 AB

between the bottom-branch curve PLPBPR and the X-axis, from the area AT between

the top-branch curve PLPTPR and the X-axis. Approximated evaluation of the areas

AB and AT by application of digital computer proceeds first with selection of a finite

number of points Pi from PL to PR. That is, to approximate a curve such as PLPTPR

by a series of linear segments. Let N be the number of points selected on the curve

PLPTPR, then the coordinates of a typical point are (Xi,Yi) for i ranges from 1 to N

and in particular (X1,Y1) = (XL,YL) and (XN,YN) = (XR,YR). The area between a

typical 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 the

inside contour C3.

© 2001 by CRC Press LLC

FIGURE 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 the

shaded strip shown in Figure 2. Obviously, the total area AT between the top branch

of contour C1, PLPTPR, and the X-axis is the sum of all strips under the N–1 linear

segments PiPi + 1 for i = 1,2,…,N. In other words, we may mathematically write:

N ?1

AT =

N ?1

? A = ? (Y + Y

i

i =1

i

i +1

i =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 be

assigned to PR and the last point to PL. Suppose that there are M points selected

along PRPBPL, then the coordinates of these points are (Xi,Yi) for i = 1,2,…,M and

in particular (X1,Y1) = (XR,YR) and (XM,YM) = (XL,YL). Consequently, the area AB

can be calculated, similar to Equation 2 as:

M ?1

AB =

?

i =1

© 2001 by CRC Press LLC

M ?1

Ai =

? (Y + Y

i

i =1

i +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 therefore

be calculated by adding AT and AB if the numbering of the points selected on the

contour follows a clockwise direction. Let the total number of points selected around

the contour C1 be denoted as K, then K = N + (M–2) because PR and PL are re-used

in consideration of the bottom branch. Hence, the area enclosed in C1 is:

K

A = AT + AB =

? (Y + Y

i

i =1

i +1

) (Xi+1 ? Xi ) 2

(4)

where the Nth point has coordinates (XN,YN) = (XR,YR) and the first and last points

have coordinates (X1,Y1) = (XK + 1,YK + 1) = (XL,YL). And it should be evident that

in case of a cut-out, such as the contour C2 shown in Figure 1(B), the subtraction

of the area enclosed by the cut-out can be replaced by an addition of the value of

the area when it is calculated by using Equation 4 but the numbering of the points

on contour C3 is ordered in counterclockwise sense.

5.2 PROGRAM NUINTGRA — NUMERICAL INTEGRATION BY

APPLICATION OF THE TRAPEZOIDAL AND SIMPSON RULES

Program NuIntGra is designed for the need of performing numerical integration

by use of either trapezoidal rule or Simpson’s rule. These two rules will be explained

later. First, let us discuss why we need numerical integration.

Figure 3 shows a number of commonly encountered cross-sectional shapes in

engineering and scientific applications. The interactive computer program NuIntGra

has an option of allowing keyboard input of the coordinates of the vertices of the

cross section and then carrying out the area computation of cross-sectional area

based on Equation 4.

The following gives some detailed printout of the results for the cross sections

shown in Figure 3. It is important to point out that the points on the contours describing

the 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 LLC

FIGURE 3. Commonly used cross sections in engineering and scientific applications.

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

By use of a Function subprogram F(X) which defines the upper branch of a

circle 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 screen

display of this interactive run is also listed below after the Function F(X).

QuickBASIC Version

FORTRAN Version

FUNCTION F(X)

F = SQR(1(X-1)^2)

END FUNCTION

FUNCTION F(X)

F = SQRT(1.-(X-1)**2)

RETURN

END

© 2001 by CRC Press LLC

More discussion on the accuracy of will be given later after we have introduced

both the trapezoidal and Simpson’s rules which have already been incorporated in

the program NuIntgra.

TRAPEZOIDAL RULE

Returning to Figure 2, we notice that if in approximating the curve by a series

of linear segments the points selected on the curve are equally spaced in X, Equation

2 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 ?1

A T = ?X

? (Y + Y

i

i =1

i +1

)2

(6)

Or, in a different form for easy interpretation, it may also be written as:

N ?1

?

?

A T = ?X? Y1 + 2

Yi + YN ? 2

?

?

i=2

?

(7)

All the in-between heights, Yi for i ranging from 2 to N–1 that is the next to the

last, are appearing twice because they are shared by two adjacent strips whereas the

first and last heights, Y1 and YN, only appear once in Equation 7. X in Equation 7

is the common width of all strips used in summing the area.

Equation 7 is the well known Trapezoidal Rule for numerical integration. In a

general case, it can be applied for approximate evaluation of an integral by the

formula:

© 2001 by CRC Press LLC

?

XR

XL

N ?1

?

?

f ( X)dX = ?X? f1 + 2

fi + 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 should

be 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,…,N

X1 = X L

and

XN = XR

(12,13)

Program NuIntGra allows the user to define the integrand function f(X) by

specifying a supporting Subprogram FUNCTION F(X) and to interactively input

the integration limits, XL and XR, along with the total number of points, N, to

determine the value of an integral based on Equation 8. As an example, we illustrate

below the estimation of a semi-circular area specifying XL = 0, XR = 2, and N = 21

and defining the integrand function f(X) in FUNCTION F(X) with a statement

F = SQRT(1.-(X–1.)*(X–1.))

This statement describes that the center of the circle is located at (1,0) and radius

is equal to 1. When N is increased from 21 to 101 with an increment of 20, the

following table shows that the accuracy of trapezoidal rule is steadily increased when

the estimated value of the semi-circle is approaching the exact value of /2.

N

Area

Error

21

1.552

1.21%

41

1.564

0.45%

61

1.567

0.26%

81

1.568

0.19%

101

1.569

0.13%

SIMPSON’S RULE

The Trapezoidal Rule approximates the integrand function f(X) in Equation 8

by a series of connected straight-line segments as illustrated in Figure 2. These

straight-line segments can be expressed as linear functions of X. A straight line

© 2001 by CRC Press LLC

which passes through two points (Xi,Yi) and (Xi + 1,Yi + 1) may be described by the

equation:

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, we

obtain two equations

Yi = m i X i + b i

and

Yi +1 = m i X i +1 + b i

The slope mi can be easily solved by subtracting the first equation from the

second equation to be

m 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 computer

programming, it is better to write the equation describing the straight line as:

1

Y = a 0 + a1X =

?a X

i

(17)

i

i=0

That is to replace mi and bi in Equation 14 by a1 and a0, respectively. From

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

and

A logical extension of Equation 17 is to express Y as a second-order polynomial

of X, namely:

2

Y = a 0 + a1X + a 2 X 2 =

?a X

i

i=0

© 2001 by CRC Press LLC

i

(20)

FIGURE 4. Three points are required on the f(X) in order to determine the three coefficients

a0, a1, and a2.

It is a quadratic equation describing a parabola. If we select Equation 20 to

approximate a segment of f(X) for numerical evaluation of the integral in Equation

8, three points are required on the f(X) curve as illustrated in Figure 4 in order to

determine the three coefficients a0, a1, and a2. For simplicity in derivation, let the

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

a0

2

and

2

As a0 is already determined in Equation 21, it can be eliminated from

Equations 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 LLC

The solutions of the above equations can be obtained by application of Cramer’s

Rule 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 the

increment in X, X, we are ready to substitute Equations 20, 21, 26, and 27 into

Equation 8 to compute the area A1?3 in Figure 4 under the parabola. That is:

A1?3 =

2 ?X

2

? ?

0

2

a i X i dX =

i=0

?

i=0

2

? a i i +1 ?

?? i + 1 X ?? ?X

0

After simplification, it can be shown that the area1?3 is related to the ordinates

Y1, 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 by

approximating 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 ?2

A=

?

i =1,3

( odd )

M ?2

A 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 must

be an odd integer, and the stepsize X = (XU-XL)/(M–1). These stations are divided

into groups of three stations. Equation 28 has been successively employed for

evaluating 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 ?1

M ?2

?X ?

?

A=

Y +4

Yi + 2

Yi + YM ?

3 ? 1

i =3,5

i = 2 ,4

?

?

?

?

( even )

( odd )

?

© 2001 by CRC Press LLC

?

(30)

which is the well-known Simpson’s rule. Program NuIntGra has the option of using

Simpson’s rule for numerical integration. It can be shown that when this program is

applied 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 in

M

Area

Error

21

1.564

0.45%

41

1.568

0.19%

61

1.569

0.13%

81

1.570

0.05%

101

1.570

0.05%

Presented below are the program listings of NuIntGra in both QuickBASIC

and FORTRAN versions.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

FORTRAN VERSION

© 2001 by CRC Press LLC

MATLAB APPLICATION

MATLAB has a file quad.m which can perform Simpson’s Rule. To evaluate

the area of a semi-circle by application of Simpson’s Rule using quad.m, we first

prepare 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 disk

drive 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 which

the integrand function is defined whereas the second and third arguments specify

the limits of integration. Since the center of the semi-circle is located at x = 1, the

limits of integration are x = 0 and x = 2. The display resulted from the execution of

the above MATLAB statement is:

Notice that warning messages have been printed but the numerical result is not

affected.

© 2001 by CRC Press LLC

MATHEMATICA APPLICATION

Mathematica 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.14159

5.3 PROGRAM VOLUME — NUMERICAL APPROXIMATION

OF DOUBLE INTEGRATION

Program 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 is

denoted as z(x,y). If the double integration is to be carried for the region xL?x?xU

and 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 the

sub-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 base

area. The volume is there-fore a column which rises above the base area and bounded

by the z(x,y) surface, assuming that z is always positive. Mathematically, the volume

can be expressed as:

V=?

yU

yL

?

xU

xL

z ( x , y )d

(1)

If we are interested in finding the volume of sphere of radius equal to R, the

bounds can be selected as xL = yL = 0 and xU = yU = R, and let z = (R2x2y2).5. Equation

1 can then be applied to find the one-eighth of spherical volume. In fact, the result

can be obtained analytically for this z(x,y) function. We are here, however, interested

in a computational method for the case when the integrand function z(x,y) is too

complex to allow analytical solution.

The trapezoidal rule for single integration discussed in the program NuIntgra

can be extended to double integration by observing from Figure 1 that the total

volume V can be estimated as a sum of all columns erected within the space bounded

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

and

z(x, y) = 0,

© 2001 by CRC Press LLC

FIGURE 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 use

Equation 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), the

right-hand side of Equation 2 is an imaginary number. Adding Equation 3 will

remedy this situation.

If we partition the base area into a gridwork by using uniform increments x

andy along the x- and y-directions, respectively. If there are M and N equally

spaced stations along the x- and y-direction, respectively, then the increments can

be 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 averaged

for calculation of the volume, Vij, of this column by the equation:

© 2001 by CRC Press LLC

Vi, j =

(

)

1

z +z +z +z

?x?y

4 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 j

ranging from 1 to N. Or,

V=

R

??

0

=

R

0

?x?y

4

(x

2

)

.5

+ y 2 dxdy

??(

i

z i, j + z i, j+1 + z i +1, j + z i +1, j+1

j

)

(8)

The two summations in Equation 8 are loosely stated. Actually, the heights

calculated at all MxN grid-points on the base area used in Equation 7 can be separated

into 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 base

area, excluding those at the corners, are needed twice because they are shared by

two adjacent columns, and (3) all heights at interior grid-points are needed four

times in Equation 8 because they are shared by four adjacent columns. That is to

say, in terms of the subscripts I and j the weighting coefficients, wi,j, for zi,j can be

summarized as follows:

wi,j = 1

= 4

= 2

for

for

for

(i,j) = (1,1),(1,N),(M,1),(M,N),

i = 2,3,…,M-1 and j = 2,3,…,N-1

other i and j combinations

Subsequently, Equation 8 can be written as:

V=

?x?y

4

M

N

i =1

j=1

??w

z

i, j i, j

(9)

A more precise way to express V in terms zi,j is to introduce a weighting

coefficient vector for Trapezoidal rule, {wt}. Since we have averaged the four heights

of each contributing column, that is linearly connecting the four heights. That is,

the trapezoidal rule has been applied twice, one in x-direction and another in ydirection. When M and N stations are employed in x- and y-directions, respectively.

we may therefore define two weighting coefficient vectors

© 2001 by CRC Press LLC

{w } = [1

2 2 … 2 2 2]

T

(10)

{w } = [1

2 2 … 2 2 1]

T

(11)

t x

and

t y

It should be noted that the subscripts x and y are added to indicate their association with the x- and y-axes, respectively, and that the orders of these two vectors

are M and N, respectively, and that the beginning and ending components in both

vectors are equal to one and the other components are equal to two. Having defined

{wt}x and {wt}y, it is now easy to show that Equation 9 can be simply written as:

V=

T

?x

{wt }x [Z]{wt }y ?2y

2

(12)

where [Z] is a matrix of order M by N having zi,j as its elements. Since {wt}x is a

vector of order M by one, its transpose is of order one by M and {wt}y is of order

N by one, the matrix multiplication of the three matrices can be carried out and the

result does agree with the requirement on wi,j spelled out in Equation 9.

Use of weighting coefficient vectors has the advantage of extending the numerical evaluation of double integrals from trapezoidal rule to Simpson’s rule where

three adjacent heights are parabolically fitted (referring to program NuIntgra for

more details). To illustrate this point, let us first introduce a weighting coefficient

vector for Simpson’s rule as:

{w } = [1

s

4 2 … repeat of 4 and 2 … 4 1]

T

(13)

If we wish to integrate by application of Simpson’s rule in both x- and ydirections and using M and N (both must be odd) stations, the formula for the volume

is simply:

V=

T

?x

{ws}x [Z]{ws}y ?3y

3

(14)

If for some reason one wants to integrate using Simpson’s rule along x-direction

by adopting M (odd) stations, and using trapezoidal rule along y-direction by adopting N (no restriction whether odd or even) stations, then:

V=

T

?x

{ws}x [Z]{wt }y ?2y

3

(15)

Let us present a numerical example to further clarify the above concept of

numerical volume integration. Consider the problem of estimating the volume

© 2001 by CRC Press LLC

between the surface z(x,y) = 2x + 3y2 + 4 and the plane z = 0 for 0?x?2 and 1?y?2

by application of trapezoidal rule along the x direction using an increment of 0.5,

and Simpson’s rule along y direction using also an increment of 0.5. The increments

of o.5 in both x- and y-directions make M = 5 and N = 3. First, we calculate the

elements of [Z] which is of order 5 by 3:

Next, the volume is calculated to be:

Program Volume has been developed for interactive specification of the integrand function z(x,y), the integration limits xL, xU, yL, and yU, the method(s) of

integration (i.e., , trapezoidal or Simpson’s rule) and number of stations in both xand y-direction. The integrand function z(x,y) needs to be individually compiled.

Both QuickBASIC and FORTRAN versions are made available. Listings are provided below along with sample applications.

FORTRAN VERSION

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Sample Application

The FUNCTION Z(X,Y) listed above is for finding the volume under the surface

z(x,y) = (x2 + y2–4).5 over the region 0?x?2 and 0?y?2. The exact solution is

volume = 4.1889. For a sample run of the program Volume using trapezoidal rule

and 21 stations along both x- and y-directions, the screen display of interactive

communication through keyboard input and the calculated result is:

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Applications

The same calculation of one-eighth of a sphere of radius equal to 2 as in the

FORTRAN version is run but here using Simpson’s rule. The screen display is:

We have presented earlier the manual calculation of the double integration for

z(x,y) = 2x + 3y2 + 4, program Volume can be applied to have the results displayed

on the monitor screen as below. The answer is exactly the same as from manual

computation.

MATLAB APPLICATION

A Volume.m file can be created and added to MATLAB m files to calculate a

double integral when the integrand function is specified by another m file. For

integrating a function Z(X,Y) over the region X1?X?X2 and Y1?Y?Y2 by either

Trapezoidal or Simpson’s rules (designated as rule 1 or 2) and with NX and NY

stations along the X and Y directions, respectively, this file may be written as

followed:

© 2001 by CRC Press LLC

For each problem, the integrand function Z(X,Y) needs to be prepared as a m

file. In case that a hemisphere of radius 2 and centered at X = 0 and Y = 0, we may

write:

In case of Z(X,Y) = 2X + 3Y2 + 4, we may write a new file as:

Once the files Volume.m, FuncZ.m, and FuncZnew.m, the following MATLAB executions can be achieved:

© 2001 by CRC Press LLC

Notice the first and second integrations of the hemisphere use Trapezoidal and

Simpson’s rule in both X and Y, respectively. Both use 21 stations in X and Y

directions. The third integration of Z = 2X + 3Y2 + 4 over the region 0?X?2 and

1?Y?2 is carried out using Trapezoidal rule along X direction with 5 stations and

Simpson’s rule along Y direction with 3 stations.

MATLAB has a mesh plot capability of generating three-dimensional hiddenline surface. For example, when the function FuncZ is used to generate a hemispherical surface of radius equal to 2 described by a square matrix [Z], a plot shown in

Figure 6 can be obtained by entering MATLAB commands as follows:

We observe from Figure 6 that the hidden-line feature is apparent but the hemisphere appears like a semiellipsoid. This is due to the aspect ratios of the display

monitor and/or of the printer. mesh is the option of specifying different scale factors

for the X-, Y-, and Z-axes. To make the three-dimensional surface to appear as a

perfect hemisphere, user has to experiment with different scale factors for the three

axes. This is left as a homework problem. Also, mesh has option for displaying the

surface by viewing it from different angles, user is again urged to try generation of

different 3D hidden-line views.

MATHEMATICA APPLICATIONS

Mathematica has a three-dimensional plot function called Plot3d which can be

applied for drawing the hemispherical surface. Figure 7 is the result of entering the

statement:

© 2001 by CRC Press LLC

FIGURE 6.

FIGURE 7.

© 2001 by CRC Press LLC

Input]: = Sphere = Plot3D[If[4X^2Y^2>0, Sqrt[4X^2Y^2],0,

{X,–2,2},{Y,–2,2},PlotPoints->{60,60}]

The If command tests the first expression inside the brackets, it the condition

is true then the statement which follows is implemented and other the last statement

inside the bracket is implemented. In this case, the surface only rises over the base

circle of radius equal to 2. The PlotPoints command specifies how many gird points

along X- and Y-directions should be taken to plot the surface. The default number

of point is 15 in both directions. The greater the number of grid points, the smoother

the surface looks.

The same result can be obtained by first defining a surface function, say sf, and

then apply Plot3d for drawing the surface using sf as follows:

Input]: = sf[X_,Y_] = If[4X^2Y^2>0, Sqrt[4X^2Y^2], 0]

Input[2]: = Plot3D[sf[X,Y],{X,–2,2},{Y,–2,2},PlotPoints->{60,60}]

5.4 PROBLEMS

NUINTGRA

1. Having learned how to apply Trapezoidal Rule for numerical integration,

how would you find the area under the line y(x) = 1 + 2x and between

x = 1 and x = 2? Do it not by direct integration, but numerically. What

should be the stepsize for x in order to ensure an accurate result?

2. Having learned how to apply Simpson’s Rule for numerical integration,

how would you find the area under the parabolic curve y(x) = 1 + 2x +

3x2 and between x = 1 and x = 2? Do it not by direct integration but

numerically! What should be the stepsize for x in order to ensure an

accurate result?

3. If Trapezoidal Rule, instead of Simpson’s Rule, is applied for Problem 2,

find out how small should be the stepsize for x in order to achieve the

same result accurate to the fifth significant digit.

4. Could Simpson’s Rule be applied for Problem 1? Would the result be

different? If the result is the same, explain why.

5. Given five points (1,1), (2,3), (3,2), (4,5), and (5,4), use a stepsize of x =

1 to compute ydx by application of Simpson’s and Trapezoidal rules.

6. Use the trapezoidal and Simpson’s rules to find the area within the ellipse

described by the equation (x/a)2 + (y/b)2 = 1. Compare the numerical

results with the exact solution of ab.

7. Implement the integration of the function f(x) = 3e–2xsinx over the interval

from x = 0 to x = 1 (in radian) by applying both the Trapezoidal and

Simpson’s rules and using an increment of x = 0.25.

8. Find the exact solution of Problem 7 by referring to an integration formula

for f(x) from any calculus book. Decrease the increment of x (i.e., ,

© 2001 by CRC Press LLC

9.

10.

11.

12.

13.

14.

15.

16.

17.

18.

increase the number of points at which the integrand function is computed)

to try to achieve this analytical result using both Trapezoidal and Simpson’s Rules.

Apply the function Quad.m of MATLAB to solve Problem 1.

Apply the function Quad.m of MATLAB to solve Problem 2.

Apply the function Quad.m of MATLAB to solve Problem 6.

Apply the function Quad.m of MATLAB to solve Problem 7.

Apply MATLAB to spline curve-fit the five points given in Problem 5

and then integrate.

Apply the function NIntegrate of Mathematica to solve Problem 1.

Apply the function NIntegrate of Mathematica to solve Problem 2.

Apply the function NIntegrate of Mathematica to solve Problem 6.

Apply the function NIntegrate of Mathematica to solve Problem 7.

Problem 13 but apply Mathematica instead.

VOLUME

1. Apply trapezoidal rule for integration along the x direction and Simpson’s

rule along the y direction to calculate the volume under the surface z(x,y) =

3x + 2y2 + 1 over the rectangular region 0?x?2 and 0?y?4 using increments x = y = 1.

2. Rework Problem 1 except trapezoidal rule is applied for both x and y

directions.

3. Find by numerical integration of the ellipsoidal volume based on the

double integral 3[1–(x/5)2–(y/4)2]1/2dxdy and for x values ranging from

2 to 4 and y values ranging from 1 to 2. Three stations (for using Simpson’s

rule) for the x integration and two stations (for using trapezoidal rule) for

the y integration are to be adopted.

4. Find the volume between the z = 0 plane and the spherical surface z(x,y) =

[4x2 – y2]1/2 for x and y both ranging from 0 to 2 by applying the Simpson’s

rule for both x and y integrations. Three stations are to be taken along the

x direction and five stations along the y direction for the specified numerical integration.

5. Specify a FUNCTION Z(x,y) for program Volume so that the volume

enclosed by the ellipsoid (x/a)2 + (y/b)2 + (z/c)2 = 1 can be estimated by

numerical integration and compare to the exact solution of 4abc/3.2

6. In Figure 8, the shape and dimensions of a pyramid are described by the

coordinates of the five points (Xi,Yi) for I = 1,2,…,5. For application of

numerical integration to determine its volume by either trapezoidal or

Simpson’s rule, we have to partition the projected plane P2P3P4P5 into a

gridwork. At each interception point of the gridwork, (X,Y), the height

Z(X,Y) needs to be calculated which requires knowing the equations

describing the planes P1P2P3, P1P3P4, P1P4P5and P1P5P2. The equation of

a plane can be written in the form of 2(X–a) + m(Y–b) + n(Z–c) = 1

where (a,b,c) is a point on the plane and (2,m,n) are the directional cosines

© 2001 by CRC Press LLC

FIGURE 8. Problem 6.

7.

8.

9.

10.

of the unit normal vector of the plane.3 Apply the equation of plane and

assign proper values for the coordinates (Xi,Yi) describing the pyramid,

and then proceed to write a FUNCTION Z(X,Y) to determine its volume

by using program Volume.

Find the volume under the surface z = 3x2–4y + 15 over the base area of

0?x?2 and 1?y?2 by applying Simpson’s Rule along the x-direction using

an increment of x = 1, and Trapezoidal Rule along the y-direction using

an increment of y = 0.25.

How do you find the volume under the plane z = 2x–0.5y and above the

rectangular area bounded by x = 0, x = 1, y = 0, and y = 2 numerically

and not by actually integrating the z function? Explain which method and

stepsizes in x and y directions you will use, give the numerical result and

discuss how accurate it is.

Use the function FuncZnew which defines the equation Z = 2X + 3Y2 +

4 and plot the Z surface for 0?X?2 and 1?Y?2 by applying mesh of

MATLAB. Experiment with different increments of X and Y.

Modify the use of mesh by defining a vector {S} = [SX SY SZ} containing

the values of scaling factors for the three coordinate axes and then enter

mesh(Z,S) to try to improve the appearance of a hemisphere, better than

the one shown in Figure 2. Referring to Figure 2, the lowest point is the

original and the X-axis is directed to the right (width), Y-axis is directed

to the left (depth), and Z-axis is pointing upward (height). Since the

hemisphere has a radius equal to 2 and by actually measuring the width,

depth, and height to be in the approximate ratios of 2 7/8”: 2 7/16”: 2

3/4”. Based on these values, slowly adjust the values for SX, SY and SZ.

© 2001 by CRC Press LLC

FIGURE 9. Problem 11.

11. Figure 9 is obtained by using mesh to plot the surface Z = 1.5Re–2R and

R = (X2 + Y2)H for –15?X,Y?15 with increment of 1 in both X and Y

directions. Try to generate this surface by interactively entering MATLAB

commands. Apply the m file volume and modify the function FuncZnew

to accommodate this new integrand function to calculate the volume of

this surface above the 30x30 base area.

12. Apply Mathematica to solve Problem 6.

13. Apply Mathematica to solve Problem 7.

14. Apply Mathematica to solve Problem 9.

15. Apply Mathematica to solve Problem 11.

5.5 REFERENCES

1. M. Abramowitz and I. A. Stegum, editors, Handbook of Mathematical Functions with

Formulas, Graphs and Mathematical Tables, National Bureau of Standards Applied

Mathematics Series 55, Washington, DC, 1964.

2. R. C. Weast, Standard Mathematical Tables, the Chemical Rubber Co. (now CRC

Press LLC), Cleveland, OH, 13th edition, 1964.)

3. H. Flanders, R. R. Korfhage, and J. J. Price, A First Course in Calculus with Analytic

Geometry, Academic Press, New York, 1973.

© 2001 by CRC Press LLC