Page - 147 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

6.4 VectorizingtheFunctions 147 Theevaluationpoints in themidpointmethodarexi =a+h2+ih, i=0,.. .,n−1. That is, n uniformly distributed coordinates between a+h/2 and b−h/2. Such coordinatescan be calculated byx = linspace(a+h/2, b-h/2, n). Given that the Python implementationf of the mathematical functionf works with an array argument,whichisveryoftenthecase inPython,f(x)willproduceall thefunction 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 producethe integralvalue.Thecomplete 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)) Thecode is found in thefileintegration_methods_vec.py. Let us test the code interactively in a Python shell by computing ∫1 0 3t 2et 3 dt. Thefile with thecodeabovehas thenameintegration_methods_vec.pyand is avalidmodule fromwhichwecan import thevectorizedfunction: 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 useexp fromnumpy: ourv function will be called withx as anarray,and theexp functionmustbe capableofworkingwith anarray. 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 thealgorithmlooksclearer than in the loop-basedimplementation. 6.4.2 VectorizingtheTrapezoidalRule 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 insteadofonehalf.Aremedyis to subtract theerror fromsum(f(x)):sum(f(x)) - 0.5*f(a) - 0.5*f(b).The vectorizedversionof the trapezoidalmethod then becomes(thecode is found inintegration_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