1
$\begingroup$

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]
$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$ Jun 7, 2023 at 17:21
  • $\begingroup$ 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. $\endgroup$ Jun 7, 2023 at 17:31
  • $\begingroup$ 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. $\endgroup$ Jun 7, 2023 at 17:41
  • $\begingroup$ 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 $\endgroup$
    – Andreas
    Jun 7, 2023 at 18:10
  • $\begingroup$ @Andreas: As did I. Takes a bit of time, though. $\endgroup$ Jun 7, 2023 at 18:12

2 Answers 2

4
$\begingroup$

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]

enter image description here

$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – Hexe
    Jun 7, 2023 at 21:14
  • $\begingroup$ the integrand within NIntegrate doesn't match the defined function kintegr... $\endgroup$
    – Andreas
    Jun 7, 2023 at 21:18
  • $\begingroup$ I wrote your code and obtained just an empty frame for a plot and no graph. $\endgroup$
    – Hexe
    Jun 7, 2023 at 21:24
  • $\begingroup$ takes about 200sec on Apple M1 $\endgroup$
    – Andreas
    Jun 7, 2023 at 21:33
  • $\begingroup$ @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. $\endgroup$ Jun 7, 2023 at 22:22
1
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ Thank you very much. Now it works great ;) $\endgroup$
    – Hexe
    Jun 8, 2023 at 7:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.