# Romberg.jl: Simple Romberg integration in Julia

**URL:** <https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313>\
**Category:** Package Announcements\
**Created:** [May 11, 2020, 11:48pm UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313 "2020-05-11T23:48:55Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![EP-Guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ep-guy/32/3442_2.png) [@EP-Guy](https://discourse.julialang.org/u/EP-Guy)\
**Post date:** [May 11, 2020, 11:48pm UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/1 "2020-05-11T23:48:55Z")

</div>

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

> **[GitHub - fgasdia/Romberg.jl: Romberg integration routine in Julia](https://github.com/fgasdia/Romberg.jl)**
>
> Romberg integration routine in Julia. Contribute to fgasdia/Romberg.jl development by creating an account on GitHub.

This package uses [Trapz.jl](https://github.com/francescoalemanno/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.

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 12, 2020, 4:04pm UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/2 "2020-05-12T16:04:45Z")

</div>

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

you can get around limitation 2 via:

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

---

<div class="post-metadata">

**Author:** ![Mason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mason/32/2423_2.png) [@Mason](https://discourse.julialang.org/u/Mason)\
**Post date:** [May 12, 2020, 6:26pm UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/3 "2020-05-12T18:26:14Z")

</div>

Very interesting! Regarding this point:

> If your integrand is in functional form, then other methods are probably a better choice, e.g. [QuadGK.jl](https://github.com/JuliaMath/QuadGK.jl) or [HCubature.jl](https://github.com/JuliaMath/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 [GitHub - JuliaMath/Richardson.jl: Richardson extrapolation in Julia](https://github.com/JuliaMath/Richardson.jl) could be used to make a robust functional version of this.

---

<div class="post-metadata">

**Author:** ![jtravs](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jtravs/32/2010_2.png) [@jtravs](https://discourse.julialang.org/u/jtravs)\
**Post date:** [May 12, 2020, 10:35pm UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/4 "2020-05-12T22:35:30Z")

</div>

How does this compare to [GitHub - dextorious/NumericalIntegration.jl: Basic numerical integration routines for presampled data.](https://github.com/dextorious/NumericalIntegration.jl)

---

<div class="post-metadata">

**Author:** ![EP-Guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ep-guy/32/3442_2.png) [@EP-Guy](https://discourse.julialang.org/u/EP-Guy)\
**Post date:** [May 13, 2020, 3:56am UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/5 "2020-05-13T03:56:13Z")

</div>

@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](https://juliapackages.com/).

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.

---

<div class="post-metadata">

**Author:** ![francesco.alemanno](https://avatars.discourse-cdn.com/v4/letter/f/e8c25b/32.png) [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)\
**Post date:** [May 17, 2020, 5:04pm UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/6 "2020-05-17T17:04:38Z")

</div>

> [@EP-Guy](#):
>
> Most of the time in Romberg.jl is spent on creating views of the data and running `trapz` on the view.

julia 1.5 will significantly improve the performance of view code

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [May 17, 2020, 7:06pm UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/7 "2020-05-17T19:06:08Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 8, 2020, 7:06pm UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/8 "2020-12-08T19:06:37Z")

</div>

Sorry to revive an old thread, but I have some updates and a question about the best place for future development.

> [@EP-Guy](#):
>
> @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 2n+12^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.

Actually, contrary to the typical presentation of Romberg integration, the method is _not_ limited to power-of-two+1 sizes, and instead can work with any factorization. Moreover, if you use the latest version of the Richardson.jl package, it can handle this for you and results in a small and elegant implementation:

```julia
import Primes, Richardson

@views function romberg(Δx::Real, y::AbstractVector; kws...)
    n = length(y)
    n > 1 || throw(ArgumentError("y array must have length > 1"))
    endsum = (y[begin]+y[end])/2
    n == 2 && return endsum * Δx
    factors = Primes.factor(n-1)
    vals = [(Δx * (sum(y[begin+1:end-1]) + endsum), Δx)]
    sizehint!(vals, sum(values(factors)))
    step = 1
    for (f,K) in factors
        for k = 1:K
            step *= f
            sΔx = step*Δx
            pushfirst!(vals, (sΔx * (sum(y[begin+step:step:end-step]) + endsum), sΔx))
        end
    end
    return Richardson.extrapolate(vals, power=2; kws...)
end

romberg(x::AbstractRange, y::AbstractVector; kws...) = romberg(step(x), y; kws...)

```

Not only is the above code quite a bit shorter than the implementations in Romberg.jl and NumericalIntegration.jl, it also has the following advantages:

- It supports arbitrary array lengths n, although it works best when n-1 is highly composite.

- It returns the estimated integral as well as an error estimate.

- It supports integrands in arbitrary normed vector spaces.

- Even for n = 2^k + 1 lengths, it seems to be more accurate than the implementation in NumericalIntegration.jl (see below for an example).

- In general, Richardson.jl implements more refined extrapolation than most implementations of this algorithm, because it compares the estimated errors from every level of the Neville tableau and picks the estimate with the smallest error.

```julia
julia> x = range(0, 10, length=65); y = cos.(x);

julia> val, err = romberg(x, y)
(-0.5440211107599316, 2.0984771309517924e-7)

julia> val - sin(10) # exact error
1.2943812688348544e-10

julia> val2 = NumericalIntegration.integrate(x, y, NumericalIntegration.RombergEven())
┌ Warning: RombergEven :: final step reached, but accuracy not: 0.5440207498733947 > 1.0e-12
└ @ NumericalIntegration ~/.julia/packages/NumericalIntegration/3mVSe/src/NumericalIntegration.jl:175
-0.5440205148227985

julia> val2 - sin(10) # worse error than above
5.960665713233837e-7

```

As mentioned above, it still works well when n-1 is any highly composite size:

```julia
julia> x = range(0, 10, length=61); y = cos.(x);

julia> val, err = romberg(x, y)
(-0.544021113675439, 4.007245680170968e-6)

julia> val - sin(10)
-2.786069264182345e-9

```

The performance typically seems a bit worse than `NumericalIntegration`, probably because it has more dynamic allocations — but could be sped up further, e.g. by replacing the `sum` calls over views with direct loops and omitting the call to `Primes.factor` in the power-of-two case.

**Question:** should I submit a PR? If so, to which package?

---

<div class="post-metadata">

**Author:** ![EP-Guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ep-guy/32/3442_2.png) [@EP-Guy](https://discourse.julialang.org/u/EP-Guy)\
**Post date:** [December 8, 2020, 8:21pm UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/9 "2020-12-08T20:21:17Z")

</div>

I’d be willing to accept the new algorithm into Romberg.jl — but I also don’t mind if others think it makes more sense to put it in NumericalIntegration.jl.

I opted not to register Romberg.jl after realizing it is outperformed by NumericalIntegration’s more classic and faster implementation, but I would register it if we decide the algorithm should go there.

---

<div class="post-metadata">

**Author:** ![EP-Guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ep-guy/32/3442_2.png) [@EP-Guy](https://discourse.julialang.org/u/EP-Guy)\
**Post date:** [December 14, 2020, 5:54pm UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/10 "2020-12-14T17:54:09Z")

</div>

With further thought, I’d encourage you to submit the PR to Romberg.jl.

The existence of Trapz.jl already gives some precedence to having Romberg integration in its own package as opposed to another method in NumericalIntegration.jl. Additionally, the owner of NumericalIntegration.jl hasn’t had public Github activity for over 5 months, so it’s unknown when a PR there would be accepted.

I’ll tag a version of Romberg.jl as it is, but with this PR, much of the current package will change.

---

<div class="post-metadata">

**Author:** ![EP-Guy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ep-guy/32/3442_2.png) [@EP-Guy](https://discourse.julialang.org/u/EP-Guy)\
**Post date:** [December 20, 2020, 6:31am UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/11 "2020-12-20T06:31:18Z")

</div>

I’m happy to announce that [Romberg.jl](https://github.com/fgasdia/Romberg.jl) is now a registered package and can be installed with

```julia
] add Romberg

```

from the REPL. `v0.1.0` is the old implementation and `v0.2.0` is the new one. Huge thanks to @stevengj for providing the new faster and more accurate implementation.

---

<div class="post-metadata">

**Author:** ![dextorious](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dextorious/32/20490_2.png) [@dextorious](https://discourse.julialang.org/u/dextorious)\
**Post date:** [December 24, 2020, 1:46am UTC](https://discourse.julialang.org/t/romberg-jl-simple-romberg-integration-in-julia/39313/12 "2020-12-24T01:46:32Z")

</div>

As the owner in question, it’s true that I haven’t been very active on GitHub lately, but I do usually respond within two weeks or so. That said, NumericalIntegration.jl has been feature complete for my own purposes for a long time, so I mostly just accept PRs when people propose new things. I also don’t mind accepting other people as maintainers or transferring the package to an appropriate organization if the community would like more active maintenance.

As for this particular algorithm, it looks quite interesting and I’ll benchmark it when I get a chance. If it’s considerably better than what NumericalIntegration.jl has now, I might just add Romberg.jl as a dependency.