Romberg.jl: Simple Romberg integration in Julia

Romberg integration combines trapezoidal integration with Richardson extrapolation for significantly more accurate estimates at slightly more cost than trapezoidal integration alone.

This package uses Trapz.jl for the trapezoidal integration, and then calculates a basic extrapolation.

Like trapz, romberg is meant to be used for integrating discrete points. For integrating a function, there are much better packages available. Additionally, there are a few limitations to Romberg:

  1. The points must be equally spaced along x
  2. The number of points - 1 must be a power of 2
  3. 1-dimensional integration only

Practically, limitation 2 is the most significant, but it is probably possible to relax this in a future update. If these three conditions can be met, romberg is much more accurate than trapz.

I originally considered making a pull request to add this functionality to Trapz.jl, but decided that it makes more sense to keep it separate as it’s a different integration technique. I plan to register it soon.

10 Likes

very happy you used my package for building a higher order method! :smiley:

you can get around limitation 2 via:

using Romberg
using Trapz

function arb_romberg(x,y)
    L=length(x)
    L<=1 && return zero(eltype(x)) * zero(eltype(x))
    pL=2^(leading_zeros(oftype(L,0)) - leading_zeros(L)-1)+1
    if pL<=L
        return romberg(x[1:pL], y[1:pL]) + arb_romberg(x[pL:end],y[pL:end])
    else
        return trapz(x[1:end], y[1:end])
    end
end
vx=0:0.05:1
vy=@. vx*vx
arb_romberg(vx,vy) # 0.33333333333333326
trapz(vx,vy) # 0.33375
romberg(vx,vy) #currently a DomainError

(this code is not optimized)

Very interesting! Regarding this point:

If your integrand is in functional form, then other methods are probably a better choice, e.g. QuadGK.jl or HCubature.jl, among others.

is this a limitation of the idea of romber integration, or is it more of an implementation problem?

I wonder if https://github.com/JuliaMath/Richardson.jl could be used to make a robust functional version of this.

How does this compare to https://github.com/dextorious/NumericalIntegration.jl

1 Like

@francesco.alemanno Thanks for writing the Trapz.jl package!

Your suggestion for arb_romberg is very clever. I had only thought about padding the vectors with zeros, like an FFT, but that wouldn’t work here because it would cause a discontinuity in the integrand. I’ve coded up a modified version of arb_romberg and built it into the package. It works, although it is slower than if the length(x)-1 is a power of 2.

@Mason the method works by applying trapezoidal integration using just 2 samples, then it applies trapezoidal integration using 3 samples, then 5, 9, etc, all the way to the 2^n+1 samples. The extrapolation then occurs between results of those integrations and they are combined in such a way that there is cancellation of errors. The formula for this combination is quite specific, so there is no advantage to using the more general Richardson.jl for extrapolation, even if a functional integrand is available.

@jtravs I was unaware of NumericalIntegration.jl. I needed Romberg integration for something I’m working on, but nothing came up when searching on discourse or Julia Packages.

NumericalIntegration.jl has a Romberg integration routine and it is significantly faster than Romberg.jl. In my effort to take advantage of Trapz.jl, I think I handicapped myself when it comes to benchmarking. Most of the time in Romberg.jl is spent on creating views of the data and running trapz on the view. NumericalIntegration.jl has coded up the entire procedure in a single function which is much more efficient.

I hate to say it after coding up Romberg.jl, but NumericalIntegration.jl should probably be preferred.

7 Likes

julia 1.5 will significantly improve the performance of view code

1 Like

you may also try to @inline the smaller functions. As long as a view does not cross a function barrier its creation will be elided and e.g. correct indexing computed by the compiler