Correct the TrapezoidalEven() method and optimize the performance a little bit. (#8)

* Make integrations work for complex numbers as well.

* Correct the formula of the TrapezoidalEven method.

* Optimize the performance by reducing some unnecessary operations.

Author: i2000s

src/NumericalIntegration.jl

@@ -28,7 +28,7 @@ end
 
 function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::TrapezoidalEven)
     @assert length(x) == length(y) "x and y vectors must be of the same length!"
-    return 0.5 * (x[end] - x[1]) / (length(y) - 1) * (y[1] + y[end] + sum(y[2:end-1]))
+    return (x[2] - x[1]) * (0.5 * (y[1] + y[end]) + sum(y[2:end-1]))
 end
 
 function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::TrapezoidalFast)
@@ -45,7 +45,7 @@ function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector
     @fastmath @simd for i in 2 : N
         @inbounds retval += y[i]
     end
-    @inbounds return (x[end] - x[1]) / N * (retval + 0.5*y[1] + 0.5*y[end])
+    @inbounds return (x[2] - x[1]) * (retval + 0.5*y[1] + 0.5*y[end])
 end
 
 function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::SimpsonEven)
@@ -54,7 +54,7 @@ function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector
     for i in 5 : length(y) - 1
         retval += y[i]
     end
-    return (x[end] - x[1]) / (length(y) - 1) * retval
+    return (x[2] - x[1]) * retval
 end
 
 function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}, ::SimpsonEvenFast)
@@ -64,7 +64,7 @@ function integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector
     @fastmath @inbounds for i in 5 : length(y)-1
         retval += y[i]
     end
-    @inbounds return (x[end] - x[1]) / (length(y) - 1) * retval
+    @inbounds return (x[2] - x[1]) * retval
 end
 
 integrate{X<:Number, Y<:Number}(x::AbstractVector{X}, y::AbstractVector{Y}) = integrate(x, y, TrapezoidalFast())

test/runtests.jl

@@ -4,7 +4,10 @@ using Base.Test
 methods = [Trapezoidal(), TrapezoidalEven(), TrapezoidalFast(), TrapezoidalEvenFast(), SimpsonEven(), SimpsonEvenFast()]
 x = collect(-π : π/1000 : π)
 y = sin.(x)
+p = collect(0 : π/1000 : π)
+q = sin.(p)
 for method in methods
     println(string("Testing method: ", typeof(method)))
     @test abs(integrate(x, y, method)) < 1e-4
+    @test abs(integrate(p, q, method)-2.0) < 1e-4
 end