Romberg.jl

A simple Julia package to perform Romberg integration over discrete 1-dimensional data.

Romberg integration combines trapezoidal integration with Richardson extrapolation for improved accuracy. This package is meant to be used for integrating discrete data sampled with equal spacing. If your integrand can be evaluated at arbitrary points, then other methods are probably a better choice, e.g. QuadGK.jl. Similarly, if you have discrete data sampled at generic unequally spaced points, you probably need to use a low-order method like Trapz.jl.

First, install with the Julia package manager:

] add Romberg

The Romberg module exports a single function, romberg(x,y), or alternatively romberg(Δx,y), that returns a tuple (I,E) of the estimated integral I and a rough upper bound E on the error.

using Romberg

x = range(0, stop=π, length=2^8+1)
y = sin.(x)

romberg(x, y)

Unlike Trapz.jl, Romberg integration is a Newton-Cotes formula which requires each element of x be equally spaced. This is indirectly enforced in Romberg by requiring x::AbstractRange, or directly by passing the spacing Δx between points.

Romberg integration works by recursively breaking the integral down into trapezoidal-rule evaluations using larger and larger spacings Δx and then extrapolating back towards Δx → 0. This works by factorizing length(x) - 1, and therefore works best when length(x) - 1 has many small factors, ideally being a power of two.

(In the event that length(x) - 1 is prime, the romberg function is nearly equivalent to the trapezoidal rule, since it extrapolates only from 2 points to length(x) points.)

Currently Romberg only allows integration over a single dimension, so y::AbstractVector.

Given the limitations of Romberg, why use it over Trapz? For discrete samples of an underlying smooth function, Romberg integration can obtain significantly more accurate estimates at relatively low additional computational cost over trapezoidal integration for a given number of samples. Moreover, unlike the trapezoidal rule, the romberg function also returns an error estimate.

Here are some examples:

using BenchmarkTools
using Trapz, Romberg

x = range(0, stop=π, length=2^6+1)
y = sin.(x)
exact_answer = 2

tans = trapz(x, y)
rans, _ = romberg(x, y)

julia> exact_answer - tans
0.0004016113599627502

julia> exact_answer - rans
4.440892098500626e-16

julia> @btime trapz($x, $y);
  340.834 ns (1 allocation: 96 bytes)

julia> @btime romberg($x, $y);
  515.078 ns (1 allocation: 192 bytes)

So romberg is ~30% slower than trapz, but achieves nearly machine-precision accuracy, ~12 digits more accurate than trapz. Even if 500 times as many samples of the function were to be used in trapz, it would still be ~7 digits less accurate than romberg.

x = range(0, stop=1, length=2^4+1)
y = x.^3
exact_answer = 0.25

tans = trapz(x, y)
rans, _ = romberg(x, y)

julia> exact_answer - tans
-0.0009765625

julia> exact_answer - rans
0.0

romberg is able to obtain the exact answer (and in general is exact for polynomials of sufficiently low degree), compared to ~3 digits of accuracy for trapz, at the cost of ~2× the run time.

m = 3
n = 4
x = range(0, stop=π, length=2^6+1)
y = sin.(m*x).*cos.(n*x)
exact_answer = 2*m/(m^2 - n^2)

tans = trapz(x, y)
rans, _ = romberg(x, y)

julia> exact_answer - tans
0.0012075513578178043

julia> exact_answer - rans
-1.2385595582475872e-7

romberg is ~4 digits better accuracy than trapz and ~50% greater run time.