Information Technology Reference
In-Depth Information
programming language. The linear algebra field in mathematics defines some vec-
tor operations, such as inner (scalar) products, matrix-vector products and so on.
Programming languages support these mathematical operations, but also offer addi-
tional functionality. Unfortunately, there is no standard naming convention for this
additional functionality. The vectorization of an algorithm is therefore more tightly
coupled to its actual implementation in a particular programming language.
The Trapezoidal Rule in Vectorized Form
Vectorization of Algorithm
6.2
requires rewriting the instructions such that the loop
is replaced by built-in operations working on vectors. Let us start with a vector
x
,
going from
a
to
b
in steps of
h D .a b/=n
:
x
1
D a;
x
2
D a C h;:::
x
j
D a C .j 1/h;:::
x
nC1
D b:
We then apply
f
to each element in
x
such that
f.x/
is a new vector
v
of the
same length as
should be performed in Fortran, C, or
C
CC
for efficiency. Languages with support for vector operations normally supply
a function for summing all the elements of the vector
f.x/
:
s D s
u
m
nC1
iD1
x
. This application of
f
to
x
f
i
.This
operation is often named
reduce
or
sum
. The quantity
s
is not exactly what we
want, since the end points
f
1
1
2
and not 1, as implied
in the sum over all elements
f
i
. The correction is easily performed by subtracting
1
2
and
f
nC1
should have weight
.f
1
Cf
nC1
/
. Finally, we need to multiply
s
by
h
. Algorithm
6.13
summarizes the
various steps (sum
.
v
/
here is our notation for adding all elements in
v
).
Algorithm 6.13
Vectorized Trapezoidal Integration.
Trapezoidal
vec (
a
,
b
,
f
,
n
)
h D
ba
n
x D .a; a C h;:::;b/
v
D f.x/
s D h .
sum
.
v
/ 0:5.
v
1
C
v
nC1
//
return
s
Vectorized Numerical Integration in Python
Translating Algorithm
6.13
into Numerical Python is straightforward. Generation
of the sequence of evaluation points
.a; aCh;:::;b/
is performed by the statement
x = arrayrange(a, b, h, Float)
Because of round-off errors the upper limit can end up as
b-h
instead of
b
.Itis
therefore smart to set the upper limit as
b+h/2
; this will guarantee that
b
becomes
