6.4 Vectorizing the Functions
147
The evaluation points in the midpoint method are x i = a +
h
2 + ih, i = 0, . . . , n− 1.
That is, n uniformly distributed coordinates between a + h/2 and b − h/2. Such
coordinates can be calculated by x = linspace(a+h/2, b-h/2, n). Given that
the Python implementation f of the mathematical function f works with an array
argument, which is very often the case in Python, f(x) will produce all the function
values in an array. The array elements are then summed up by sum, when calling
sum(f(x)). The resulting sum is to be multiplied by the rectangle width h to
produce the integral value. The complete function is listed below.
from numpy import linspace, sum
def midpoint(f, a, b, n):
h = (b-a)/n
x = linspace(a + h/2, b - h/2, n)
return h*sum(f(x))
The code is found in the file integration_methods_vec.py.
Let us test the code interactively in a Python shell by computing
1
0 3t 2 e t 3 dt.
The file with the code above has the name integration_methods_vec.py and is
a valid module from which we can import the vectorized function:
In [1]: from integration_methods_vec import midpoint
In [2]: from numpy import exp
In [3]: v = lambda t: 3*t**2*exp(t**3)
In [4]: midpoint(v, 0, 1, 10)
Out[4]: 1.7014827690091872
Note the necessity to use exp from numpy: our v function will be called with x as
an array, and the exp function must be capable of working with an array.
The vectorized code performs all loops very efficiently in compiled code,
resulting in much faster execution. Moreover, many readers of the code will also
say that the algorithm looks clearer than in the loop-based implementation.
6.4.2 Vectorizing the Trapezoidal Rule
We can use the same approach to vectorize the trapezoidal function. However,
the trapezoidal rule performs a sum where the end points have different weight. If
we do sum(f(x)), we get the end points f(a) and f(b) with a weight of unity
instead of one half. A remedy is to subtract the error from sum(f(x)): sum(f(x))
- 0.5*f(a) - 0.5*f(b). The vectorized version of the trapezoidal method then
becomes (the code is found in integration_methods_vec.py)
def trapezoidal(f, a, b, n):
h = (b-a)/n
x = linspace(a, b, n+1)
s = sum(f(x)) - 0.5*f(a) - 0.5*f(b)
return h*s
147
The evaluation points in the midpoint method are x i = a +
h
2 + ih, i = 0, . . . , n− 1.
That is, n uniformly distributed coordinates between a + h/2 and b − h/2. Such
coordinates can be calculated by x = linspace(a+h/2, b-h/2, n). Given that
the Python implementation f of the mathematical function f works with an array
argument, which is very often the case in Python, f(x) will produce all the function
values in an array. The array elements are then summed up by sum, when calling
sum(f(x)). The resulting sum is to be multiplied by the rectangle width h to
produce the integral value. The complete function is listed below.
from numpy import linspace, sum
def midpoint(f, a, b, n):
h = (b-a)/n
x = linspace(a + h/2, b - h/2, n)
return h*sum(f(x))
The code is found in the file integration_methods_vec.py.
Let us test the code interactively in a Python shell by computing
1
0 3t 2 e t 3 dt.
The file with the code above has the name integration_methods_vec.py and is
a valid module from which we can import the vectorized function:
In [1]: from integration_methods_vec import midpoint
In [2]: from numpy import exp
In [3]: v = lambda t: 3*t**2*exp(t**3)
In [4]: midpoint(v, 0, 1, 10)
Out[4]: 1.7014827690091872
Note the necessity to use exp from numpy: our v function will be called with x as
an array, and the exp function must be capable of working with an array.
The vectorized code performs all loops very efficiently in compiled code,
resulting in much faster execution. Moreover, many readers of the code will also
say that the algorithm looks clearer than in the loop-based implementation.
6.4.2 Vectorizing the Trapezoidal Rule
We can use the same approach to vectorize the trapezoidal function. However,
the trapezoidal rule performs a sum where the end points have different weight. If
we do sum(f(x)), we get the end points f(a) and f(b) with a weight of unity
instead of one half. A remedy is to subtract the error from sum(f(x)): sum(f(x))
- 0.5*f(a) - 0.5*f(b). The vectorized version of the trapezoidal method then
becomes (the code is found in integration_methods_vec.py)
def trapezoidal(f, a, b, n):
h = (b-a)/n
x = linspace(a, b, n+1)
s = sum(f(x)) - 0.5*f(a) - 0.5*f(b)
return h*s