Taking advantage of Julia's integration abilities

Question

One of the main reasons I wanted to use Julia for my project is because of its speed, especially for calculating integrals.

I would like to integrate a 1-d function f(x) over some interval [a,b]. In general Julia's quadgk function would be a fast and accurate solution. However, I do not have the function f(x), but only its values f(xi) for a discrete set of points xi in [a,b], stored in an array. The xi's are regularly spaced, and I can get the spacing to be however small I like.

Naively, I could simply define a function f which interpolates using the values f(xi) and feed this to quadgk, (and make the spacing as small as possible), however then I won't know what my error is, which is a shame because QuadGK tells you the error in its estimation.

Another solution is to write a function myself to integrate the array (with trapezoid rule for example), but that would defeat the purpose of using Julia...

What is the easiest way to accurately integrate a function only given discrete values using Julia?

Given discrete values and no other information about the function, you can't do better than trapezoid. — Oscar Smith
@OscarSmith : I am new to Julia. For now, I prefer to use already existing integration functions in Julia than write my own. Apparently there is a function QuadGK but again, same problem I guess this only accepts functions, not an array of values? — Joshua Benabou
When you say that you can get the spacing as small as you like, that indicates that you have access to the underlying function, is that correct? Or is it measured from some physical process? If you have any knowledge about the underlying signal (or is it function), you may be able to improve on trapz. Is it smooth? Is it oscillatory? Trapz has some systematic issues with oscillations if you are actually integrating f(x)^2. — DNF

Oscar Smith Oscar Smith · Accepted Answer · 2020-06-05T03:02:11

Since you only have values, not the function itself, trapezoid will be your best bet probably. The package Trapz provides this (https://github.com/francescoalemanno/Trapz.jl). However, I think it is worth seeing how easy writing a pretty good implementation yourself would be.

function trap(A)
    return sum(A) - (A[begin] + A[end])/2
end

This takes 2.9ms for an array of 10 million floats. If they're Int, then 2.9ms. If they were complex numbers, it would still work (and take 8.9 ms)

A method like this is a good example to show how simple it can be to write pretty fast code in Julia that is still fully generic

Taking advantage of Julia's integration abilities

1 Answers