Singularity in NIntegrate

I am trying to do integration with NIntegrate command. I know that the function I am integrating (the integrand) has a singularity at x=0. I am trying to exclude this point from the integration, but still, I am getting this error

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

This error causes me not to trust the final result of the integration. Also, it takes a long time to give the answer, maybe more than a second. Since I need to do this integral for thousands of points, one second for each point is too much. Can anyone help me to get rid of this error and make the integration faster?

The real code is a little complicated. Here is the simple version of the code:

a=4.1*10^6;
b=26;
d=8.25;
f[x_]=a/(x/b(1+x/b)^2)
NIntegrate[f[Sqrt[s^2+d^2-2*d*s*Cos[l]Cos[j]]],{s,0,d+100},{l,-Pi/9,Pi/9},{j,-Pi/9,Pi/9}]


Also, because of the nature of the problem I have, I can ignore the singularity part of my function and write it like this:

f[x_]=Piecewise[{{0, x < 1}, {a/(x/b (1 + x/b)^2), x >= 1}}]


This time not only do I get the previous error, but I also get a new one which is below. Also, the result takes more longer time.

 NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 4.906936408348816*^15 and 1.0014591983189883*^13 for the integral and error estimates.

• Let us consider the integrand a = 4.1*10^6; b = 26; d = 8.25; f[x_] = a/(x/b (1 + x/b)^2);f[Sqrt[s^2 + d^2 - 2*d*s*Cos[l] Cos[j]]], i.e. $$\frac{1.066\times 10^8}{\sqrt{-16.5 s \cos (j) \cos (l)+s^2+68.0625} \left(\frac{1}{26} \sqrt{-16.5 s \cos (j) \cos (l)+s^2+68.0625}+1\right)^2} .$$ If we put s s.t. Solve[(68.0625 + s^2)/(16.5 s ) == 1, s] which results in s -> 8.25 , then the integrand Sqrt[1-Cos[j]*Cos[l]] up to a constant multiplier has the non-integrable singularity at l==0&&j==0 so the integral under consideration diverges. Feb 28, 2022 at 8:15
• I was not right. The singularity 1/Sqrt[1-Cos[j]*Cos[l]] at l==0&&j==0 is integrable because of $$\frac{1}{\sqrt{1-\left(1-\frac{i^2}{2}\right) \left(1-\frac{l^2}{2}\right)}}\approx \frac {\sqrt2} {\sqrt{i^2+l^2}}$$ around the origin. @kcr, please undelete your answer. Feb 28, 2022 at 12:21
• @user64494 just saw this. Not sure why I did not get a notification, or maybe I missed it. I will add an edit to point out to your comments as they complete the answer in a nice analytic way. Thanks!
– user49048
Feb 28, 2022 at 17:37
• Note NIntegrate::slwcon is a warning about the behavior of the integration, not an error per se, but NIntegrate::eincr is an error with an estimate of the approximation error. However, the two together often means the error estimate is inaccurate. Setting MinRecursion -> 3 (or to 2 or 4) sometimes helps. If it does not reduce the error estimate much, then it's probably one of the problems indicated by slwcon (which is the case in this integral). Mar 1, 2022 at 0:00
• Thank you all for your time and helps. Your comments all are very informative and useful. Mar 3, 2022 at 22:30

If you try to execute

a = 4.1*10^6;
b = 26;
d = 8.25;
f[x_] = a/(x/b (1 + x/b)^2);
NIntegrate[
f[Sqrt[s^2 + d^2 - 2*d*s*Cos[l] Cos[j]]], {s, 0, d + 100}, {l, -Pi/9,
Pi/9}, {j, -Pi/9, Pi/9}, WorkingPrecision -> 10]


you get

2.030657119*10^8

without any errors. Does this resolve your issue? What I mean is if this works for the more complicated case of interest you mention in the OP.

Edit: if you try

NIntegrate[
Rationalize[f[Sqrt[s^2 + d^2 - 2*d*s*Cos[l] Cos[j]]]], {s, 0,
d + 100}, {l, -Pi/9, Pi/9}, {j, -Pi/9, Pi/9}, WorkingPrecision -> 5,
MaxRecursion -> 20]
NIntegrate[
Rationalize[f[Sqrt[s^2 + d^2 - 2*d*s*Cos[l] Cos[j]]]], {s, 0,
d + 100}, {l, -Pi/9, Pi/9}, {j, -Pi/9, Pi/9},
WorkingPrecision -> 10, MaxRecursion -> 20]
NIntegrate[
Rationalize[f[Sqrt[s^2 + d^2 - 2*d*s*Cos[l] Cos[j]]]], {s, 0,
d + 100}, {l, -Pi/9, Pi/9}, {j, -Pi/9, 0, Pi/9},
WorkingPrecision -> 15, MaxRecursion -> 20]


you get

2.0315*10^8

2.030657119*10^8

2.03068474164791*10^8

and an error from the final one. If you see the the documentation for the error I think that the result

$$\sim$$ 2.03*10^8

can be trusted, unless I am missing something.

Edit: make sure to check the comments under the OP by @user64494 for nice, analytic explanation of why the singularity is integrable.

• Let us consider the integrand a = 4.1*10^6; b = 26; d = 8.25; f[x_] = a/(x/b (1 + x/b)^2);f[Sqrt[s^2 + d^2 - 2*d*s*Cos[l] Cos[j]]], i.e. $\frac{1.066\times 10^8}{\sqrt{-16.5 s \cos (j) \cos (l)+s^2+68.0625} \left(\frac{1}{26} \sqrt{-16.5 s \cos (j) \cos (l)+s^2+68.0625}+1\right)^2}$. If we put s s.t. Solve[(68.0625 + s^2)/(16.5 s ) == 1, s] which results in s -> 8.25 , then the integrand Sqrt[1-Cos[j]*Cos[l]] up to a constant multiplier has the non-integrable singularity at l==0&&j==0 so the integral under consideration diverges. Feb 28, 2022 at 8:09
• @user64494 thanks for pointing this out. Totally missed this one. I will delete my answer shortly
– user49048
Feb 28, 2022 at 8:12
• +1. This is my vote. Feb 28, 2022 at 18:45
• @user64494 thanks again for the analysis and the comments!!!
– user49048
Feb 28, 2022 at 19:52
• Thank you all for your effort to help me. I am so appreciative. So, by increasing the recession and max recursion, the errors disappear in some points but make the code slower. So, I decided to keep the recession and max recursion the default value. I learned a lot of things in this discussion. Mar 3, 2022 at 22:32

If we identify the singularity, the default singularity handler handles it:

a = 4.1*10^6;
b = 26;
d = 8.25;
f[x_] = a/(x/b (1 + x/b)^2);
NIntegrate[f[Sqrt[s^2 + d^2 - 2*d*s*Cos[l] Cos[j]]],
{s, 0, d, d + 100}, {l, -Pi/9, 0, Pi/9}, {j, -Pi/9, 0, Pi/9}]

(*  2.03069*10^8  *)


For the other integrand, using the "LocalAdaptive" strategy helps, as does identifying the singularities, include the boundaries of the Piecewise elements (use Solve with rationalized parameters):

a = 41*10^5;
b = 26;
d = 825/100;
f[x_] = Piecewise[{{0, x < 1}, {a/(x/b (1 + x/b)^2), x >= 1}}];
NIntegrate[f[Sqrt[s^2 + d^2 - 2*d*s*Cos[l] Cos[j]]],
{l, -Pi/9, 0, Pi/9}, {j, -Pi/9, 0, Pi/9},
{s,  (* switch order of integration *)
0,
Piecewise[{{1/4 (33 Cos[j] Cos[l] -
2 Sqrt[-(1073/4) + 1089/4 Cos[j]^2 Cos[l]^2]), -(1073/4) +
1089/4 Cos[j]^2 Cos[l]^2 >= 0}}, d],
Piecewise[{{1/2 (33/2 Cos[j] Cos[l] +
Sqrt[-(1073/4) + 1089/4 Cos[j]^2 Cos[l]^2]), -(1073/4) +
1089/4 Cos[j]^2 Cos[l]^2 >= 0}}, d],
d + 100},

(*  {1.10917, 1.93676*10^8}  *)


It takes longer without the singularities identified, but it computes a similar answer.

NIntegrate[f[Sqrt[s^2 + d^2 - 2*d*s*Cos[l] Cos[j]]],
{l, -Pi/9, Pi/9}, {j, -Pi/9, Pi/9}, {s, 0, d + 100},

• (+1) from me. I just wanted to point out that I tried replacing LocalAdaptive with GlobalAdaptive since in the documentation it is explicitly stated "GlobalAdaptive" Is Generally Better than "LocalAdaptive". However, in this case GlobalAdaptive runs with errors. Quit interesting. Relevant link for the interested reader: reference.wolfram.com/language/tutorial/…
• @kcr Thanks. I don't grok "LocalAdaptive" as well as "GlobalAdaptive". I went the other way: When "GlobalAdaptive" proved to be a bear, it occurred to me to try "LocalAdaptive". Maybe it's the piecewise boundary from the integrand, x >= 1, which is not a straight. Maybe something else. I often forget to try LA, in fact, but it is usually worse. Mar 1, 2022 at 4:01
• @kcr Usually identifying the singularities solves the "GlobalAdaptive" problems (unless the integrand is highly oscillatory). I'm not sure what's the problem with the second integral. Mar 1, 2022 at 5:07