# Integration over regions with varying boundaries, triangular domains and more

**URL:** <https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157>\
**Category:** Numerics\
**Tags:** question, package, symbolic, integral\
**Created:** [March 15, 2021, 1:46am UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157 "2021-03-15T01:46:30Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![Aleksandr\_Mikheev](https://avatars.discourse-cdn.com/v4/letter/a/b4bc9f/32.png) [@Aleksandr\_Mikheev](https://discourse.julialang.org/u/Aleksandr_Mikheev)\
**Post date:** [March 15, 2021, 1:46am UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/1 "2021-03-15T01:46:30Z")

</div>

This topic may be considered, in some sense, as a follow-up of my [previous thread](https://discourse.julialang.org/t/sympy-symbolics-and-hcubature/56513). Thanks to [stevengj](https://discourse.julialang.org/u/stevengj) and [ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas) I believe I more or less solved the “symbolic part” of my problem and could start working on the actual thing. In short, I need to numerically compute integrals of the kind  
f\_1(t,t',x) = \int\_{t\_0}^t \mathrm{d}\tau \int\_{t\_0}^{t'} \mathrm{d}\tau' \int\_0^{\infty} \mathrm{d}y \,F\_1(\tau,\tau',y,t,t',x)  
and  
f\_2(t,t',x) = \int\_{t}^{t\_2} \mathrm{d}\tau \int\_{t}^{t'} \mathrm{d}\tau' \int\_0^{\infty} \mathrm{d}y \,F\_2(\tau,\tau',y,t,t',x),  
with analytical expressions for F\_1 and F\_2 being given, on some grid \lbrace t\_i,t\_j',x\_k \rbrace. For concreteness, let’s say  
F\_1(\tau,\tau',y,t,t',x) = \cos[x(t - \tau)] \cos[(x + y)(\tau - \tau')]\cos[y(\tau' - t')] (x + y)^{-2} (\tau + t)^{-1} = F\_2(\ldots)   
As explained in the previous thread, actual expressions are way too lengthy (and also “more singular”, by the way) and are given by some product of other functions and derivatives of thereof, so I use `Symbolics.jl` to first get the expression and then convert it to a Julia-type function in order to perform the integration. For the latter, I decided to use the `HCubature.jl` package. So here’s my current attempt of computing f\_1 (with some particular choice of t\_0 and “infinity”):

```julia
using Symbolics
using HCubature

t_grid = collect(0.1:0.1:2);
t_pr_grid = collect(0.1:0.1:2);
x_grid = collect(0.1:0.1:2);

f1 = Array{Float64}(undef,length(t_grid),length(t_pr_grid),length(x_grid));

@variables tau, tau_pr, t, t_pr, x, y; 

F_1 = cos(x*(tau - t)) * cos((x + y)*(tau - tau_pr)) * cos(y*(tau_pr - t_pr)) * (x + y)^(-2) * (tau + t)^(-1);

function f1_eval(F_1::Num,t_cur::Real,t_pr_cur::Real,x_cur::Real)
    xmin = [1e-6,1e-6,1e-8];                                                      
    xmax = [t_cur,t_pr_cur,100];

    # F_1(tau,tau',y,t,t',x) function at the given grid point (t,t',x):
    F1 = substitute(F_1,Dict((t,t_pr,x) .=> (t_cur,t_pr_cur,x_cur)));   
    F1_func = build_function(F1,[tau,tau_pr,y],expression=Val{false});
    return hcubature(F1_func,xmin,xmax,rtol=sqrt(1e-14));
end

for i in 1:length(t_grid)
    for j in 1:length(t_pr_grid)
        for k in 1:length(x_grid)
            f1[i,j,k] = f1_eval(F_1,t_grid[i],t_pr_grid[j],x_grid[k])[1];
        end
    end
end

```

Here are my questions:

1. This sort of does the job, but I am pretty sure there is a room for improvement. In particular, instead of converting F\_1 at each grid point I could maybe convert F\_1 to a vector-type function which for a given argument `(tau,tau_pr,y)` returns an array F\_1(\tau,\tau',y,t\_i,t\_j',x\_k). I am sure that this can be done with `Symbolics.jl`, but

- I got a bit confused with two types of functions that `build_function` returns in this case: [“…one for out-of-place AbstractArray output and a second which is a mutating function”](https://symbolics.juliasymbolics.org/dev/manual/build_function/#Symbolics._build_function-Tuple%7BSymbolics.JuliaTarget,AbstractArray,Vararg%7BAny,N%7D%20where%20N%7D). I didn’t really understand how to work with either of them, to be honest;
- I don’t have fixed boundaries but rather, again, an array of integration limits. Is there any effective way of integrating an array of functions over an array of integration domains? Note that there is a one-to-one correspondence between the integrands and the integration limits.

1. So far I don’t have any good idea on what to do with the second integral, where the “temporal integration domain” is no longer a rectangle but a triangle, albeit a particularly simple one. Unfortunately, `HCubature.jl` does not support integration domains of such kind (unless I miss something), but maybe there is some trick to still do this using `HCubature.jl`? The obvious one would be to extend the integration domain to a full rectangle and compensate this by modifying F\_2 with a bunch of Heaviside functions. But I reckon that’s a very bad idea. Perhaps there is a better one? If not, are there any good Julia packages that support integration over triangular domains?

Thank you in advance!

---

<div class="post-metadata">

**Author:** ![Perrin\_Meyer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/perrin_meyer/32/15686_2.png) [@Perrin\_Meyer](https://discourse.julialang.org/u/Perrin_Meyer)\
**Post date:** [March 15, 2021, 4:07am UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/2 "2021-03-15T04:07:52Z")

</div>

Here are a lot of references on triangle schemes, in the (excellent) Python code Quadpy

[https://github.com/nschloe/quadpy#triangle-t2](https://github.com/nschloe/quadpy#triangle-t2)

this appears to be a Julia implementation of one of the Quadpy routines

[https://github.com/eschnett/GrundmannMoeller.jl](https://github.com/eschnett/GrundmannMoeller.jl)

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 15, 2021, 6:25am UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/3 "2021-03-15T06:25:34Z")

</div>

These integrals look a lot like convolution, are you sure you can’t reformulate tthem as such and then use FFTs?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 15, 2021, 9:19am UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/4 "2021-03-15T09:19:51Z")

</div>

> [@Aleksandr\_Mikheev](#):
>
> the “temporal integration domain” is no longer a rectangle but a triangle, albeit a particularly simple one. Unfortunately, `HCubature.jl` does not support integration domains of such kind (unless I miss something), but maybe there is some trick to still do this using `HCubature.jl` ?

Generally, a change of variables may be used to transform to a rectangular domain, but I did not look at the details of your particular problem to see if it makes sense.

---

<div class="post-metadata">

**Author:** ![Aleksandr\_Mikheev](https://avatars.discourse-cdn.com/v4/letter/a/b4bc9f/32.png) [@Aleksandr\_Mikheev](https://discourse.julialang.org/u/Aleksandr_Mikheev)\
**Post date:** [March 15, 2021, 1:16pm UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/5 "2021-03-15T13:16:22Z")

</div>

Unfortunately, a change of variables here would be a big pain in the neck.

---

<div class="post-metadata">

**Author:** ![Aleksandr\_Mikheev](https://avatars.discourse-cdn.com/v4/letter/a/b4bc9f/32.png) [@Aleksandr\_Mikheev](https://discourse.julialang.org/u/Aleksandr_Mikheev)\
**Post date:** [March 15, 2021, 8:55pm UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/6 "2021-03-15T20:55:38Z")

</div>

Hi! First of all, yes, the time integral is indeed a (double) convolution, albeit in a bit more general sense. Omitting the “spatial part”, it can be brought to F(t,t') = \int\_{-\infty}^{\infty} \mathrm{d}\tau \int\_{-\infty}^{\infty}\mathrm{d} \tau' A(t,\tau) B(\tau,\tau') C(\tau',t') by means of absorbing a couple of Heaviside functions into A and C.

Now, if A, B, and C only depend on difference of their arguments, then F(t,t') is also only a function of t - t' and one can easily show that F(\omega) = A(\omega) B(\omega) C(\omega). In that case, I agree that it would be sensible to compute Fourier images using FFT and then do the inverse Fourier transform to get F(t,t').

In my case, however, all the functions depend, in general, on both arguments (as I tried to emphasize in my “toy example” by adding the (\tau + t)^{-1} part). One can then invoke Wigner-Weyl calculus and introduces the so-called Wigner(-Weyl) transform, F\_W(T,\omega) = \int \mathrm{d} s\, \mathrm{e}^{-\mathrm{i}\omega s} F(T + s/2, T - s/2), to get F(T,\omega) = A(T,\omega) \star B(T,\omega) \star C(T,\omega), where \star denotes the [Moyal product](https://en.wikipedia.org/wiki/Moyal_product). Therefore, in principle, the time integration can be performed by first computing Wigner transforms of all the functions in the integrand, then taking their Moyal products (which also has an integral form, etc.) and finally performing the inverse Wigner transform. In theory, that might indeed be computationally less costly thing to do, but I am not sure if that will be more effective in practice given my poor programming skills. From what I see, I would then have to treat “spatial” and “temporal” integrals on a different footing, try to adapt existing libraries to compute things like Wigner transforms and Moyal products, etc. So unless I miss something and there is a much simpler solution, I would refrain myself from that approach if possible. Hope what I have written kind of makes sense.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [March 16, 2021, 9:55am UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/7 "2021-03-16T09:55:51Z")

</div>

Oh that’s interesting, I never saw the Wigner-Weyl calculus in time rather than in space. I guess it only is advantageous if you have a reason to believe you’re in a kind of semiclassical regime. Maybe you can try to differentiate wrt t and t’ and see if you can setup a differential equation. Otherwise, if F1 is completely non-separable you have no choice but to do the full nested integrals I guess: not sure you gain a lot by going to cubature.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 16, 2021, 10:36am UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/8 "2021-03-16T10:36:42Z")

</div>

This is a bit off-topic, but I cannot leave it be:

> [@Aleksandr\_Mikheev](#):
>
> ```julia-auto
> t_grid = collect(0.1:0.1:2);
> t_pr_grid = collect(0.1:0.1:2);
> x_grid = collect(0.1:0.1:2);
> 
> ```

Don’t `collect` ranges, just leave them as they are.

---

<div class="post-metadata">

**Author:** ![Aleksandr\_Mikheev](https://avatars.discourse-cdn.com/v4/letter/a/b4bc9f/32.png) [@Aleksandr\_Mikheev](https://discourse.julialang.org/u/Aleksandr_Mikheev)\
**Post date:** [March 16, 2021, 11:15am UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/9 "2021-03-16T11:15:25Z")

</div>

> Oh that’s interesting, I never saw the Wigner-Weyl calculus in time rather than in space.

Wigner-Weyl calculus in time domain is actually extremely common in nonequilibrium physics where there is no time-translation symmetry.

> I guess it only is advantageous if you have a reason to believe you’re in a kind of semiclassical regime.

Indeed, one option would be to essentially consider a “semicalssical regime” (in my case, it would have a different physical meaning, but I understand what you mean) by essentially expanding the Moyal product to, e.g., the first order – the Poisson bracket. But I would prefer to refrain myself from such approximations.

> Otherwise, if F1 is completely non-separable you have no choice but to do the full nested integrals I guess: not sure you gain a lot by going to cubature.

My integrands may be rather singular close to certain boundaries, so I thought it would be beneficial to use some adaptive scheme, and `HCubature` was essentially the first Julia package I saw that realizes an adaptive integration scheme. If there are other nice Julia libraries that you feel might be suitable for my purposes, I would be glad to hear about them, to be honest. 🙂

But even if I end up sticking with `HCubtautre`, I am almost certain that the code that I have written is far from being optimized.

---

<div class="post-metadata">

**Author:** ![Aleksandr\_Mikheev](https://avatars.discourse-cdn.com/v4/letter/a/b4bc9f/32.png) [@Aleksandr\_Mikheev](https://discourse.julialang.org/u/Aleksandr_Mikheev)\
**Post date:** [March 23, 2021, 4:53pm UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/10 "2021-03-23T16:53:27Z")

</div>

Sorry for such a late response. I was actually aware of the Python library by Nico, though didn’t know that there was also a Julia package for integrating over n-dimensional simplices. In both cases, however, I don’t immediately see how I can use it to compute, e.g., \int\_0^1 \mathrm{d}\tau \int\_0^{\tau} \mathrm{d} \tau' \int\_a^b \mathrm{d} x f(\tau,\tau',x). In theory, I could, e.g., first integrate over x using `HCubature.jl`, thus obtaining g(\tau,\tau') \equiv \int\_a^b \mathrm{d} x f(\tau,\tau',x), and then integrate over g using `GrundmannMoeller.jl` (or the other way around). The problem is that neither `integrate` from `GrundmannMoeller.jl` nor `hquadrature` from `HCubtarue.jl` returns functions, of course. They return numbers/vectors, which makes it tricky to perform subsequent integrations (unless I miss something).

---

<div class="post-metadata">

**Author:** ![Perrin\_Meyer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/perrin_meyer/32/15686_2.png) [@Perrin\_Meyer](https://discourse.julialang.org/u/Perrin_Meyer)\
**Post date:** [March 24, 2021, 5:08am UTC](https://discourse.julialang.org/t/integration-over-regions-with-varying-boundaries-triangular-domains-and-more/57157/11 "2021-03-24T05:08:06Z")

</div>

I’m an applied mathematician, who happens to be knee deep in trying to understand adaptive numerical quadrature routines for high order triangles for a boundary element code…

The BEM code uses Vioreanu Rokhlin triangles, and Koorwinder polynomials on said high order triangles for interpolation.

Perhaps you could use Koorwinder polynomials akin to how Trefethen et al use Chebyshev polynomials in Chebfun but for triangles, which would allow you to treat your functions more “symbolically” …

[GitHub - JuliaApproximation/ApproxFun.jl: Julia package for function approximation](https://github.com/JuliaApproximation/ApproxFun.jl)

More refs with possible insight…

SIAM J. SCI. COMPUT. c 2014  
Vol. 36, No. 1, pp. 267-288  
SPECTRA OF MULTIPLICATION OPERATORS AS A NUMERICAL TOOL  
B. VIOREANU AND V. ROKHLIN

BEM Code that has good discussion of Koorwinder polynomials and interpolation on triangles  
Fast multipole methods for evaluation of layer potentials with locally-corrected quadratures

> **[Fast multipole methods for evaluation of layer potentials with...](https://arxiv.org/abs/2006.02545v3)**
>
> While fast multipole methods (FMMs) are in widespread use for the rapid evaluation of potential fields governed by the Laplace, Helmholtz, Maxwell or Stokes equations, their coupling to high-order quadratures for evaluating layer potentials is still...