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

160 6 ComputingIntegralsandTestingCode for i in range(0, n, 1): x = (a + h/2.0) + i*h f_sum = f_sum + f(x) return h*f_sum fromSect.6.3.2“twice”?Theanswer isyes, ifwe thinkaswedid in themathemat- ics: compute the double integral as a midpoint rule for integratingg(x) and define g(xi) in terms of a midpoint rule over f in the y coordinate. The corresponding functionhasveryshort code: def midpoint_double2(f, a, b, c, d, nx, ny): def g(x): return midpoint(lambda y: f(x, y), c, d, ny) return midpoint(g, a, b, nx) The important advantage of this implementation is that we reuse a well-tested functionfor thestandardone-dimensionalmidpoint ruleand thatwe apply the one- dimensional ruleexactlyas in themathematics. Verification via Test Functions How can we test that our functions for the double integral work? The best unit test is to find a problem where the numerical approximation error vanishes because then we know exactly what the numerical answer should be. The midpoint rule is exact for linear functions, regardless of howmanysubintervalweuse.Also,anylinear two-dimensionalfunctionf(x,y)= px+qy+ r will be integrated exactly by the two-dimensional midpoint rule. We maypickf(x,y)=2x+y andcreateaproper test function that canautomatically verify our two alternative implementations of the two-dimensional midpoint rule. To compute the integral of f(x,y) we take advantage of SymPy to eliminate the possibilityoferrors inhandcalculations.The test functionbecomes def test_midpoint_double(): """Test that a linear function is integrated exactly.""" def f(x, y): return 2*x + y a = 0; b = 2; c = 2; d = 3 import sympy x, y = sympy.symbols(’x y’) I_expected = sympy.integrate(f(x, y), (x, a, b), (y, c, d)) # Test three cases: nx < ny, nx = ny, nx > ny for nx, ny in (3, 5), (4, 4), (5, 3): I_computed1 = midpoint_double1(f, a, b, c, d, nx, ny) I_computed2 = midpoint_double2(f, a, b, c, d, nx, ny) tol = 1E-14 #print I_expected, I_computed1, I_computed2 assert abs(I_computed1 - I_expected) < tol assert abs(I_computed2 - I_expected) < tol