160
6 Computing Integrals and Testing Code
for i in range(0, n, 1):
x = (a + h/2.0) + i*h
f_sum = f_sum + f(x)
return h*f_sum
from Sect. 6.3.2 “twice”? The answer is yes, if we think as we did in the mathematics: compute the double integral as a midpoint rule for integrating g(x) and define
g(x i ) in terms of a midpoint rule over f in the y coordinate. The corresponding
function has very short 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
function for the standard one-dimensional midpoint rule and that we apply the onedimensional rule exactly as in the mathematics.
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
how many subinterval we use. Also, any linear two-dimensional function f (x, y) =
px + qy + r will be integrated exactly by the two-dimensional midpoint rule. We
may pick f (x, y) = 2x + y and create a proper test function that can automatically
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
possibility of errors in hand calculations. The test function becomes
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
6 Computing Integrals and Testing Code
for i in range(0, n, 1):
x = (a + h/2.0) + i*h
f_sum = f_sum + f(x)
return h*f_sum
from Sect. 6.3.2 “twice”? The answer is yes, if we think as we did in the mathematics: compute the double integral as a midpoint rule for integrating g(x) and define
g(x i ) in terms of a midpoint rule over f in the y coordinate. The corresponding
function has very short 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
function for the standard one-dimensional midpoint rule and that we apply the onedimensional rule exactly as in the mathematics.
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
how many subinterval we use. Also, any linear two-dimensional function f (x, y) =
px + qy + r will be integrated exactly by the two-dimensional midpoint rule. We
may pick f (x, y) = 2x + y and create a proper test function that can automatically
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
possibility of errors in hand calculations. The test function becomes
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