146
6 Computing Integrals and Testing Code
Remark
The trapezoidal and midpoint methods are just two examples in a jungle of
numerical integration rules. Other famous methods are Simpson’s rule and
Gauss quadrature. They all work in the same way:
b
a
f (x)dx ≈
n−1
i=0
w i f (x i ) .
That is, the integral is approximated by a sum of function evaluations, where
each evaluation f (x i ) is given a weight w i . The different methods differ in
the way they construct the evaluation points x i and the weights w i . Higher
accuracy can be obtained by optimizing the location of x i .
6.4 Vectorizing the Functions
The functions midpoint and trapezoidal usually run fast in Python and compute
an integral to satisfactory precision within a fraction of a second. However, long
loops in Python may run slowly in more complicated implementations. To increase
speed, the loops can be replaced by vectorized code. The integration functions offer
simple and good examples on how to vectorize loops.
We have already seen simple examples on vectorization in Sect. 1.5, when we
evaluated a mathematical function f (x) for a large number of x values stored in an
array. Basically, we can write
def f(x):
return exp(-x)*sin(x) + 5*x
from numpy import exp, sin, linspace
x = linspace(0, 4, 101) # coordinates from 100 intervals on [0, 4]
y = f(x)
# all points evaluated at once
The result y is an array that, alternatively, could have been computed by running a
for loop over the individual x values and called f for each value. Vectorization essentially eliminates this explicit loop in Python (i.e., the looping over x and application of f to each x value are instead performed in a library with fast, compiled code).
6.4.1 Vectorizing the Midpoint Rule
We start by vectorizing the midpoint function, since trapezoidal is not equally
straightforward to vectorize. In both cases, our vectorization will remove the explicit
loop. The fundamental ideas of the vectorized algorithm are to
1. compute and store all the evaluation points in one array x
2. call f(x) to produce an array of corresponding function values
3. use the sum function to sum up the f(x) values
6 Computing Integrals and Testing Code
Remark
The trapezoidal and midpoint methods are just two examples in a jungle of
numerical integration rules. Other famous methods are Simpson’s rule and
Gauss quadrature. They all work in the same way:
b
a
f (x)dx ≈
n−1
i=0
w i f (x i ) .
That is, the integral is approximated by a sum of function evaluations, where
each evaluation f (x i ) is given a weight w i . The different methods differ in
the way they construct the evaluation points x i and the weights w i . Higher
accuracy can be obtained by optimizing the location of x i .
6.4 Vectorizing the Functions
The functions midpoint and trapezoidal usually run fast in Python and compute
an integral to satisfactory precision within a fraction of a second. However, long
loops in Python may run slowly in more complicated implementations. To increase
speed, the loops can be replaced by vectorized code. The integration functions offer
simple and good examples on how to vectorize loops.
We have already seen simple examples on vectorization in Sect. 1.5, when we
evaluated a mathematical function f (x) for a large number of x values stored in an
array. Basically, we can write
def f(x):
return exp(-x)*sin(x) + 5*x
from numpy import exp, sin, linspace
x = linspace(0, 4, 101) # coordinates from 100 intervals on [0, 4]
y = f(x)
# all points evaluated at once
The result y is an array that, alternatively, could have been computed by running a
for loop over the individual x values and called f for each value. Vectorization essentially eliminates this explicit loop in Python (i.e., the looping over x and application of f to each x value are instead performed in a library with fast, compiled code).
6.4.1 Vectorizing the Midpoint Rule
We start by vectorizing the midpoint function, since trapezoidal is not equally
straightforward to vectorize. In both cases, our vectorization will remove the explicit
loop. The fundamental ideas of the vectorized algorithm are to
1. compute and store all the evaluation points in one array x
2. call f(x) to produce an array of corresponding function values
3. use the sum function to sum up the f(x) values