I am working on a problem (solid state physics, I am stripping all the details for brevity but if more details can help I'll elaborate) where I need to numerically calculate an integral of the form: $$\vec{A}(\vec{r}) = \int_{\mathbb{R}^{3}\setminus E}\vec{B}(\lambda, \vec{r},\vec{w}) d^{3}w, \quad \vec{r} \in E $$

where E is volume of rotational ellipsoid centered in the origin with z-axis as its axis of symmetry. My questions are:

a) If the function $\vec{B}$ is well behaved (finite, no sharp peaks or wild oscillations, decays exponentially far from the origin so convergence should not be a problem) should one try to integrate this with Monte Carlo or classical integration methods? I am considering the cuba4.2 library for the implementation

b) Would you expect any problems when integrating over domain without some interior, line-connected volume? The implementation would be the obvious one - integrand function would simply return 0.0 whenever $\vec{w}$ is within $E$. Cuba allows for some preexisting knowledge to be supplied (like location of peaks) but not flexible boundaries.

c) The real problem here is amount of work. Basically, $\lambda$ describes parameter space which I need to parse fairly densely - $10^{4}$-$10^{5}$ points and for each $\lambda$ I need to calculate integral for $10^{5}$ different $\vec{r}$, so considerable amount of work - preliminary runs give estimate around 45 days straight on my PC. Cuba itself is written for parallel execution, but if anyone has used cuba library - are there room for any speedups, such as using "state" files to provide an preliminary mash that can be further refined or something like that? Idea would be, since I got to do $10^5$ integrals for each $\lambda$, a reasonable starting mash could help achieve the precision goal faster. Any other advice?

  • 1
    $\begingroup$ Can you find the integral of $B$ over $\mathbb{R}^3$? $\endgroup$
    – nicoguaro
    Commented Dec 10, 2017 at 0:11
  • $\begingroup$ Can you use the equations that describe $\vec A$ and $\vec B$ to express the integral via integrals over the boundary of $E$? This would reduce the awkwardness of the infinite integration domain, and bring things down to a lower-dimensional manifold. $\endgroup$ Commented Dec 10, 2017 at 6:57
  • $\begingroup$ @nicoguaro: $\overrightarrow{B}$ is defined everywhere but within E has nonintegrable poles. $\endgroup$ Commented Dec 10, 2017 at 9:33
  • 1
    $\begingroup$ If you return 0 inside $E$ that kills most classical integration methods, because they need to assume the function is sufficiently differentiable, and you'd introduce a discontinuity. If you change variables to elliptical coordinates, it could change the domain to something regular that a normal method could handle. Also, classical methods handle fast decaying functions by treating them as integrals with weights, which is more efficient. $\endgroup$
    – Kirill
    Commented Dec 10, 2017 at 17:09
  • 1
    $\begingroup$ I'd still try to rewrite things as a boundary integral. This takes away one source of error, and the dimensional reduction should yield a much more accurate result. $\endgroup$ Commented Dec 10, 2017 at 17:20

1 Answer 1


Ok, just a small followup - maybe someone would find it useful. So I took a test volume to integrate which is, in prolate spheroidal coordinate system, per @Kirill suggestion

$\frac{z^{2}}{\mu ^{2}}+\frac{\rho^{2}}{\mu ^{2}-1} = L^{2} $

$\frac{z^{2}}{\nu ^{2}}+\frac{\rho^{2}}{1-\nu ^{2}} = L^{2} $

$L = 1$

bounded with two curves, namely

$2 \leqslant \mu \leqslant 3$

Curves bounding calculation volume

So, the volume integral is:

$V = \int_{0}^{2\pi}\int_{-1}^{1}\int_{2}^{3} (\mu ^{2}-\nu ^{2})d\mu d\nu d\phi$

with exact value of $V = 24\pi \approx 75.398223686 $. When translated and rescaled into unit hypercube as required by cube4.2, the integral to be calculated is:

$V = 4\pi \int_{0}^{1}\int_{0}^{1}\int_{0}^{1} (p-2q+3)(p+2q+1)dp dq dw$

Results of all algorithms except Suave (cant get it to run for some reason) of cuba4.2 is:

-------------------- Vegas test -------------------
VEGAS RESULT:   neval 4500      fail 0
VEGAS RESULT:   75.40052478 +- 0.06988593       p = 0.000

-------------------- Divonne test -------------------
DIVONNE RESULT: neval 2134      fail 0  nregions 14
DIVONNE RESULT: 75.40013735 +- 0.03957392       p = 0.000

-------------------- Cuhre test -------------------
CUHRE RESULT:   neval 381       fail 0  nregions 2
CUHRE RESULT:   75.39816000 +- 0.00000000       p = 0.099

My original proposal was trying to integrate the characteristic function $\chi $ of the subset for which @Kirill suggested would break classical integration methods. Here, we'd have:

$V = 6^{3} \int_{0}^{1}\int_{0}^{1}\int_{0}^{1} \chi dx dy dz$

$\chi = \left\{\begin{matrix} & (2z-1)^{2}+1.125((2x-1)^{2}+(2y-1)^{2}) \leqslant 1\\ 1 & \wedge\\ & 2.25(2z-1)^{2}+3((2x-1)^{2}+(2y-1)^{2}) > 1\\ 0 & \end{matrix}\right.$

Numerical integration results are:

-------------------- Vegas test -------------------
VEGAS RESULT:   neval 52000     fail 1
VEGAS RESULT:   75.47990417 +- 0.36806625       p = 0.005

-------------------- Divonne test -------------------
DIVONNE RESULT: neval 31419     fail 0  nregions 24
DIVONNE RESULT: 75.24224767 +- 0.07119902       p = 0.915

-------------------- Cuhre test -------------------
CUHRE RESULT:   neval 50165     fail 1  nregions 198
CUHRE RESULT:   75.54317090 +- 19.77503339      p = 0.000

Ok, so in conclusion I can say two things - expressing integral function in prolate system makes considerable difference in precision and speed of calculation - entire order of magnitude. Trying to integrate characteristic function of a set doesn't break classical function in a sense that they can't deliver, it breaks them in as they need to do 100x evaluation more to get a much more imprecise result. Divonne has best performance in characteristic function integration, as it seems its importance mash adjust to borders way better than VEGAS.

Here are speed tests made in python 3.6 with pycuba 4.2

Characteristic function integration:

  • VEGAS: 237 ms ± 5.28 ms per loop
  • DIVONNE: 145 ms ± 933 µs per loop
  • CUHRE: 214 ms ± 2.87 ms per loop

Prolate coordianate system:

  • VEGAS: 37.6 ms ± 217 µs per loop
  • DIVONNE: 29.6 ms ± 266 µs per loop
  • CUHRE: 24.5 ms ± 260 µs per loop
  • $\begingroup$ @WolfgangBangerth Actually that was quite a neat idea with dimension reduction, if I manage to work something out, I'll post. $\endgroup$ Commented Dec 29, 2017 at 17:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.