Simpsons Rule

New to Julia
, ,
21

I have found an optimum solution with the easiest and reasonable syntax. Thanks to @lmiq and @Jeff_Emanuel for great insights.

# integral from the function and range
@views function simps(f::Function, x::AbstractRange) #A range with odd number
    h = step(x)
    I= h/3*(f(x[1])+2*sum(f,x[3:2:end-2])+4*sum(f,x[2:2:end-1])+f(x[end]))
    return I
end

# from the function and range definition
function simps(f::Function, a::Real, b::Real, n::Integer) #n as an even number
    return simps(f, range(a, b, length=n+1))
end

@views function simps(fx::AbstractVector, h::Real)
    if length(fx)%2==1
        I= h/3*(fx[1]+2*sum(fx[3:2:end-2])+4*sum(fx[2:2:end-1])+fx[end])
    else
        I=h/3*(fx[1]+2*sum(fx[3:2:end-5])+4*sum(fx[2:2:end-4])+fx[end-3])+
        (3*h/8)*(fx[end-3]+3fx[end-2]+3fx[end-1]+fx[end])
    end
    return Float64(I)
end

# from the function values and range
function simps(fx::AbstractVector, a::Real, b::Real)
    return simps(fx, (b-a)/(length(fx)-1))
end

simps (generic function with 5 methods)

With respect to the Simpson’s rule that N intervals should be “even”, thus the vector values length odd.
I have mixed both 1/3 rule and 3/8 rule for the case where N intervals are odd.

The memory allocations both tested (in a i5 4740H with 8Gb RAM):

y=LinRange(0,1,1000).^3
g(x)=x^3

using BenchmarkTools;
@btime simps($y,0,1)
@btime simps($g,0,1,1000)

  913.158 ns (0 allocations: 0 bytes)
  2.867 μs (0 allocations: 0 bytes)

I hope this will be useful to future research works on modelling FVM or FDM PDE’s…as in my case it gives an accuracy of 1e-8.

2 Likes