6.7 Double and Triple Integrals
157
Demonstrating Correct Convergence Rates Computing convergence rates requires somewhat more tedious programming than for the previous tests, but it can
be applied to more general integrands. The algorithm typically goes like
• for i = 0, 1, 2, . . . , q
– n i = 2 i+1
– Compute integral with n i intervals
– Compute the error E i
– Estimate r i from (6.28) if i > 0
The corresponding code may look like
def convergence_rates(f, F, a, b, num_experiments=14):
from math import log
from numpy import zeros
expected = F(b) - F(a)
n = zeros(num_experiments, dtype=int)
E = zeros(num_experiments)
r = zeros(num_experiments-1)
for i in range(num_experiments):
n[i] = 2**(i+1)
computed = trapezoidal(f, a, b, n[i])
E[i] = abs(expected - computed)
if i > 0:
r_im1 = -log(E[i]/E[i-1])/log(n[i]/n[i-1])
# Truncate to two decimals:
r[i-1] = float(’{:.2f}’.format(r_im1))
return r
Making a test function is a matter of choosing f, F, a, and b, and then checking
the value of r i for the largest i:
def test_trapezoidal_conv_rate():
"""Check empirical convergence rates against the expected value 2."""
from math import exp
v = lambda t: 3*(t**2)*exp(t**3)
V = lambda t: exp(t**3)
a = 1.1; b = 1.9
r = convergence_rates(v, V, a, b, 14)
print(r)
tol = 0.01
msg = str(r[-4:]) # show last 4 estimated rates
assert (abs(r[-1]) - 2) < tol, msg
Running the test shows that all r i , except the first one, equal the target limit 2
within two decimals. This observation suggests a tolerance of 10 −2 .
6.7 Double and Triple Integrals
6.7.1 The Midpoint Rule for a Double Integral
Given a double integral over a rectangular domain [a, b] × [c, d],
b
a
d
c
f (x, y)dydx,
157
Demonstrating Correct Convergence Rates Computing convergence rates requires somewhat more tedious programming than for the previous tests, but it can
be applied to more general integrands. The algorithm typically goes like
• for i = 0, 1, 2, . . . , q
– n i = 2 i+1
– Compute integral with n i intervals
– Compute the error E i
– Estimate r i from (6.28) if i > 0
The corresponding code may look like
def convergence_rates(f, F, a, b, num_experiments=14):
from math import log
from numpy import zeros
expected = F(b) - F(a)
n = zeros(num_experiments, dtype=int)
E = zeros(num_experiments)
r = zeros(num_experiments-1)
for i in range(num_experiments):
n[i] = 2**(i+1)
computed = trapezoidal(f, a, b, n[i])
E[i] = abs(expected - computed)
if i > 0:
r_im1 = -log(E[i]/E[i-1])/log(n[i]/n[i-1])
# Truncate to two decimals:
r[i-1] = float(’{:.2f}’.format(r_im1))
return r
Making a test function is a matter of choosing f, F, a, and b, and then checking
the value of r i for the largest i:
def test_trapezoidal_conv_rate():
"""Check empirical convergence rates against the expected value 2."""
from math import exp
v = lambda t: 3*(t**2)*exp(t**3)
V = lambda t: exp(t**3)
a = 1.1; b = 1.9
r = convergence_rates(v, V, a, b, 14)
print(r)
tol = 0.01
msg = str(r[-4:]) # show last 4 estimated rates
assert (abs(r[-1]) - 2) < tol, msg
Running the test shows that all r i , except the first one, equal the target limit 2
within two decimals. This observation suggests a tolerance of 10 −2 .
6.7 Double and Triple Integrals
6.7.1 The Midpoint Rule for a Double Integral
Given a double integral over a rectangular domain [a, b] × [c, d],
b
a
d
c
f (x, y)dydx,