Simpsons Rule

[simonluzara]

Hello everyone, am a new enthusiastic member of julia community. I have spent quite a fair share of time watching tutorials from Chris and all the others.

I have a matlab simpsons rule code which was very generic for numerical integration for both vector of a given length and a function as below;

function I = simpsons(f,a,b,n)

% This function computes the integral "I" via Simpson's rule in the interval [a,b] with n+1 equally spaced points
% 
% Syntax: I = simpsons(f,a,b,n)
% 
% Where,
%  f= can be either an anonymous function (e.g. f=@(x) sin(x)) or a vector
%  containing equally spaced values of the function to be integrated
%  a= Initial point of interval
%  b= Last point of interval
%  n= # of sub-intervals (panels), must be integer
% 
%  Written by Juan Camilo Medina  - The University of Notre Dame
%  09/2010 (copyright Dr. Simpson)
% 
% 
% Example 1:
% 
% Suppose you want to integrate a function f(x) in the interval [-1,1].
% You also want 3 integration points (2 panels) evenly distributed through the
% domain (you can select more point for better accuracy).
% Thus:
%
% f=@(x) ((x-1).*x./2).*((x-1).*x./2);
% I=simpsons(f,-1,1,2)
% 
% 
% Example 2:
% 
% Suppose you want to integrate a function f(x) in the interval [-1,1].
% You know some values of the function f(x) between the given interval,
% those are fi= {1,0.518,0.230,0.078,0.014,0,0.006,0.014,0.014,0.006,0}
% Thus:
%
% fi= [1 0.518 0.230 0.078 0.014 0 0.006 0.014 0.014 0.006 0];
% I=simpsons(fi,-1,1,[])
%
% note that there is no need to provide the number of intervals (panels) "n",
% since they are implicitly specified by the number of elements in the
% vector fi

if numel(f)>1 % If the input provided is a vector
    n=numel(f)-1; h=(b-a)/n;
    I= h/3*(f(1)+2*sum(f(3:2:end-2))+4*sum(f(2:2:end))+f(end));
else % If the input provided is an anonymous function
    h=(b-a)/(2*n); xi=a:h:b;
    I= h/3*(f(xi(1))+2*sum(f(xi(3:2:end-2)))+4*sum(f(xi(2:2:end)))+f(xi(end)));
end

I have tried to write in in Julia as a training task but have got quite few errors now and then.

function simps(f,a,b,n)
    if isvector(f)
        n=length(f)-1
        h=(b-a)/n
        I= h/3*(f[1]+2*sum(f[3:2:end-2])+4*sum(f[2:2:end])+f[end])
    else
        h=(b-a)/(2*n)
        xi=a:h:b
        I= h/3*(f(xi[1])+2*sum(f(xi[3:2:end-2]))+4*sum(f(xi[2:2:end]))+f(xi[end]));
    end
    return I
end

I tried it with g(x)=x^3, then simply writing simps(g(),0,1) but doesn’t seem to work.

I will be very grateful for your insights.

[Jeff_Emanuel]

@simonluzara

sum(f(xi[3:2:end-2]))

There’s a method for sum that takes a function and an array. Then your function doesn’t need to know how to operate on the array. Or you could broadcast f over xi.

[lmiq]

I didn’t check for correctness, but this is one way that could be written in more Julian style:

@views function simps(x::AbstractVector,a,b)
    n=length(x)-1
    h=(b-a)/n
    I = h/3*(x[1]+2*sum(x[3:2:end-2])+4*sum(x[2:2:end])+x[end])
    return I
end 
    
function simps(f::Function,a,b,n)
    h=(b-a)/(2*n)
    fx=f.(a:h:b)
    return simps(fx,a,b)
end 

The first method receives the vector, and the second method receives the function and the additional parameter, and computes the range, which is passed then to the first method.

ps: added @views to the first function because otherwise the slices like x[3:2:end-2] allocate intermediate arrays.

ps2: edited, adding the computation of fx = f.(a:h:b), which computes the array of the results of applying the function, thanks @Jeff_Emanuel. But then this can be optimized out (soon here).

[Jeff_Emanuel]

@lmiq

function simps(f::Function,a,b,n)
    h=(b-a)/(2*n)
    x=a:h:b
    return simps(x,a,b)
end 

It appears you dropped f

[lmiq]

Ok, now this is better:

@views function simps(f::Function, x::AbstractVector)
    n = length(x) - 1
    h = (x[end]-x[begin])/n
    I= h/3*(f(x[1])+2*sum(f,x[3:2:end-2])+4*sum(f,x[2:2:end])+f(x[end]))
    return I
end
simps(f::Function, a::Real, b::Real, n::Integer) = simps(f, a:((b-a)/n):b)
simps(x::AbstractVector) = simps(identity, x)

The first method receives the function to be integrated and the domain of points, as a vector. It can be called with:

julia> g(x) = x^3
g (generic function with 1 method)

julia> simps(g, 0:1/10:1)
0.25

The second method receives the function, the interval and the number of bins. It just defines the number of steps and calls the first one. It can be used with:

julia> simps(g, 0, 1, 10)
0.25

The third method just calls the first passing the identity function as parameter:

julia> x = collect((0:1/10:1).^3);

julia> simps(x)
0.25

None is allocating intermediates here:

julia> @btime simps($g, 0:1/10:1)
  317.585 ns (0 allocations: 0 bytes)
0.25

julia> @btime simps($g, 0, 1, 10)
  318.183 ns (0 allocations: 0 bytes)
0.25

julia> @btime simps($x)
  30.869 ns (0 allocations: 0 bytes)
0.25

Note that I’ve used sum(f, range) as pointed by @Jeff_Emanuel .

Note, @simonluzara, that we could have used a structure more similar to that of the Matlab code, with a an “if vector else …” construct. However, that is not the most interesting way to write this code in Julia, thanks to multiple dispatch, which is what I have done there. With multiple dispatch we decide what to do from the type of the variable being provided to the code, and that eliminates the conditional.

You can, then, easily, add additional methods to the same function, for instance if you like to be able to provide a range explicitly to the integration, we can just define:

simps(f::Function, r::UnitRange, n::Integer) = simps(f, r[begin], r[end], n)

which will call the first method when used with:

julia> simps(g, 0:1, 10)
0.25

You could also define now simps(f, x::AbstractMatrix) to integrate on a greater dimension (although probably the method itself would be generalized in this case).

[stevengj]

@lmiq

h = (x[end]-x[begin])/n

This seems wrong to me — it assumes all of the entries in the x array are equally spaced. It seems like you need to restrict yourself to AbstractRange subtypes, in which case you can use step(x).

@lmiq

a:((b-a)/n):b

I would recommend using range(a, b, length=n+1) to ensure that you have exactly n+1 points (n intervals), rather than being potentially susceptible to rounding errors.

@lmiq

simps(x::AbstractVector) = simps(identity, x)

This also seems like not what is intended by the original poster (who wanted to pass an array of function values for equally-spaced points).

[lmiq]

Third try

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

# from the function and range definition
function simps(f::Function, a::Real, b::Real, n::Integer) 
    return simps(f, range(a, b, length=n))
end

# from the function values and the step length
@views function simps(fx::AbstractVector, h::Real)
    I= h/3*(fx[1]+2*sum(fx[3:2:end-2])+4*sum(fx[2:2:end])+fx[end])
    return I
end

# from the function values and range
simps(fx::AbstractVector, x::AbstractRange) = simps(fx, step(x))

Well, I only intended to exemplify multiple dispatch, ended learning new things.

julia> g(x) = x^5
g (generic function with 1 method)

julia> simps(g, 0:1/1000:1)
0.166666666667

julia> simps(g, 0, 1, 1000)
0.16700450199899838

julia> fx = (0:1/1000:1).^5;

julia> simps(fx, 1/1000)
0.166666666667

julia> simps(fx, 0:1/1000:1)
0.166666666667

[simonluzara]

That’s rather not straight forward, a bit of jargon there.

[Jeff_Emanuel]

Methods and dynamic dispatch are probably the most distinguishing feature of Julia:
https://docs.julialang.org/en/v1/manual/methods/

Docs for the sum method that accepts a function:
https://docs.julialang.org/en/v1/base/collections/#Base.sum

Broadcasting is convenient syntax to iteratively apply some operation:
https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting

HTH

[simonluzara]

Thank you very much, it looks like I have got still a long cliff to climb getting my expert level as I was in matlab.
What does @views does?

[lmiq]

A slice of an array like x[1:3] allocates, in Julia, a new array. With views you get only a reference to the slice of the original array, without new allocations.

Your original code was almost functional, and it could be a more or less direct translation of the Matlab code. It is just not the most natural and perhaps elegant and efficient use of Julia.

[simonluzara]

I am quite grateful for this, let’s work on Simpsons 8th rule which is more accurate to that of 3rd.

Second thing, I think I might have been wrong all the way from the start. The Simpson 1/3 rule indicates that N must be even number, making the function values uneven. This is to say we are restricted to some cases when let’s say we have the velocities of 100 points, and each value of the point. I thought putting 2*N would force it to make the spacing even, but the 100 points vectors stay the same hence f[1]+4*sum(f[2:2:end-1])+2*sum(f[3:2:end-2])+f[end] misses a point if f[:] is even.

[lmiq]

Just a final comment: Your original code could “work” just with two small changes, to fix how to test if the input is a vector, and how to broadcast the application of a function to all the elements of an array:

julia> function simps(f,a,b,n)
           if f isa AbstractVector # >>>>  test here if f is a vector
               n=length(f)-1
               h=(b-a)/n
               I= h/3*(f[1]+2*sum(f[3:2:end-2])+4*sum(f[2:2:end])+f[end])
           else
               h=(b-a)/(2*n)
               xi=a:h:b
               # >>> in the next line, use the notation f.(x[range]) (note the dot)
               I= h/3*(f(xi[1])+2*sum(f.(xi[3:2:end-2]))+4*sum(f.(xi[2:2:end]))+f(xi[end])) 
           end
           return I
       end
simps (generic function with 1 method)

julia> g(x) = x^5;

julia> simps(g, 0, 1, 1000)
0.1666666666666875

julia> fx = (0:1/1000:1).^5;

julia> simps(fx, 0, 1, 1000)
0.166666666667

Yet, particularly because of the lack of @views, this is slower than the above alternatives:

julia> @btime simps($g, 0, 1, 1000)
  152.772 μs (2 allocations: 15.88 KiB)
0.1666666666666875

vs, with the above alternative:

julia> @btime simps($g, 0, 1, 1000)
  75.962 μs (0 allocations: 0 bytes)
0.16700450199899838

If you add @views to the corrected version of your function, the performance is essentially the same, for a single call. If you were calling this from a loop millions of times, the other options are better, because at least you don’t need to check if the input is a vector every time (although the compiler may well deal with that in some cases). The greatest advantage of the dispatch based versions is that you can add methods to your function without changing the already working code.

[simonluzara]

Thank you very much. This is indeed helpful. isvector written in a different way.

[simonluzara]

I got this solution form the math exchange. I think if we are able to add this, then the most correct Simpson’s Integration non bias code would be created in Julia thanks to us.

[lmiq]

Note that there are very robust numerical integration packages in Julia already, for instance:

GitHub - SciML/Integrals.jl: A common interface for quadrature and numerical...

https://juliamath.github.io/QuadGK.jl/latest/

[simonluzara]

I have seen this one already, quads are used mainly on integrating FEM. And the guy didn’t make the quadratures standard.
If you are interested, we could work on this Simpson’s (most of it is actually done).

[lmiq]

I’m not particularly interested. But go for it! It will be a nice learning exercise, and you can write a small package with the results.

[stevengj]

@lmiq

at least you don’t need to check if the input is a vector every time (although the compiler may well deal with that in some cases).

No, type checks like f isa AbstractVector are optimized out by the compiler, since types are known at compile-time. (Remember, a specialized version of each function is compiled for each argument type that is passed.)

Defining separate methods for different argument types is more idiomatic and more extensible, however.

[simonluzara]

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.