Page - 80 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

80 3 Computing Integrals If this function is stored in amodulefilemidpoint_double.py,we can compute someintegral,e.g., R3 2 R2 0 .2xCy/dxdyD9 inaninteractiveshellanddemonstrate that the functioncomputes the rightnumber: >>> from midpoint_double import midpoint_double1 >>> def f(x, y): ... return 2*x + y ... >>> midpoint_double1(f, 0, 2, 2, 3, 5, 5) 9.0 Reusing code for one-dimensional integrals It is very natural to write a two- dimensionalmidpointmethodaswedid in functionmidpoint_double1whenwe have the formula (3.25). However,we could alternatively ask,much aswe did in themathematics, canwe reuse awell-tested implementation for one-dimensional integrals tocomputedouble integrals?That is, canweuse functionmidpoint 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 fromSect.3.3.2“twice”?Theanswer isyes, ifwethinkaswedid in themathemat- ics: compute thedouble integral as amidpoint rule for integratingg.x/anddefine g.xi/ in terms of amidpoint rule overf in they 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) Theimportantadvantageof this implementation is thatwereuseawell-tested func- tion for the standard one-dimensional midpoint rule and that we apply the one- dimensional ruleexactlyas in themathematics. Verification via test functions How canwe test that our functions for the dou- ble integral work? The best unit test is to find a problemwhere the numerical approximation error vanishes because then we know exactly what the numerical answershouldbe. Themidpointrule isexact for linear functions, regardlessofhow many subinterval we use. Also, any linear two-dimensional function f.x;y/ D pxCqyCrwill be integratedexactlyby the two-dimensionalmidpoint rule.We maypickf.x;y/D2xCy andcreateaproper test function thatcanautomatically verify our two alternative implementations of the two-dimensionalmidpoint rule. To compute the integral off.x;y/we take advantage of SymPy to eliminate the possibilityoferrors inhandcalculations. The test functionbecomes