5

I am trying to do integrals of the type $$ \int d^3\vec{p} \int d^3\vec{p}' e^{-p^2} e^{-{p'}^2}f(\vec{p}, \vec{p}') $$ where $\vec{p}$ and $\vec{p}'$ are three dimensional vectors represented using spherical coordinates, $\vec{p} = \{p,\theta,\phi\}$, and $f$ is a non-trivial, potentially complex, function. The integrals over $\phi$ and $\phi'$ can be done analytically even though the answers are rather complicated. However that is not true for the other integrals.

So I was wondering what would be the best method to approach this problem or if there are any packages (preferably for python) that do this kind of integrals. I plan to try SciPy's nquad but I hear that it is not suggested for integrals weighted by $e^{-p^2}$.

Daniel Shapero
  • 10,263
  • 1
  • 28
  • 59
e-eight
  • 163
  • 6
  • 1
    How cheap (relatively) is the evaluation of $f$? – Anton Menshov Aug 14 '18 at 18:54
  • Not quite sure, but $f$ is a product of rational functions of the $p$'s, Laguerre polynomials involving $p$'s and trig functions involving the angles. So I would think that evaluation is not too costly. However the function is quite oscillatory. – e-eight Aug 14 '18 at 19:03
  • I assume that your integral extend to infinity, am I right? – nicoguaro Aug 14 '18 at 19:48
  • The integral over $p$ and $p'$ are from 0 to $\infty$. The integral over $\theta$ and $\theta'$ are from 0 to $\pi$ and the integrals over $\phi$ and $\phi'$ are from 0 to $2\pi$. – e-eight Aug 14 '18 at 19:54
  • Then, I think that your integral is two dimensional, since the integrals over $\theta$ and $\phi$ are simple enough to do it analytically. – nicoguaro Aug 14 '18 at 22:43
  • 1
    Like I said the integrals over $\phi$ and $\phi'$ can be done analytically but the integrals over $\theta$ and $\theta'$ cannot be. The function is such that the variables cannot be seperated, for example $f$ can be $F(p,p')/(p^2 + p'^2 + 2pp'\cos\theta\cos\theta' + 2pp'\sin\theta\sin\theta'\cos(\phi-\phi') + a^2)$. – e-eight Aug 14 '18 at 23:10
  • 1
    You might find the answers to this question useful too. – Daniel Shapero Aug 15 '18 at 01:22
  • Could you provide a more general example of what $f$ looks like, for instance is it Laguerre polynomials multiplied by spherical harmonics of the same degree, etc.? Structure like this should inform your choice of integration rule. – smh Aug 16 '18 at 15:02
  • @smh It will be something like this $$L_n^{l+1/2}(p^2)L_{n'}^{l'+1/2}(p'^2)e^{-(p'2+p^2)/2}(pp')^{l+2}\sin\theta'\sin\theta \frac{1}{p^2p'^2+2pp'\cos\theta\cos\theta'+2pp'\sin\theta\sin\theta'\cos(\phi-\phi')+a^2}Y^*{l'm'}(\theta',\phi')Y{lm}(\theta, \phi)$$. – e-eight Aug 16 '18 at 18:56

1 Answers1

5

The Genz-Malik algorithm [1], as implemented in the cubature library, works well for computing 6-dimensional integrals.

[1] A. C. Genz and A. A. Malik, “Remarks on algorithm 006: An adaptive algorithm for numerical integration over an N-dimensional rectangular region,” Journal of Computational and Applied Mathematics, vol. 6, no. 4, pp. 295–302, Dec. 1980.

Juan M. Bello-Rivas
  • 3,994
  • 16
  • 24
  • Looks interesting. I will try it and see if it works. – e-eight Aug 15 '18 at 01:49
  • But, how could be used to integrate complicated boundaries? Is it possible? – Zythos Aug 18 '18 at 18:38
  • 1
    My understanding is that it can only work on hypercubes or regions that can be transformed into such. – Juan M. Bello-Rivas Aug 18 '18 at 19:01
  • 1
    @Zythos In most cases of interest you can convert the integrals to hypercubes by a change of variables. For example in the integral that I posted here, using the transformations $p = q/(1-q)$, $\theta = \pi t$, $\phi = 2\pi u$, and likewise for their primed counterparts will convert the integration region to a hypercube. Does not mean that the Genz-Malik will always work (as I am finding out), but it can still be applied. – e-eight Aug 20 '18 at 16:16
  • @monstergroup42 Yes, I have seen the overview in github. But if the boundary is only known numerically (a set of points that defines a closed volume) is there a way to do it? – Zythos Aug 20 '18 at 18:54
  • 1
    @Zythos If you have a mesh defining the surface, you can use a surface parametrization method (see https://doc.cgal.org/latest/Surface_mesh_parameterization/index.html#Surface_mesh_parameterizationFixedBorder ) to turn the mesh into a more manageable domain (like a circle that can, in turn, be transformed into a rectangle via polar coordinates). – Juan M. Bello-Rivas Aug 21 '18 at 01:56
  • @JuanM.Bello-Rivas Thanks for the reference! I am thinking of a particular problem that has surface points distributed non-uniformly so a mesh has difficulties representing correctly the true surface. Nevertheless, I will inspect your suggestion. I didn't knew that approach. – Zythos Aug 22 '18 at 06:24