# 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.

## Usage

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)
```

## Limitations

`x`

Equally spaced 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.

`length(x) - 1`

a power of 2

Most effective if 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.)

### 1-dimensional

Currently `Romberg`

only allows integration over a single dimension, so
`y::AbstractVector`

.

`trapz`

Comparison to 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:

### 1)

```
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`

.

### 2)

```
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.

### 3)

```
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.