nquad: Best way to use quadpy to compute matrix integral on simplex

[DanPuzzuoli]

Hi,

I'd like to write a function that: given an integer k >= 1, an array-valued function f on k real variables, and a float T > 0, computes the nested integral of f(t1, t2, ..., tk) for:

• t1 from 0 to T,
• t2 from 0 to t1,
• ...
• tk from 0 to t{k-1}

I believe this domain of integration is a simplex (hence the title).

As an initial attempt, I've been using scipys nquad to evaluate these integrals term by term, using a change of variables to transform the integral into one for which all variables are integrated from 0 to T. Evaluating things term by term however is very slow.

I am wondering: what is the most convenient tool in quadpy for performing these integrals (preferably something adaptive)? I see that there are methods labeled for n-simplex and n-cube (the latter of which would be useable using the same change of variables I used with nquad), but I'm unsure of what will be the most natural and convenient.

[nschloe]

This indeed is a simplex with corners

[0, 0, ...]
[T, 0, ...]
[T, T, 0, ...]
...
[T, T, T, ...]

quadpy does have a few integration schemes for n-simplices, https://github.com/nschloe/quadpy#n-simplex-tn.

Note that you don't get adaptivity here. Adaptivity in quadpy's routines is based on the assumption that you can subdivide an integration domain into multiple smaller blocks. You can do that for intervals, and rectangles, triangles. It is, however, not possible to canonically subdivide a tetrahedron into a number of smaller tetrahedra the same way, and much less for for general n-simplices. (I say canonically here because you can of course split a tetrahedron into many small tetrahedra, but not in an obvious way.)

[DanPuzzuoli]

Ah I see, interesting, thank you.

Are there adaptive methods for n-cubes?

[nschloe]

Not right now, but this should actually be in quadpy. I'll take a look.