3D plot of a triple integral in polar coordinates

Question

I have a triple integral with 2 changing angles (a and t). I need to integrate it and obtain a 3DPlot in polar coordinates (R,a,t). But Mathematica does not make it. What is wrong here?

s = 3;
n = 1;
p = 1;
r = 1;
L = Sqrt[(4*p + r^2)/3];
b = Sqrt[2*n*p];
R[a_, t_] := NIntegrate[((k^2)*Exp[-1.5*k^2])*((u^2)*(1 - u^2)*Exp[-(b*u)^2]*
  (Cos[s*k*u*Cos[a]/L]))*((Cos[c])^2)*(Cos[k*Sqrt[1 - u^2]*
  (s*Sin[a]*Cos[t]*Cos[c] + s*Sin[a]*Sin[t]*Sin[c])/L]), 
  {k, 0, Infinity}, {u, -1, 1}, {c, 0, 2*Pi}]

Plot3D[R[a, t]*B, {a, 0, 2*Pi}, {t, 0, 2*Pi}, PolarAxes -> True, 
  PolarGridLines -> Automatic, PlotRange -> All, AxesLabel -> {"a", "R", "t"}, 
  ImageSize -> 500]

Answer 1

A few notes:

• In your triple integral, the integral over k can be solved analytically. Note that your integrand is basically of the form $$ X e^{-3k^2/2} k^2 \cos (Y k) \cos (Z k) = \frac{X}{2} e^{-3k^2/2} k^2 \left[ \cos ((Y-Z) k) + \cos (Y+Z) k \right] $$ where $X, Y, Z$ are quantities that are independent of $k$. Mathematica can do the integral of this quantity from $k = 0$ to $\infty$; it takes about half a minute but does arrive at a result:$$ \frac{1}{18} \sqrt{\frac{\pi }{6}} e^{-\frac{1}{6} (Y+Z)^2} \left(-(Y+Z)^2-e^{\frac{2 Y Z}{3}} \left((Y-Z)^2-3\right)+3\right) $$ Doing this integral analytically allows you to reduce your numerical integral from three dimensions down to two, which leads to a huge speedup in the amount of time required to calculate each point.

• In your original code, B is not defined. You'll need to actually define it to get a plot; the argument to plotting functions must always return a number (not a symbolic quantity) when it's evaluated at points in the plotting range. In what's below, I have simply set B = 1, but you can define it as you wish.

• PolarAxes and PolarGridLines are not valid options for Plot3D. To create a spherical plot in 3D, you need to use the function SphericalPlot3D instead. If a is your "polar" angle and t is your "azimuthal" angle, then you want to issue the command SphericalPlot3D[(*function*), a, t]

Revised code

The code below takes about three minutes to yield the plot at the end.

kintegr[u_, c_, a_, t_] = 
          Integrate[((k^2)*Exp[-3*k^2/2])*((u^2)*(1 - u^2)*
          Exp[-(b*u)^2]*(Cos[s*k*u*Cos[a]/L]))*((Cos[c])^2)
          *(Cos [k*Sqrt[1 - u^2]*(s*Sin[a]*Cos[t]*Cos[c]+s*Sin[a]*Sin[t]*Sin[c])/L]), {k, 0, Infinity}, 
          Assumptions -> {-1 <= u <= 1, 0 <= c <= 2 Pi, 0 <= a <= Pi, 0 <= t <= 2 Pi}]
R2[a_?NumericQ, t_?NumericQ] := 
          NIntegrate[kintegr[u, c, a, t], {u, -1, 1}, {c, 0, 2*Pi}]
SphericalPlot3D[R2[a, t], a, t]

Answer 2

Comments break the code formatting. So I put the code into this answer:

s = 3; n = 1; p = 1; r = 1; L = Sqrt[(4*p + r^2)/3]; b = Sqrt[2*n*p];
kintegr[u_, c_, a_, t_] = 
Integrate[((k^2)*Exp[-3*k^2/2])*((u^2)*(1 - u^2)*
Exp[-(b*u)^2]*(Cos[s*k*u*Cos[a]/L]))*((Cos[c])^2)*
(Cos[k*Sqrt[1 - u^2]*(s*Sin[a]*Cos[t]*Cos[c] + s*Sin[a]*Sin[t]*Sin[c])/L]),
{k, 0, Infinity}, 
Assumptions -> {-1 <= u <= 1, 0 <= c <= 2 Pi, 0 <= a <= Pi, 0 <= t <= 2 Pi}]
R2[a_?NumericQ, t_?NumericQ] := 
NIntegrate[kintegr[u, c, a, t], {u, -1, 1}, {c, 0, 2*Pi}]
SphericalPlot3D[R2[a, t] 10., a, t]

with B set to 10. The credit goes to Michael...