162
6 Computing Integrals and Testing Code
For each of these one-dimensional integrals we apply the midpoint rule:
p(x, y) =
f
e
g(x, y, z)dz ≈
n z −1
k=0
g(x, y, z k ),
q(x) =
d
c
p(x, y)dy ≈
n y −1
j =0
p(x, y j ),
b
a
d
c
f
e
g(x, y, z)dzdydx =
b
a
q(x)dx ≈
n x −1
i=0
q(x i ),
where
z k = e +
1
2
h z + kh z , y j = c +
1
2
h y + j h y x i = a +
1
2
h x + ih x .
Starting with the formula for
b
a
d
c
f
e g(x, y, z)dzdydx and inserting the two
previous formulas gives
b
a
d
c
f
e
g(x, y, z) dzdydx ≈
h x h y h z
n x −1
i=0
n y −1
j =0
n z −1
k=0
g(a +
1
2
h x + ih x , c +
1
2
h y + j h y , e +
1
2
h z + kh z ) .
(6.30)
Note that we may apply the ideas under Direct derivation at the end of Sect. 6.7.1
to arrive at (6.30) directly: divide the domain into n x × n y × n z cells of volumes
h x h y h z ; approximate g by a constant, evaluated at the midpoint (x i , y j , z k ), in each
cell; and sum the cell integrals h x h y h z g(x i , y j , z k ).
Implementation We follow the ideas for the implementations of the midpoint rule
for a double integral. The corresponding functions are shown below and found in
the file midpoint_triple.py.
def midpoint_triple1(g, a, b, c, d, e, f, nx, ny, nz):
hx = (b - a)/nx
hy = (d - c)/ny
hz = (f - e)/nz
I = 0
for i in range(nx):
for j in range(ny):
for k in range(nz):
xi = a + hx/2 + i*hx
yj = c + hy/2 + j*hy
zk = e + hz/2 + k*hz
I = I + hx*hy*hz*g(xi, yj, zk)
return I