Integrate will not evaluate three-dimensional gaussian

Questions
, , , , , , ,
1

Hello, I have constructed a function which is simply a three-dimensional gaussian:

As you can see, these are both gaussian terms, so they should be easy to integrate. I have tried many methods of NIntegrate, which do not converge or give a consistent answer. This function is very simple and well-behaved, it is a ball at r=0 and a shell at r=3, both of gaussian shape. How can I produce the definite integral of this function?

Conv[x_, y_, z_, xp_, yp_, zp_] := (1/4)*1/Sqrt[Pi]*Exp[-((xp - x)^2 + (yp - y)^2 + (zp - z)^2)] + (3/4)*(1/Sqrt[Pi])*Exp[-4*(Sqrt[(xp - x)^2 + (yp - y)^2 + (zp - z)^2] - 3.5)^2]
2

Hello, I am not exactly sure what you mean by r in this case, or if you are integrating in a ball centered at {xp,yp,zp}.

However,

gauss = (1/4)*1/Sqrt[Pi]*
   Exp[-((xp - x)^2 + (yp - y)^2 + (zp - z)^2)] + (3/4)*(1/Sqrt[Pi])*
   Exp[-4*(Sqrt[(xp - x)^2 + (yp - y)^2 + (zp - z)^2] - 7/2)^2]

But, I am supposing now you are going to integrate in a ball centered at {xp,yp,zp}, so we can just move the ball to the origin:

gauss = Simplify[gauss /. {xp -> 0, yp -> 0, zp -> 0}]

Integrate doesn’t give up, but I didn’t wait for it to finish; so I gave it a little help:

gauss = Simplify[gauss /. x^2 -> r^2 - y^2 - z^2, 
  Assumptions :> r > 0]

This integrates quickly:

int = Integrate[4 Pi r^2 gauss, {r, 0, rBall}]

Plot[int, {rBall, 0, 3}]

gauss.png
gauss.png432Ă—279 7.59 KB

3

You have removed some information from the function in your steps. The integral is a sum of two functions. As I mentioned, if you set xp,yp,zp to zero, you get a ball at the origin, and a shell at 3.5. So if you transform to r, you will get a quick rise at 0, nothing in between, and another rise at 3-4.

This is a 3D integral, which has infinite bounds in both directions. The gaussians diverge so this is OK. I am trying to find the numerical value, symbolically if possible, of integrating over all space.

4

Hmmm,
Your original post’s formula wasn’t fully displayed, I thought I was able to copy it all.
I changed your 3.5 to 7/2. I didn’t see where you discussed {xp,yp,zp} or what you meant by r. You can let rBall->Infinity…

Does this not look correct:

gauss = (1/4)*1/Sqrt[Pi]*
   Exp[-((xp - x)^2 + (yp - y)^2 + (zp - z)^2)] + (3/4)*(1/Sqrt[Pi])*
   Exp[-4*(Sqrt[(xp - x)^2 + (yp - y)^2 + (zp - z)^2] - 7/2)^2]

equation.png

gauss = Simplify[gauss /. {xp -> 0, yp -> 0, zp -> 0}]
gauss = Simplify[gauss /. x^2 -> r^2 - y^2 - z^2, 
  Assumptions :> r > 0]
Plot[gauss, {r, 0, 4}]

4120gauss.png
4120gauss.png432Ă—274 8.83 KB

5

That does look like the correct integral. The simplification to zero I already did, which was the simplest way to evaluate the integral. I see the transformation to spherical coordinates. That works. I am surprised Mathematica cannot do these in cartesian fashion, though.

I did not know how to use the “Assumptions” option, thanks! Answer is…

1/32 ((42 Sqrt[\[Pi]])/E^49 + \[Pi] (305 + 297 Erf[7]))
6

Hmmm, I wouldn’t conclude that it can’t do the integral (I would be surprised if it couldn’t). It may just take longer than one has patience.

Limit[int, rBall -> Infinity]
(*gives your answer*)
7

I have waited a while without an answer. It would be very useful if I could do them in Cartesian so I could convolute them with other functions. I would like to integrate this with some cartesian exponentials that do not have any easy transformation to sphericals. I have a pretty fast computer, so I am not sure if Integrate will ever return an answer if I wait longer than I already have.

8

For example, I’d like to integrate the product of these two functions over prime coordinates:

\[Rho][xp_, yp_, zp_] := 
 1/(12*Pi^2)*Exp[-xp^2/16]*Exp[-yp^2/16]*Exp[-s*zp]*(1/(4*Pi))
Conv[x_, y_, z_, xp_, yp_, 
  zp_] := (1/4)*1/Sqrt[Pi]*
   Exp[-((xp - x)^2 + (yp - y)^2 + (zp - z)^2)] + (3/4)*(1/Sqrt[Pi])*
   Exp[-4*(Sqrt[(xp - x)^2 + (yp - y)^2 + (zp - z)^2] - (7/2))^2]

(with s some small real number, and 7/2 here being used for mu).

Transforming the density function to sphericals makes a huge mess. But this is just gaussian functions after the product, so I am surprised it doesn’t give definite integrals.

9

It will eventually return something, I’ve waited over a day in some cases. Definite integrals tend to take longer than indefinite ones. Using Assumptions often helps because it doesn’t need to find all the conditionals.

If I were doing your nonspherical convolutions, I’d probably go with a NIntegrate.