# 3D plot of a triple integral in polar coordinates

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]

• What is B in your Plot3D command? It isn't defined in your code snippet. Also, note that PolarAxes and PolarGridLines are not valid options for Plot3D. Jun 7 at 17:21
• Actually, wait; what do you mean by a "3D plot in polar coordinates"? Can you give us an example of such a plot? There's the SphericalPlot3D function, but you wouldn't be plotting your two angles from 0 to 2π if that's what you were trying to do. Jun 7 at 17:31
• Also note that your integral over k can be performed analytically. If you do that and feed the result into NIntegrate to do the u and c integrals, it'll save you some computer cycles. Jun 7 at 17:41
• If you implement all the hints given by @Michael and amend the AxesLabel -> {"a", "t", "R"}, you'll get (as I did) a nice 3D output Jun 7 at 18:10
• @Andreas: As did I. Takes a bit of time, though. Jun 7 at 18:12

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] • Thank you much, but it is interesting, because in my case the process starts and is going to long and does not finish. Now it took 15 minutes and is still going on. I am not sure that it will be finished at all.
– Hexe
Jun 7 at 21:14
• the integrand within NIntegrate doesn't match the defined function kintegr... Jun 7 at 21:18
• I wrote your code and obtained just an empty frame for a plot and no graph.
– Hexe
Jun 7 at 21:24
• takes about 200sec on Apple M1 Jun 7 at 21:33
• @Hexe: I had a small typo in the code which I've now fixed. Try running it again. Also note that you'll need to do run your own code first (everything up to the Plot3D command) before my code will work. Jun 7 at 22:22

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...

• Thank you very much. Now it works great ;)
– Hexe
Jun 8 at 7:16