# Error in numerical integration using Wolfram Mathematica

I hope you are doing well.

I am having serious problems while calculating a numerical integral in Wolfram Mathematica, as the following error messages appear:

1- "Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small";

2- "NIntegrate failed to converge to prescribed accuracy after 27 recursive bisections in z near {[Rho],[Phi],z} = {0.935392,5.30379,-1.13233333152924497966296257800422608852386474609375000000000000000}.NIntegrate obtained 0.2509415135962599 and 0.0037048596190567802 for the integral and error estimates".

The idea of the code is: I define Q, the upper limit of integration (limit), parameter W and p=W-d. For each "d" within the range {0.001, W, 0.01}, a list of energies is generated for each value of "a" within the range {1/2, 2, 0.1}. That is, at the end of each "d", I have a list of "Energies" and then I store the minimum value of the list in a new list. In the end, I have a list of minimum energies for each d. The ListPlot is the graph that answers the following question: "For each d, what is the minimized energy value?"

Please find the code and plot below.

c2 = 2 \[Sqrt](E^(((2 (d + p))/
a))/(a^3 (18 a^4 (-E^(((2 d)/a)) - E^((2 p)/a) +
2 E^((2 (d + p))/a)) + 8 d^2 E^((2 (d + p))/a) p^2 -
2 a (d + p)^2 (d E^((2 p)/a) + E^((2 d)/a) p) -
3 a^3 (5 d E^((2 d)/a) + 7 d E^((2 p)/a) +
7 E^((2 d)/a) p + 5 E^((2 p)/a) p) +
2 a^2 (d^2 (-2 E^((2 d)/a) - 5 E^((2 p)/a) +
4 E^((2 (d + p))/a)) -
d (7 E^((2 d)/a) + 7 E^((2 p)/a) +
16 E^((2 (d + p))/a)) p + (-5 E^((2 d)/a) -
2 E^((2 p)/a) + 4 E^((2 (d + p))/a)) p^2)))) Sqrt[2/\[Pi]];

\[Psi][\[Rho]_, zz_, aa_, dd_, pp_] := c2*Exp[-Sqrt[\[Rho]^2 + zz^2]/aa]*(zz + dd)*(zz - pp);

Pot4[aj_, dj_, q_, lim_, pj_] := NIntegrate[-q*\[Rho]*\[Psi][\[Rho], z, aj, dj, pj]^2/(2*(z + pj)), {\[Rho], 0, lim}, {\[Phi], 0, 2*Pi}, {z, -dj, pj}];

Q=-0.5;

limit=1000;

W=7/3;

Emin={};

Do[list = {}; p = W - d; Do[
result = Pot4[a, d, Q, limit, p];
AppendTo[list, result] ,
{a, 1/2, 2, 0.1}];
AppendTo[Emin, {d, Min[list]}],
{d, 0.001, W, 0.1}]
ListLinePlot[Emin]

• One of the offending integrands seems to be -((-0.5 \[Rho] \[Psi][\[Rho], z, 0.5, 1.201, 1.1323333333333334]^2)/(2 (z + 1.1323333333333334))), which has a singularity in the integration region. Either that, or I mis-edited out the image link from the NIntegrate[] code. Please check. May 13 at 0:24
• Okay, I fixed what I wrote in the question. Thank you for the correction, it was really an error when I was adding the image May 13 at 0:33
• Your code worked (more) correctly before you accepted the edit by Reda.Kebbaj (I and another user had rejected the edit, but the OP can override any other user). Now there are new errors that make it difficult to work with. It would be good to rollback to the previous version. If you click on the "edited..." link above Reda.Kebbaj's user flair, you'll find "rollback" buttons under every version of the post. May 14 at 1:03
• Your integrand does not depend on $\phi$ and so you should perform the $\phi$-integration analytically by multiplying the integrand by $2\pi$ instead of integrating numerically. Also, the $\rho$-integration can be done analytically. This is not just a matter of speed but also of precision. May 14 at 10:27
• Also, why don't you integrate $\rho$ all the way to infinity? In Mathematica there is no problem with specifying Infinity` as an integration limit, and there is no need for artificial cutoffs. May 14 at 10:35