# Discretization and Dirac Delta function

Consider the following integral: \begin{align} I= \int_{\mathbb{R}^9}\frac{d^3\mathbf{x_1}\,d^3\mathbf{x_2}\,d^3\mathbf{x_3}}{|\mathbf{x}_1|\sqrt{|\mathbf{x}_2||\mathbf{x}_3|}}f(\mathbf{x_2})f(\mathbf{x_3})\delta(\Omega+|\mathbf{x}_1|-|\mathbf{x}_2|)\delta(\Omega+|\mathbf{x}_1|-|\mathbf{x}_3|) \end{align} where $$\mathbf{x_j} = (x_j,y_j,z_j)$$. This is a nine-dimensional integral. I would typically choose $$f(\mathbf{x})$$ to be some Gaussian of the form \begin{align} f(\mathbf{x}-\mathbf{x}_0) = \exp -\sum_{j=1}^3\frac{-(x_j-x_{0,j})^2}{2\sigma_j} \end{align} where $$\mathbf{x}_{0}$$ is some constant vector. Now, this integral is most likely not doable analytically, so I would need NIntegrate. NIntegrate does not understand DiracDelta[x] so the way I would try to do this is using $$\text{ImplicitRegion}\left[\Omega+|\mathbf{x}_1|-|\mathbf{x}_2| = 0\, \wedge \, \Omega+|\mathbf{x}_1|-|\mathbf{x}_3| = 0\right]$$.

Problem: it turns out that this does not work (i.e. Mathematica does not compute this), and the reason (as far as I understand) is this: at the backend, Mathematica needs to perform DiscretizeRegion which is only supported for three-dimensional domain (the error messages show this).

Question: is there a way out of this problem, i.e. perform this numerically without facing problem with DiscretizeRegion?

This is mainly frustrating because if $$\sigma_j=\sigma$$ for $$j=1,2,3$$, I could in fact do this integral by hand but for my problem at hand I need at least $$\sigma_3\neq \sigma_1,\sigma_2$$ which definitely makes it not analytically solvable.

Remark: I tried breaking the problem into a few steps, namely first doing the $$\mathbf{x_2}$$ integral (integral over $$\mathbf{x}_3$$ is identical) over region $$R$$ and then hoping that integrating over $$\mathbf{x}_1$$ will work. This means

R[y1_, y2_, y3_, s1_, s2_, s3_, gap_] :=
DiscretizeRegion[ImplicitRegion[{Sqrt[x1^2 + x2^2 + x3^2] - Sqrt[y1^2 + y2^2 + y3^2] == gap},
{{x1, -5 s1, 5 s1}, {x2, -5 s2, 5 s2}, {x3, -5 s3, 5 s3}}]]

where I chose $$s_j=\sigma_j$$ just to make the integral domain bounded (ideally I want $$x_j$$ to be integrated over the real line since the original integral is over $$\mathbb{R}^9$$ instead of $$\pm 5s_j$$). Then

G[y1_, y2_, y3_, gap_, x0_, y0_, z0_, s1_, s2_, s3_] := NIntegrate[f[x1, x2, x3, x0, y0, z0, s1, s2, s3]/Sqrt[Sqrt[x1^2 + x2^2 + x3^2]],
Element[{x1, x2, x3}, R[y1, y2, y3, s1, s2, s3, gap]]]

where $$G$$ essentially is the result after integrating over $$\mathbf{x}_2$$. The original integral is obtained by integrating $$G(\mathbf{x_1})^2/|\mathbf{x_1}|$$. Now, $$G(\mathbf{x}_1)$$ actually gives a definite result for fixed $$\mathbf{x}_1:=(y_1,y_2,y_3)$$ but I cannot integrate over $$y_j$$.

• I think the real problem here is that you need to be very careful with your deltas in higher dimensions. The delta you use in the top equation is 1D, since its arguments are scalar. I presume you want to integrate over a spherical shell. This is quite different from a nD delta that picks out a single point from an nD volume under an integral. I think it's easier to just try and define your integral over the surface you're trying to pick out. – Sjoerd Smit Jul 31 '19 at 7:17
• Comments are not for extended discussion; this conversation has been moved to chat. – Kuba Jul 31 '19 at 12:55

In spherical coordinates,

x1 = r1 {Sin[θ1] Cos[φ1], Sin[θ1] Sin[φ1], Cos[θ1]};
x2 = r2 {Sin[θ2] Cos[φ2], Sin[θ2] Sin[φ2], Cos[θ2]};
x3 = r3 {Sin[θ3] Cos[φ3], Sin[θ3] Sin[φ3], Cos[θ3]};

and with the definition (modified a bit from yours because I think you had some typos)

f[x_, y_, z_] = Exp[-((x - x0)^2/(2 σx^2))]*
Exp[-((y - y0)^2/(2 σy^2))]*
Exp[-((z - z0)^2/(2 σz^2))];

Integrate[((f @@ x2) (f @@ x3))/(r1 Sqrt[r2 r3]) *
r1 Sin[θ1] * r2 Sin[θ2] * r3 Sin[θ3] /. {r2 -> r1 + Ω, r3 -> r1 + Ω},
{r1, Max[0,-Ω], ∞}, {θ1, 0, π}, {φ1, 0, 2 π},
{θ2, 0, π}, {φ2, 0, 2 π},
{θ3, 0, π}, {φ3, 0, 2 π}]

In this way you don't have any Dirac distributions left and you have a regular integral. Further, as $$r_2=r_3\ge0$$, the denominator is cancelled by the Jacobian and the whole integral becomes non-singular.

To continue, here's the explicit numerical integral, where $$\theta_1$$ and $$\phi_1$$ have been integrated exactly (thanks @yarchik!), and we are left with five integration variables, four of which are finite:

4π*NIntegrate[
f[(r1+Ω)Sin[θ2]Cos[φ2], (r1+Ω)Sin[θ2]Sin[φ2], (r1+Ω)Cos[θ2]]*
f[(r1+Ω)Sin[θ3]Cos[φ3], (r1+Ω)Sin[θ3]Sin[φ3], (r1+Ω)Cos[θ3]]*
(r1+Ω)*Sin[θ2]*Sin[θ3],
{r1, Max[0, -Ω], ∞}, {θ2, 0, π}, {φ2, 0, 2π}, {θ3, 0, π}, {φ3, 0, 2π}]

In this integral there is no cross-talk between $$\vec{x}_2$$ and $$\vec{x}_3$$, so we can decompose it further, noticing that the integrals over $$\vec{x}_2$$ and $$\vec{x}_3$$ give the same result $$F(r_1+\Omega)$$:

F[r_?NumericQ] := NIntegrate[f[r*Sin[θ]Cos[φ], r*Sin[θ]Sin[φ], r*Cos[θ]]*Sin[θ], {θ,0,π}, {φ,0,2π}]

4π*NIntegrate[(r1+Ω)*F[r1+Ω]^2, {r1, Max[0, -Ω], ∞}]

In this form, the 9-dimensional integral has been reduced to three dimensions.

• Integral over $x_1$ can be performed beforehand. It reduces the whole thing to a 6d integral with $\delta(|\vec x_2|-|\vec x_3|)$. This one can be integrated one time in spherical coordinates. – yarchik Jul 31 '19 at 8:36
• @yarchik as $\vec x_2$ and $\vec x_3$ depend on $\vec x_1$ through the Dirac distribution, I don't see how you can do the $\vec x_1$ integral beforehand. Could you please elaborate? – Roman Jul 31 '19 at 9:59
• $I=\iint d\vec x_2d\vec x_3\int x_1 dx_1d\Omega_1 F(\vec x_2,\vec x_3)\delta(\Omega+x_1-x_2)\delta(\Omega+x_1-x_3) = 4\pi\times\iint d\vec x_2d\vec x_3 F(\vec x_2,\vec x_3)(x_2-\Omega)\delta(x_2-x_3)$. – yarchik Jul 31 '19 at 11:38
• Nice, you managed to go even one step further ! – yarchik Jul 31 '19 at 19:40
• Thanks Roman and @yarchik, I will have to convince myself for a while in a bit. While I am at it, I should check something: in my implementation above, is ImplicitRegion I imposed the same as using Dirac delta over the hypersurface? You both removed Dirac delta in the intermediate step but I was wondering if brute force way like mine would work in principle? – Everiana Aug 1 '19 at 0:12