Seite - 83 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

3.7 DoubleandTriple Integrals 83 Implementation Wefollowthe ideasfor the implementationsof themidpoint rule for a double integral. The corresponding functions are shownbelowand found in thefilemidpoint_triple.py. def midpoint_triple1(g, a, b, c, d, e, f, nx, ny, nz): hx = (b - a)/float(nx) hy = (d - c)/float(ny) hz = (f - e)/float(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 += hx*hy*hz*g(xi, yj, zk) return I def midpoint(f, a, b, n): h = float(b-a)/n result = 0 for i in range(n): result += f((a + h/2.0) + i*h) result *= h return result def midpoint_triple2(g, a, b, c, d, e, f, nx, ny, nz): def p(x, y): return midpoint(lambda z: g(x, y, z), e, f, nz) def q(x): return midpoint(lambda y: p(x, y), c, d, ny) return midpoint(q, a, b, nx) def test_midpoint_triple(): """Test that a linear function is integrated exactly.""" def g(x, y, z): return 2*x + y - 4*z a = 0; b = 2; c = 2; d = 3; e = -1; f = 2 import sympy x, y, z = sympy.symbols(’x y z’) I_expected = sympy.integrate( g(x, y, z), (x, a, b), (y, c, d), (z, e, f)) for nx, ny, nz in (3, 5, 2), (4, 4, 4), (5, 3, 6): I_computed1 = midpoint_triple1( g, a, b, c, d, e, f, nx, ny, nz) I_computed2 = midpoint_triple2( g, a, b, c, d, e, f, nx, ny, nz) tol = 1E-14 print I_expected, I_computed1, I_computed2 assert abs(I_computed1 - I_expected) < tol assert abs(I_computed2 - I_expected) < tol if __name__ == ’__main__’: test_midpoint_triple()