Message Boards Message Boards

Integrate will not evaluate three-dimensional gaussian

Posted 11 years ago
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]
POSTED BY: C. E. Coppola
8 Replies
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}]
POSTED BY: W. Craig Carter
Posted 11 years ago
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.
POSTED BY: C. E. Coppola
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]
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}]
POSTED BY: W. Craig Carter
Posted 11 years ago
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]))
POSTED BY: C. E. Coppola
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*)
POSTED BY: W. Craig Carter
Posted 11 years ago
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.
POSTED BY: C. E. Coppola
Posted 11 years ago
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.
POSTED BY: C. E. Coppola
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.
POSTED BY: W. Craig Carter
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract