Simpson's rule in Julia

An approximation to the integral of a function $f(x)$ over an interval $[a,b]$ can be approximated by the Simpson's rule as follows:

Using the composite Simpson's rule, the interval $[a,b]$ is divided in subintervals and Simpson's rule is applied to each interval and the results are summed. The composite Simpson's rule is given by

where $n$ is the number (even) of subintervals and $h=(b-<!--\; -\; -->a)/n$.

The composite Simpson's rule of a function f over the interval [a,b] taking n intervals can be implemented in Julia (0.4.2) as

function simps(f::Function, a::Number, b::Number, n::Number)
    n % 2 == 0 || error("`n` must be even")
    h = (b-a)/n
    s = f(a) + f(b)
    s += 4sum(f(a + collect(1:2:n) * h))
    s += 2sum(f(a + collect(2:2:n-1) * h))
    return h/3 * s
end

As an example, ${\int <!--\; \int \; -->}_0^{\mathrm{\pi <!--\; \pi \; -->}}3\sin <!--\; \; -->(x{)}^{3}dx=4$, in Julia:

julia> simps(x -> 3sin(x).^3,0,pi,4)
3.792237795874079

julia> simps(x -> 3sin(x).^3,0,pi,6)
3.9783434011671583

julia> simps(x -> 3sin(x).^3,0,pi,12)
3.998978466021562

Scipy provides a Simpson's rule for Python. It takes an array to integrate and an array with the sampling points or the spacing of the integrating points (the parameter $h$). Let's implement the function following a similar structure.

The number of points of the array corresponds with the number of intervals and must be odd, being the number of subdivisions even. Although Scipy deals with an even number of intervals, here we are going to keep it simple and just implement the mathematical definition. Some code can be simplified, for example using auxiliar functions, but it is kept to be more explicit.

Writing the composite Simpson's rule as

can be implemented in Julia for samples y and sampling spacing h as

function simps(y::Vector, h::Number)
    n = length(y)-1
    n % 2 == 0 || error("`y` length (number of intervals) must be odd")
    s = sum(slice(y,1:2:n) + 4slice(y,2:2:n) + slice(y,3:2:n+1))
    return h/3 * s
end

And using the same example as before for a vector of 201 elements:

julia> (y=3sin(linspace(0,pi,201)).^3; h=pi/200; simps(y,h))
3.9999999878202863

Similarly, knowing the samples array y and the sampling points array x:

function simps(y::Vector, x::Union{Vector,Range})
    n = length(y)-1
    n % 2 == 0 || error("`y` length (number of intervals) must be odd")
    length(x)-1 == n || error("`x` and `y` length must be equal")
    h = (x[end]-x[1])/n
    s = sum(slice(y,1:2:n) + 4slice(y,2:2:n) + slice(y,3:2:n+1))
    return h/3 * s
end

julia> (x=linspace(0,pi,201); y=3sin(x).^3; simps(y,x))
3.9999999878202863

Remark that the vectors are sliced using the function slice and not the direct indexing (array[start:step:stop]), avoiding the creation of new vectors. Also, x can be a Vector or a Range, like a LinSpace.

With the Scipy implementation in mind, let's add support for matrices by adding the parameter axis to integrate the matrix y by row (axis=2, default) or column (axis=1).

The implementation of the Simpson's rule integration for the samples y, knowing the sampling spacing h:

function simps(y::Matrix, h::Number, axis::Integer=2)
    n = size(y)[axis]-1
    n % 2 == 0 || error("`y` length (number of intervals) must be odd")
    axis in [1,2] || error("`axis` must be 1 or 2")
    if axis == 2
        inds = [(:,1:2:n),(:,2:2:n),(:,3:2:n+1)]
    else
        inds = [(1:2:n,:),(2:2:n,:),(3:2:n+1,:)]
    end
    s = sum(slice(y,inds[1]...) + 4slice(y,inds[2]...) + slice(y,inds[3]...),axis)
    return h/3 * s
end

The implementation of the Simpson's rule integration for the samples y, knowing the sampling points array x:

function simps(y::Matrix, x::Union{Matrix,Vector,Range}, axis::Integer=2)
    n = size(y)[axis]-1
    n % 2 == 0 || error("`y` length (number of intervals) must be odd")
    axis in [1,2] || error("`axis` must be 1 or 2")
    if length(size(x)) == 1
        length(x)-1 == n || error("`x` and `y` length must be equal in axis `axis`")
        h = (x[end]-x[1])/n
    else
        size(x)[axis]-1 == n || error("`x` and `y` length must be equal in axis `axis`")
        h = nothing
    end
    if axis == 2
        h != nothing || (h = (x[1,end]-x[1,1])/n)
        inds = [(:,1:2:n),(:,2:2:n),(:,3:2:n+1)]
    else
        h != nothing || (h = (x[end,1]-x[1,1])/n)
        inds = [(1:2:n,:),(2:2:n,:),(3:2:n+1,:)]
    end
    s = sum(slice(y,inds[1]...) + 4slice(y,inds[2]...) + slice(y,inds[3]...),axis)
    return h/3 * s
end

As we can see, x might be a Matrix but also a Vector or a Range, like in the previous case. The requirement is that it needs to have the same length as y in the axis axis. In both cases, the axis defines the order of the indices to slice and the axis of the summation.

Adding to the previous example that ${\int <!--\; \int \; -->}_0^{\mathrm{\pi <!--\; \pi \; -->}}\sin <!--\; \; -->(x)dx=2$ and using it for testing the last two functions gives:

julia> (y=[3sin(linspace(0,pi,201)).^3 sin(linspace(0,pi,201))]'; h=pi/200; simps(y,h))
2x1 Array{Float64,2}:
 4.0
 2.0

julia> (y=[3sin(linspace(0,pi,201)).^3 sin(linspace(0,pi,201))]; h=pi/200; simps(y,h,1))
1x2 Array{Float64,2}:
 4.0  2.0

julia> (x=linspace(0,pi,201); y=[3sin(x).^3 sin(x)]'; simps(y,x,2))
2x1 Array{Float64,2}:
 4.0
 2.0

julia> (x=linspace(0,pi,201); y=[3sin(x).^3 sin(x)]; simps(y,x,1))
1x2 Array{Float64,2}:
 4.0  2.0