144
6 Computing Integrals and Testing Code
6.3.2 A General Implementation
We follow the advice and lessons learned from the implementation of the trapezoidal
method. Thus, we make a module midpoint.py with a general implementation
of (6.20) and a function application, just like we did with the trapezoidal
function:
def midpoint(f, a, b, n):
h = (b-a)/n
f_sum = 0
for i in range(0, n, 1):
x = (a + h/2.0) + i*h
f_sum = f_sum + f(x)
return h*f_sum
def application():
from math import exp
v = lambda t: 3*(t**2)*exp(t**3)
n = int(input(’n: ’))
numerical = midpoint(v, 0, 1, n)
# Compare with exact result
V = lambda t: exp(t**3)
exact = V(1) - V(0)
error = abs(exact - numerical)
print(’n={:d}: {:.16f}, error: {:g}’.format(n, numerical, error))
if __name__ == ’__main__’:
application()
In midpoint, observe how the x values in the loop start out at x = a +
h
2 (when
i is 0), which is in the middle of the first rectangle. The x values then increase by
h for each iteration, meaning that we repeatedly “jump” to the midpoint of the next
rectangle as i increases. This is consistent with the formula in (6.20), as is the final
x value of x = a +
h
2 + (n − 1)h. To convince yourself that the first, intermediate
and final x values are correct, look at a case with only three rectangles, for example.
When application is called, the particular problem
1
0 3t 2 e t 3 dt is computed,
i.e., the same integral that we handled with the trapezoidal method. Running the
program with n = 4 gives the output
n=4: 1.6189751378083810, error: 0.0993067
The magnitude of this error is now about 0.1 in contrast to 0.2, which we got with
the trapezoidal rule. This is in fact not accidental: one can show mathematically that
the error of the midpoint method is a bit smaller than for the trapezoidal method.
The differences are seldom of any practical importance, and on a laptop we can
easily use n = 10 6 and get the answer with an error of about 10 −12 in a couple of
seconds.
Précédent

- 165/350

Suivant