Taking advantage of Julia's integration abilities

user_7077761 0
Viewed 395

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?

julia integration
1 Answers

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

user_51413280
Related