Just being curious, I compared these two version on accuracy and execution
time, having grown up in Matlab as well. Let's call the first, vectorized
version by Tomas as 'trapz1', the second version by Cameron as 'trapz2',
see below. Then I get the following results:
julia> x = linspace(0, pi, 1000000);
julia> y = sin(x);
julia> trapz1(x, y)
1.999999999998355
julia> trapz2(x, y)
1.9999999999984162
julia> @time for i in 1:100; trapz1(x, y); end
elapsed time: 4.858258575 seconds (5600105600 bytes allocated)
julia> @time for i in 1:100; trapz2(x, y); end
elapsed time: 0.41866061 seconds (1600 bytes allocated)
and in general the unvectorized version is more than 10 times faster. But
why are the two results different (though only by 6.1e-14)? It looks like
they are performing exactly the same computations. Does the sequence of
elementary operations play a role?
Is it to be expected that a version with a localized length is faster? I
didn't see significantly different running times.
PS: The function definitions I used here are:
function trapz1{T<:Number}(x::Vector{T}, y::Vector{T})
local n = length(x)
if length(y) != n
error("Vectors 'x', 'y' must be of same length.")
end
if n == 1; return 0.0; end
sum((x[2:end]-x[1:end-1]) .* (y[2:end]+y[1:end-1])) / 2.0
end
function trapz2{T<:Number}(x::Vector{T}, y::Vector{T})
local n = length(x)
if (length(y) != n)
error("Vectors 'x', 'y' must be of same length")
end
if n == 1; return 0.0; end
r = 0.0
for i in 2:n
r += (x[i] - x[i-1]) * (y[i] + y[i-1])
end
r / 2.0
end
On Wednesday, April 23, 2014 11:20:50 PM UTC+2, Tomas Lycken wrote:
>
>
> On Wednesday, April 23, 2014 11:10:15 PM UTC+2, Cameron McBride wrote:
>>
>> Or you can use the non-vectorized version and save the overhead of the
>> temporary arrays being created by the addition and multiplication steps.
>>
>
> There's really no way I can hide that I learnt scientific computing in
> Matlab, is there? :P
>
>>
>>
