enter image description hereI 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}];





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}]
  • $\begingroup$ 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. $\endgroup$
    – Michael E2
    May 13 at 0:24
  • $\begingroup$ 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 $\endgroup$ May 13 at 0:33
  • 1
    $\begingroup$ 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. $\endgroup$
    – Michael E2
    May 14 at 1:03
  • 1
    $\begingroup$ 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. $\endgroup$
    – Roman
    May 14 at 10:27
  • 1
    $\begingroup$ 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. $\endgroup$
    – Roman
    May 14 at 10:35


Your Answer

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