Information Technology Reference
In-Depth Information
the last element in
x
. The sum operation on a vector
v
is called
add.reduce(v)
in Numerical Python. Knowing about
arrayrange
and
add.reduce
, it should be
quite easy to understand the implementation of the vectorized algorithm in Python:
def trapezoidal(a, b, f, n):
h = (b-a)/float(n)
x = arrayrange(a, b+h/2, h, Float)
v = f(x)
r=h
*
(sum(v) - 0.5
*
(v[0] + v[-1]))
return r
Indexing Python arrays is a slow process, but we do it only twice, which is not
critical if
n
is large. Arrays in Python start with index
0
,and
-1
is a notation for the
last index.
Numerical Python redefines the multiplication operator
*
such that
h
*
f(x)
is a
scalar
h
multiplied by a vector
f(x)
. The expression
f(x)
means applying a func-
tion
f
to a vector
x
, element by element. However, we must avoid implementing
f(x)
as an explicit loop in Python. To this end, Numerical Python supports the effi-
cient application of standard mathematical functions and operators on vectors. For
example,
sin(x)
computes the sine of each element in a vector
x
. We can therefore
implement a function
f.x/ D e
x
2
ln
.1 C x
sin
x/
applied to each element of a
vector
x
as
def f1(x):
f = exp(-x
*
x)
*
log(1+x
*
sin(x))
return f
This is exactly the same syntax as that used for applying a scalar function on a scalar
value. Sending, e.g., a number 3.1 to
f1
, results in the computation
e
3:1
2
ln
.1 C
3:1
sin
3:1/
. Sending a vector
x
to
f1
results in more complicated computations
“behind the curtain”:
1.
temp1 = sin(x)
, i.e., apply the sine function to each entry in
x
2.
temp2 = 1 + temp1
, i.e., add 1 to each element in
temp1
3.
temp3 = log(temp2)
, i.e., compute the natural logarithm of all elements in
temp2
4.
temp4 = x
*
x
, i.e., square each element in
x
5.
temp5 = exp(temp4)
, i.e., apply the exponential function to each element in
temp4
6.
f = temp5
*
temp3
, i.e., multiply
temp5
and
temp3
element by element, as in a
scalar multiplication
As we see,
f
is built of several vector operations, requiring temporary arrays, and
hence additional storage, compared with a straightforward loop. However, a straight
loop computing
f
element by element required 40 times as long CPU time as the
vectorized version (on an Intel laptop computer with Linux and Python compiled by
GNU's C compiler
gcc
).
The main program calling up the vectorized version of the
trapezoidal
function
is exactly the same program as we used to test the scalar version of
Trapezoidal
:
