# Integration of functions over geometries

**URL:** <https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994>\
**Category:** Specific Domains\
**Tags:** question, fem, mesh\
**Created:** [May 31, 2022, 2:04pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994 "2022-05-31T14:04:04Z")\
**Posts on this page:** 11\
**Page:** 2

<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:** [January 30, 2024, 4:42pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/21 "2024-01-30T16:42:01Z")

</div>

> [@mike.ingold](#):
>
> It looks like the implicit caveat to the “16-digit accuracy and in `O(n)` time” bumper sticker is that this is only precisely true for functions that can be exactly represented by a polynomial of order `2n-1` over the specified interval.

No, that’s not the case. For polynomials Gaussian quadrature is _exact_, but for analytic functions Gaussian quadrature [converges exponentially fast](https://people.maths.ox.ac.uk/trefethen/publication/PDF/2008_127.pdf) with the number of quadrature points n, so you can indeed get to nearly machine precision (nearly 16 digits) rapidly, and if you need to go to large n then FastGaussQuadrature.jl computes the quadrature rule in O(n) time.

A downside of FastGaussQuadrature.jl is that it just gives you an n-point rule, it doesn’t tell you what n should be for a given tolerance and a given function, and the rule is not a “nested” rule — it doesn’t give you an error estimate. So you are only your own in finding the correct n for your problem, e.g. by repeatedly doubling n until it converges. (Repeatedly increasing the order is also called p-adaptive quadrature.)

QuadGK is h-adaptive, in that it repeatedly subdivides the integration interval until a tolerance is achieved, and it uses an [embedded rule](https://en.wikipedia.org/wiki/Gauss%E2%80%93Kronrod_quadrature_formula) to get an error estimate. This is convenient, and it is also fairly robust if your integrand is badly behaved in some localized region because it will only refine the quadrature around that region.

If your function is analytic, then using high-order quadrature is theoretically more efficient, but on the other hand you rarely have to go to very high order in practice, and the time for the computation of the rule is less important if you are doing integration over and over for many different functions and/or domains. So for smooth functions you can simply pass `order=n` for a larger `n` to QuadGK.jl to use a higher-order rule (which is computed once and cached). The computation of the rule is O(n^2), so it will be much slower than FastGaussQuadrature.jl for large n, but since in practice you rarely need n larger than a few hundred it doesn’t matter too much in my experience.

The other big advantage of FastGaussQuadrature.jl is that it has built-in rules for weighted quadrature with many different weight functions, which allow you to efficiently handle _known_ endpoint singularities. See the discussion of [weighted Gauss rules in the QuadGK manual](https://juliamath.github.io/QuadGK.jl/stable/weighted-gauss/), and [integrands with singularities and discontinuities](https://juliamath.github.io/QuadGK.jl/stable/quadgk-examples/#Integrands-with-singularities-and-discontinuities).

PS. Although QuadGK is not asymptotically exponentially convergent for a fixed order, unlike repeatedly increasing the order, it often gives the _illusion_ of exponential convergence up to a given precision. See also this post: [(Insanely-) Higher Order Derivatives - #17 by stevengj](https://discourse.julialang.org/t/insanely-higher-order-derivatives/104455/17)

---

<div class="post-metadata">

**Author:** ![mike.ingold](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mike.ingold/32/203749_2.png) [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)\
**Post date:** [January 30, 2024, 10:26pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/22 "2024-01-30T22:26:28Z")

</div>

It’s for exactly these reasons that I generally tend toward h-adaptive quadrature/cubature methods. These days I find myself writing a lot of tooling that will be used to run calculations/models, moreso than direct data-in-hand problems. That being the case, I find it hugely beneficial to have an interface where I can specify some acceptable `atol` or `rtol` and get a result within those specs reliably since I can’t know at time of writing what function/signal will even be integrated.

---

<div class="post-metadata">

**Author:** ![mike.ingold](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mike.ingold/32/203749_2.png) [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)\
**Post date:** [January 30, 2024, 10:53pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/23 "2024-01-30T22:53:08Z")

</div>

Updates on [LineIntegrals](https://github.com/mikeingold/LineIntegrals.jl):

- I’ve re-implemented all of the existing `integral` methods using FGQ.
- All `QuadGK.quadgk(f, <:Meshes.Geometry)` methods remain defined. Again, this is type piracy so would probably need to be PR’d into another package.
- I’m planning to implement an internal function with in-place memory re-use for geometries with multiple segments.
- The README has a toy example integrating `f(x,y) = abs(x+y)` over a Bezier curve approximation of the unit circle using both techniques. The `integrate` methods default to `n=100` for Gauss-Legendre nodes, and `quadgk` uses all default `kwargs`. Both are solved in similar time with similar allocations, but QuadGK.jl seems to be more accurate. Increasing `n` with `integral` causes its result to converge toward `quadgk`’s, but at the expense of a significant increase in allocations. Not a scientific test by any means; your mileage may vary.

---

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [January 30, 2024, 11:45pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/24 "2024-01-30T23:45:34Z")

</div>

Hi,

I’m new to the thread, so I may have missed some details, and I wanted to share a package that allows integration over polytopes using nested integration as an algorithm, which is convenient for integration in low-dimensions (up to 3 or 4), and even ideal for nearly singular integrands. In two (and more) dimensions the nested integration algorithm evaluates the multivariate integral in terms of recursive 1d integrals, i.e:

\int\_{a\_x}^{b\_x} dx\ \int\_{a\_y(x)}^{b\_y(x)} dy\ f(x,y) = \int\_{a\_x}^{b\_x} dx\ I\_y(x), \quad I\_y(x) = \int\_{a\_y(x)}^{b\_y(x)} dy\ f(x,y)

What is nice about this algorithm is that it is compatible with any quadrature scheme and any geometry (polytopes, manifolds) as long as either a parametrization of the geometry is available, or optimization on that geometry is possible and fast compared to evaluation of the integrand.  
The example below is with the package [AutoBZCore.jl](https://github.com/lxvm/AutoBZCore.jl), which emulates the Integrals.jl interface, and it allows for integration over polytopes from Polyhedra.jl

```julia
julia> using AutoBZCore, Polyhedra

# make convex hull in 3d
julia> domain = load_limits(polyhedron(vrep(rand(10, 3))))
PolyhedraExt.PolyhedralLimits{3, Float64, DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}}(convexhull([0.900543, 0.738892, 0.00729183], [0.165838, 0.909174, 0.316936], [0.552774, 0.774557, 0.908199], [0.906346, 0.161918, 0.456898], [0.979034, 0.303522, 0.0330761], [0.146472, 0.13728, 0.183881], [0.673617, 0.694496, 0.738391], [0.380986, 0.0822024, 0.346046], [0.894882, 0.748271, 0.0540279], [0.40436, 0.03022, 0.151078]))

# setup integrand
julia> prob = AutoBZCore.IntegralProblem((x,p) -> sin(prod(x)), domain)
AutoBZCore.IntegralProblem{var"#5#6", PolyhedraExt.PolyhedralLimits{3, Float64, DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}}, AutoBZCore.NullParameters}(var"#5#6"(), PolyhedraExt.PolyhedralLimits{3, Float64, DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}}(convexhull([0.900543, 0.738892, 0.00729183], [0.165838, 0.909174, 0.316936], [0.552774, 0.774557, 0.908199], [0.906346, 0.161918, 0.456898], [0.979034, 0.303522, 0.0330761], [0.146472, 0.13728, 0.183881], [0.673617, 0.694496, 0.738391], [0.380986, 0.0822024, 0.346046], [0.894882, 0.748271, 0.0540279], [0.40436, 0.03022, 0.151078])), AutoBZCore.NullParameters())

# call quadgk recursively
julia> AutoBZCore.solve(prob, NestedQuad(QuadGKJL()))
AutoBZCore.IntegralSolution{Float64, Float64}(0.019284903620528718, 1.7385351242394514e-18, true, -1)

julia> using FastGaussQuadrature

# use a product rule of 10-point Gauss-Legendre quadratures
julia> AutoBZCore.solve(prob, NestedQuad(QuadratureFunction(fun=gausslegendre, npt=10)))
AutoBZCore.IntegralSolution{Float64, Nothing}(0.01928490362052873, nothing, true, -1)

# pick quadgk for two variables and Gauss-Legendre for the last
julia> AutoBZCore.solve(prob, NestedQuad(QuadGKJL(), QuadratureFunction(fun=gausslegendre, npt=10), QuadGKJL()))
AutoBZCore.IntegralSolution{Float64, Float64}(0.01928490362052872, 5.454892180317694e-19, true, -1)

```

My hope is to integrate this functionality with Integrals.jl (no pun intended) when I find the time, and the main idea of the implementation is in  
[IteratedIntegration.jl](https://github.com/lxvm/IteratedIntegration.jl/blob/main/ext/PolyhedraExt.jl).

---

<div class="post-metadata">

**Author:** ![mike.ingold](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mike.ingold/32/203749_2.png) [@mike.ingold](https://discourse.julialang.org/u/mike.ingold)\
**Post date:** [January 31, 2024, 12:05am UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/25 "2024-01-31T00:05:56Z")

</div>

Nice work! I’d actually been wondering earlier, among the discussions about the different integration packages, if the better interface wouldn’t be implement this functionality into Integrals.jl. The core functionality that I’m really looking for, personally, is the ability to integrate functions over some standardized 1D/2D/3D geometry definitions, ideally with an easy way to derive the geometry from a CAD mesh model.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [January 31, 2024, 11:21am UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/26 "2024-01-31T11:21:32Z")

</div>

> [@mike.ingold](#):
>
> The core functionality that I’m really looking for, personally, is the ability to integrate functions over some standardized 1D/2D/3D geometry definitions, ideally with an easy way to derive the geometry from a CAD mesh model.

Same here @mike.ingold. The ability to integrate functions over CAD-derived geometries is really missing, and you are closing this gap 🙏

As soon as we finish the new coordinate system infrastructure in Meshes.jl, we will be able to address very practical problems such as computing areas and volumes of complicated domains, which are not made of polytopes.

---

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [February 3, 2024, 4:59am UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/27 "2024-02-03T04:59:28Z")

</div>

@juliohm Out of curiosity, does Meshes.jl implement elimination, i.e. finding a lower-dimensional hull from the intersection of a polytope and a hyperplane, and double description, i.e. computing half-spaces from vertices and vice-versa? Both of these are all you need for applying iterated integration to convex polytopes.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [February 3, 2024, 12:34pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/28 "2024-02-03T12:34:04Z")

</div>

@lxvm we don’t have those. It would be nice to review PRs with these additions 🙂

---

<div class="post-metadata">

**Author:** ![aplavin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aplavin/32/222056_2.png) [@aplavin](https://discourse.julialang.org/u/aplavin)\
**Post date:** [February 3, 2024, 12:55pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/29 "2024-02-03T12:55:00Z")

</div>

Does “integration of functions over geometries” include 2d integration over arbitrarily defined regions, through automatic variable change?  
Like, computing \int\_a^b dx \int\_{lo(x)}^{hi(x)} f(x, y) dy.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [February 3, 2024, 1:02pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/30 "2024-02-03T13:02:12Z")

</div>

@aplavin as discussed above, we are considering all sorts of integrations over geometries. People are working on different fronts to make these more readily available.

---

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [February 3, 2024, 4:07pm UTC](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994/31 "2024-02-03T16:07:20Z")

</div>

Yes, that is the goal. In this direction, [IteratedIntegration.jl](https://lxvm.github.io/IteratedIntegration.jl/stable/methods/) defines an interface for that form of limits of integration and uses Polyhedra.jl to implement it for integration over polytopes. The package is still experimental in its algorithms, but that interface on integration domains will remain stable and could be extended to more geometries that are convex and support optimization for finding the extrema of the coordinate axes on that manifold.

[Previous page](https://discourse.julialang.org/t/integration-of-functions-over-geometries/81994.md?page=1)