Page - 76 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

76 3 Computing Integrals 3.5 Vectorization The functionsmidpoint andtrapezoidusually run fast in Python and compute an integral to a satisfactory precisionwithin a fractionof a second. However, long loops inPythonmayrunslowly inmorecomplicated implementations. To increase the speed, the loops canbe replacedbyvectorizedcode. The integration functions constitute a simpleandgoodexample to illustratehowtovectorize loops. We have already seen simple examples on vectorization in Sect. 1.4whenwe couldevaluateamathematical functionf.x/ for a largenumberofx values stored inanarray.Basically,wecanwrite 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 the array that would be computed if we ran a for loop over the individualxvaluesandcalledf foreachvalue.Vectorizationessentially eliminates this loop inPython(i.e., the loopingoverxandapplicationoff to eachxvalueare insteadperformedina librarywith fast, compiledcode). Vectorizing the midpoint rule The aim of vectorizing the midpoint and trapezoidal functions is also to remove the explicit loop in Python. We start with vectorizing themidpoint function since trapezoid is not equally straight- forward. The fundamental ideasof thevectorizedalgorithmare to 1. computeall theevaluationpoints inonearrayx 2. callf(x) toproduceanarrayofcorrespondingfunctionvalues 3. use thesum function tosumthef(x)values Theevaluationpoints in themidpointmethodarexi DaC.iC12/h, i D0;:: :;n 1. That is,nuniformlydistributedcoordinatesbetweenaCh=2andb h=2. Such coordinatescanbecalculatedbyx = linspace(a+h/2, b-h/2, n). Given that thePython implementationfof themathematical functionf workswith an array argument,whichisveryoftenthecase inPython,f(x)willproduceall thefunction values in an array. The array elements are then summed up bysum: sum(f(x)). This sumis tobemultipliedby the rectanglewidthh toproduce the integralvalue. Thecomplete function is listedbelow. from numpy import linspace, sum def midpoint(f, a, b, n): h = float(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.