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
that returns a tuple
(I,E) of the estimated integral
I and a rough upper bound
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
Romberg by requiring
x::AbstractRange, or directly by passing the
Δx between points.
Most effective if
length(x) - 1 a power of 2
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
Romberg only allows integration over a single dimension, so
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)
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
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
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.