search

dextorious / NumericalIntegration.jl

Basic numerical integration routines for presampled data.
Other
62 stars 21 forks source link

improved implementation of Romberg integration #30

Open stevengj opened 3 years ago

stevengj commented 3 years ago

I posted an implementation on discourse (https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/8?u=stevengj) using Richardson.jl that is much more general (is not restricted to power-of-two sizes), seems to be more accurate than your current code, and is vastly simpler. I'm not sure if you would be interested in adopting it here?

My implementation could be sped up by inlining the sum calls and a few other tweaks, if necessary

dextorious commented 3 years ago

I responded on Discourse, but I do like the algorithm quite a bit. I'll leave this open for now as a reminder to benchmark it later.

stevengj commented 3 years ago

Note that the priority of my implementation is more about accuracy and generality than performance, though the performance of the code now in Romberg.jl isn't too bad (but still a bit slower than yours, maybe a factor of 2 IIRC).

(My sense is that usually the computation of the integrand is the limiting factor in performance for quadrature routines, not the integration code itself.)

dextorious commented 3 years ago

Your sense is no doubt correct in general, but I wrote this package to explicitly address the case where you have presampled data and can't cheaply re-evaluate it at different points.

I haven't looked at the final Romberg.jl code yet, but from what I saw in the Discourse thread it I thought it could be made quite a bit faster still. It could be quite interesting if it could match or surpass the classic implementation. And even if that's not the case, I might just add a RombergAccurate or RombergGeneral method that calls Romberg.jl.

stevengj commented 3 years ago

I wrote this package to explicitly address the case where you have presampled data and can't cheaply re-evaluate it at different points.

Even in that case, obtaining the data (e.g. from parsing a CSV file, or measurements, or from a finite-difference simulation) is probably more expensive than the quadrature. I find it hard to imagine a situation where you are generating discretely sampled data so rapidly that the performance of the quadrature is the rate-limiting step.

I haven't looked at the final Romberg.jl code yet, but from what I saw in the Discourse thread it I thought it could be made quite a bit faster still.

The final code is quite a bit faster.