6.2 The Composite Trapezoidal Rule
141
v_sum = 0
for i in range(1, n, 1):
t = a + i*dt
v_sum = v_sum + v(t)
numerical = dt*(0.5*v(a) + v_sum + 0.5*v(b))
V = lambda t: exp(t**3)
exact_value = V(b) - V(a)
error = abs(exact_value - numerical)
rel_error = (error/exact_value)*100
print(’n={:d}: {:.16f}, error: {:g}’.format(n, numerical, error))
Unfortunately, the two other problems (1. and 3.) remain, and they are fundamental.
Computing Another Integral Suppose you next want to compute another integral,
say
1.1
−1 e −x 2 dx, using the previous specific implementation as your “starting
point”. What changes are required in the code then?
First of all, an anti-derivative can not (easily) be found 2,3 for this new integrand,
so we drop computing the integration error, and must remove the corresponding
code lines. In addition,
• the notation should be changed to fit the new problem. Thus, t and dt should be
replaced by x and h. Also, the integrand is (most likely) not a velocity any more,
so the name v should be changed with, e.g., f. Similarly, v_sum should rather be
f_sum then.
• the formula for v (or f) must be replaced by a new formula
• the limits a and b must be changed
These changes are straightforward to implement, but they are scattered around in
the program, a fact that requires us to be very careful so we do not introduce new
programming errors while we modify the code. It is also very easy to forget one or
two of the required changes.
For the sake of comparison, we might see how easy it is to rather use our general
implementation in trapezoidal.py for the task. With the following interactive
session, it should be clear that this implementation allows us to compute the new
integral
1.1
−1 e −x 2 dx without touching the implemented mathematical algorithm! We
can simply do:
In [1]: from trapezoidal import trapezoidal # ...general implementation
In [2]: from math import exp
In [3]: trapezoidal(lambda x: exp(-x**2), -1, 1.1, 400)
Out[3]: 1.5268823686123285
2 You cannot integrate e −x 2 by hand, but this particular integral is appearing so often in so many
contexts that the integral is a special function, called the Error function and written erf(x). In a
code, you can call erf(x). The erf function is found in the math module.
3 http://en.wikipedia.org/wiki/Error_function.
141
v_sum = 0
for i in range(1, n, 1):
t = a + i*dt
v_sum = v_sum + v(t)
numerical = dt*(0.5*v(a) + v_sum + 0.5*v(b))
V = lambda t: exp(t**3)
exact_value = V(b) - V(a)
error = abs(exact_value - numerical)
rel_error = (error/exact_value)*100
print(’n={:d}: {:.16f}, error: {:g}’.format(n, numerical, error))
Unfortunately, the two other problems (1. and 3.) remain, and they are fundamental.
Computing Another Integral Suppose you next want to compute another integral,
say
1.1
−1 e −x 2 dx, using the previous specific implementation as your “starting
point”. What changes are required in the code then?
First of all, an anti-derivative can not (easily) be found 2,3 for this new integrand,
so we drop computing the integration error, and must remove the corresponding
code lines. In addition,
• the notation should be changed to fit the new problem. Thus, t and dt should be
replaced by x and h. Also, the integrand is (most likely) not a velocity any more,
so the name v should be changed with, e.g., f. Similarly, v_sum should rather be
f_sum then.
• the formula for v (or f) must be replaced by a new formula
• the limits a and b must be changed
These changes are straightforward to implement, but they are scattered around in
the program, a fact that requires us to be very careful so we do not introduce new
programming errors while we modify the code. It is also very easy to forget one or
two of the required changes.
For the sake of comparison, we might see how easy it is to rather use our general
implementation in trapezoidal.py for the task. With the following interactive
session, it should be clear that this implementation allows us to compute the new
integral
1.1
−1 e −x 2 dx without touching the implemented mathematical algorithm! We
can simply do:
In [1]: from trapezoidal import trapezoidal # ...general implementation
In [2]: from math import exp
In [3]: trapezoidal(lambda x: exp(-x**2), -1, 1.1, 400)
Out[3]: 1.5268823686123285
2 You cannot integrate e −x 2 by hand, but this particular integral is appearing so often in so many
contexts that the integral is a special function, called the Error function and written erf(x). In a
code, you can call erf(x). The erf function is found in the math module.
3 http://en.wikipedia.org/wiki/Error_function.