Page - 162 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

162 6 ComputingIntegralsandTestingCode Foreachof theseone-dimensional integralswe apply themidpoint rule: p(x,y)= ∫ f e g(x,y,z)dz≈ nz−1∑ k=0 g(x,y,zk), q(x)= ∫ d c p(x,y)dy≈ ny−1∑ j=0 p(x,yj), ∫ b a ∫ d c ∫ f e g(x,y,z)dzdydx= ∫ b a q(x)dx≈ nx−1∑ i=0 q(xi), where zk = e+ 1 2 hz+khz, yj = c+ 1 2 hy+jhy xi =a+ 1 2 hx+ ihx . Starting with the formula for ∫b a ∫d c ∫f e g(x,y,z)dzdydx and inserting the two previousformulasgives ∫ b a ∫ d c ∫ f e g(x,y,z)dzdydx≈ hxhyhz nx−1∑ i=0 ny−1∑ j=0 nz−1∑ k=0 g(a+ 1 2 hx+ ihx,c+ 1 2 hy+jhy,e+ 1 2 hz+khz). (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 intonx ×ny ×nz cells of volumes hxhyhz; approximategbyaconstant,evaluatedat themidpoint(xi,yj,zk), ineach cell; andsumthecell integralshxhyhzg(xi,yj,zk). Implementation Wefollowthe ideas for the implementationsof themidpoint rule for a double integral. The corresponding functions are shown below and found in thefilemidpoint_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