Integration over regions with varying boundaries, triangular domains and more

This topic may be considered, in some sense, as a follow-up of my previous thread. Thanks to stevengj and 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)
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”):

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

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];

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”. 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!

1 Like

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

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

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


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.

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

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


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.

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

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


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. :slight_smile:

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.

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

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

More refs with possible insight…

Vol. 36, No. 1, pp. 267-288

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