2D integration over non-rectangular domain using cubature


Hello everyone,

I have a very basic question.
I’m looking for the most efficient way to perform the following very classical integral:

For now I’m using hquadrature of the cubature package in two steps.

This works, but I’m not sure whether this is the most efficient way.
More generally, I would like to know, how to perform multidimensional integration over non-rectangular domains.
It seems cubature does not allow bounds to be functions. Am I correct?

Many thanks in advance


You can use a change of variables in general. However, since your problem is separable in x and u, you should take advantage of that, eg approximate the integral of g(u)du from 0 to x as some G(x), similarly for f, and then F(a)G(x) is your answer.

(BTW, does the forum allow LaTeX now? or was that inserted as an image?)


Thanks for your answer,
Ok, I will stick with two-step approach then.

BTW, yes, it’s an image. :slight_smile:


Why is the problem separable? The integration is over the x variable which is the limit of integration of the second integral.

In any case, your suggestion has been accepted.


The way I read this notation is that the dx closes the scope for the x, so the second x is a free variable. I am aware that some people follow a different convention, but I find that confusing.


It’s quite natural, once you get used to it, to think of as “\int dx” as the operator “integrate the following function with respect to x”


I would recommend trying a change of variables to turn the integration region into a rectangle. For example, let y=u/x, so that the integral becomes

In general, a 2d cubature code can be more efficient than nested 1d quadratures, assuming you are using adaptive quadrature/cubature. This is because the inner 1d integral, if it is performed as an independent 1d adaptive quadrature, can waste a lot of integrand evaluations trying to refine the 1d integral to high relative accuracy even if its overall contribution to the integral is small.

(There are also specialized cubature algorithms for triangular domains, but I’ve never seen an implementation of a nested/adaptive one.)


Thank you very much @stevengj for your answer.
That was exactly the information I was looking for.

I know in matlab you can actually use a function handle to define an integration region.
I was wondering if something similar was possible with one of the julia packages.

Thanks again


I’m not sure what you’re referring to — as far as I can tell, Matlab doesn’t include multi-dimensional integration routines.


This is what I found for Matlab R2017a


In older versions I used quad2d
In both versions you could pass a function handle (anonymous function) as an integration bound.

I quickly checked the publications cited at the end of the page, but I have no idea whether matlab actually implements the algorithms descibed in these or relies on something else.