Information Technology Reference
In-Depth Information
a = 0; b = 2; n = 1000
for i in range(10000):
result = trapezoidal(a, b, f1, n)
import time
t1 = time.clock()
print result, t1
We note that in Python you can say
print a
,where
a
is any built-in object (includ-
ing an array) and get nicely formatted output. This is very convenient during
program development and debugging.
The reader should notice how convenient a language using dynamic typing is;
many of the code segments developed for scalar computations can be directly reused
for vectorized computations. (The same functionality is obtained in C
CC
by the use
of templates and overloaded operators, though at the cost of increased language and
technical complexity.)
The vectorized version of the trapezoidal method uses one-seventh of the CPU
time required by the loop-based counterpart. This means that vectorized Python runs
slightly faster than Java and at about half the speed of Fortran 77, C, or C
CC
.
Vectorized Numerical Integration in Matlab
Algorithm
6.13
can be implemented directly in Matlab:
function r = Trapezoidal_vec(a, b, f, n)
% TRAPEZOIDAL Numerical integration from a to b
% by the Trapezoidal rule
f = fcnchk(f);
h = (b-a)/n;
x = [a:h:b];
v = feval(f, x);
r=h
*
(sum(v) - 0.5
*
(v(1) + v(length(v))));
The construction of the sequence of evaluation points having a uniform distance
h
between
a
and
b
is created by the expression
[a:h:b]
.The
feval(f,x)
expression
we know from the scalar (non-vectorized) version of the function can also be used
when
x
is a vector. As in Python, the return value from
f
is then a vector. Summing
up all elements in a vector
v
is performed with
sum(v)
, and indexing the last element
in
v
can be written as
v(length(v))
(if the first index in
v
is set to 1, which is the
default value).
The function to be integrated must now work for a vector. In the previous
Python example we could reuse a function implemented for a scalar input and
return value. Such reuse is not possible in Matlab when the function expression
f.x/ D e
x
2
log
.1 C x
sin
x/
. The reason is because the expression
exp(-x
*
x)
for a vector
x
implies a matrix multiplication between
x
and itself. What we need
is the
.
*
operator such that an element in the answer is the product of the corre-
sponding two elements in the operands. Hence, the vectorized version of the scalar
f1
function is
function y = f1_vec(x)
y = exp(-x .
*
x) .
*
log(1 + x .
*
sin(x));
