# Which method should I use for NIntegrate near a singularity?

I am looking at the following integral:

$$I= \int_{-\infty}^\infty d\tau_3 \int_{-\infty}^\infty d\tau_4 \frac{1}{(1+\tau_3^2)^2} \Phi \left(\frac{1+\tau_4^2}{1+\tau_3^2},\frac{(\tau_3-\tau_4)^2}{1+\tau_3^2} \right) \tag{1}$$

where $$\Phi(r,s)$$ is a complicated dimensionless function given in the code at the end of the question. Importantly, the integrand is finite except when $$\tau_3 = \tau_4$$ (see left plot below). By defining $$(\tau_3 - \tau_4)^2 > \epsilon^2$$ (point-splitting regularization), we can extract the divergence of $$(1)$$ to get:

$$\left. I \right|_\text{div} = - \frac{\pi^2}{2} \log \epsilon^2 \tag{2}$$

Now I would like to reproduce that result numerically, i.e. I redefine the integration limits as follows:

$$\left. I \right|_\text{reg} = \left(\int_{-\infty}^\infty d\tau_3 \int_{-\infty}^{\tau_3-\epsilon/2} d\tau_4 + \int_{-\infty}^\infty d\tau_3 \int_{\tau_3+\epsilon/2}^\infty d\tau_4 \right) \frac{1}{1+\tau_3^2} \Phi \left(\frac{1+\tau_4^2}{1+\tau_3^2},\frac{(\tau_3-\tau_4)^2}{1+\tau_3^2} \right) \tag{3}$$

My idea is the following: use NIntegrate to collect numerical data for different values of $$\epsilon$$ near $$0$$, then fit the data to a model $$I(\epsilon) = a \cdot \log \epsilon^2 + b$$. I have imagined this method myself, and thus can not back it up with references. Neither can I say that this should work for sure. However comparing the 3D plot of the integrand of $$(1)$$ and a plot for my fit function ($$a=-\pi^2/2, b=0$$) gives me hope:

The problem comes by the practical part. I thought LocalAdaptive would be a good method, but honestly that is just because of the name, I have no clue how it is evaluating the integral really. Anyhow I get the following data for $$\epsilon = 0.001, 0.002, ..., 0.010$$:

Not only I did not manage to reproduce the $$\log$$ behavior, but also when I increase WorkingPrecision, I find that NIntegrate fails to converge, which may suggest that the values I find are too inaccurate near the singularity. And if I include $$\epsilon=0$$ in the dataset, the integral converges, which obviously should not be the case, or at least it should have a large value compared to the rest of the dataset (I get $$32.4685$$, in the linear continuity of the plot).

So all in all the question is: how should I perform this numerical integration to obtain reliable data, and is this approach worth pursuing in the first place?

Here is my code so far:

x1 = 1;
R[\[Tau]3_, \[Tau]4_] := (x1^2 + \[Tau]4^2)/(x1^2 + \[Tau]3^2);
S[\[Tau]3_, \[Tau]4_] := (\[Tau]3 - \[Tau]4)^2/(x1^2 + \[Tau]3^2);
a[\[Tau]3_, \[Tau]4_] := 1/4 Sqrt[4*R[\[Tau]3, \[Tau]4]*S[\[Tau]3, \[Tau]4] - (1 - R[\[Tau]3, \[Tau]4] - S[\[Tau]3, \[Tau]4])^2];
F[\[Tau]3_, \[Tau]4_] := I Sqrt[-((1 - R[\[Tau]3, \[Tau]4] - S[\[Tau]3, \[Tau]4] - 4 I*a[\[Tau]3, \[Tau]4])/(1 - R[\[Tau]3, \[Tau]4] - S[\[Tau]3, \[Tau]4] + 4 I*a[\[Tau]3, \[Tau]4]))];
Phi[\[Tau]3_, \[Tau]4_] := 1/a[\[Tau]3, \[Tau]4] Im[PolyLog[2, F[\[Tau]3, \[Tau]4] Sqrt[R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]] + Log[Sqrt[R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]]*Log[1 - F[\[Tau]3, \[Tau]4] Sqrt[R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]]];
Integrand[\[Tau]3_, \[Tau]4_] := 1/(x1^2 + \[Tau]3^2)^2 Phi[\[Tau]3, \[Tau]4];
data = Table[{\[Epsilon], NIntegrate[Integrand[\[Tau]3, \[Tau]4], {\[Tau]3, -\[Infinity], \\[Infinity]}, {\[Tau]4, -\[Infinity], \[Tau]3 - \[Epsilon]/2},Method -> "LocalAdaptive"] + NIntegrate[Integrand[\[Tau]3, \[Tau]4], {\[Tau]3, -\[Infinity], \\[Infinity]}, {\[Tau]4, \[Tau]3 + \[Epsilon]/2, \[Infinity]}, Method -> "LocalAdaptive"]}, {\[Epsilon], 0.001, 0.01, 0.001}];
ListPlot[data]

• It is exactly the place where science transforms into art. There is a tutorial at Menu/Help/WolframDocumentation/tutorial/NIntegrateOverview/NIntegrateIntroduction/Automatic Singularity Handling and the next section "Special Strategies". Also two next parts of the /tutorial/NIntegrateOverview: "NIntegrate Integration Strategies" and "NIntegrate Integration Rules" may be useful. There are several strategies of how to cope with the singularities in the multidimensional case. One should try them one by one. Nobody can say which one is better. Commented Mar 26, 2020 at 11:23
• Continuation: according to what you write, I would first try the regularization that you have already proposed with a sequence of decreasing epsilons, and then check if this would converge. There has recently been a useful discussion on this subject here: mathematica.stackexchange.com/a/215891/4999 Commented Mar 26, 2020 at 11:27
• @AlexeiBoulbitch Thanks for all the useful links, I will go through them now! I am not sure I understand the suggestion that you make in your second comment: when I decrease the epsilons, the integral does converge, even at $0$, but I would like it to diverge. Or what do you mean?
– Pxx
Commented Mar 26, 2020 at 11:36
• Please have a look at my answer. Commented Mar 26, 2020 at 12:12

Generally, LocalAdaptive is considered to be less good as GlobalAdaptive. Try these approaches:

NIntegrate[
Integrand[τ3, τ4], {τ3, -100,
100}, {τ4, -100, τ3 - ϵ/2},
"SingularityHandler" -> "DuffyCoordinates"}, AccuracyGoal -> 3,
WorkingPrecision -> 10] // Timing

NIntegrate[
Integrand[τ3, τ4], {τ3, -100,
100}, {τ4, -100, τ3 - ϵ/2},
Method -> {"GlobalAdaptive", "SingularityHandler" -> "IMT"},
AccuracyGoal -> 3, WorkingPrecision -> 10] // Timing


yielding

(*
{1.85938, 15.74479851}
{1.65625, 15.74484120}
*)


The first figure is the time of computation and the second is the value. We see that the estimates of the integral are close to one another. The timing seems to be a bit better with IMT. The messages that you get on the way only indicate that the convergence is slow. They do not warn about any incorrectness of the calculation. Your condsturction:

data = Table[{ϵ = 10^-n,
NIntegrate[
Integrand[τ3, τ4], {τ3, -100,
100}, {τ4, -100, τ3 - ϵ/2},
"SingularityHandler" -> "DuffyCoordinates"},
AccuracyGoal -> 3, WorkingPrecision -> 10] +
NIntegrate[
Integrand[τ3, τ4], {τ3, -100,
100}, {τ4, τ3 + ϵ/2, 100},
"SingularityHandler" -> "DuffyCoordinates"},
AccuracyGoal -> 3, WorkingPrecision -> 10]}, {n, 2, 8}] //
Quiet;


gives

Show[{
ListLogLinearPlot[data, PlotRange -> All,
AxesLabel -> {Style["ϵ", 16, Black],
Style["int", 16, Black]}],
ListLogLinearPlot[data, PlotRange -> All, Joined -> True]
}]


looks as if the result converges to 31.68 or so.

Have fun!

• Hi, and thanks for your answer! Shouldn't the integral diverge as $\epsilon \to 0$? That is the result that I found in $(2)$ and would like to reproduce. Why does it appear finite numerically?
– Pxx
Commented Mar 26, 2020 at 12:15
• The numeric estimate does not show it. If you are right and it, indeed, diverges, this means that you cannot calculate it. No way. The SingularityHandler supresses the singularity during the numeric estimate of the integral, and I guess the weak singularity, such as the logarithmic one, can be erroneously removed. In this case you may obtain the value of the integral when in reality it does not exist. So, be careful. Commented Mar 26, 2020 at 12:19
• Hmm well, the idea was to reproduce the $-\log \epsilon^2$ divergence in the plot that you produced. What about disabling the SingularityHandler?
– Pxx
Commented Mar 26, 2020 at 12:29
• Try to put "SingularityHandler" -> None. Commented Mar 26, 2020 at 12:42
• @AlexeiBoulbitch Try to put "SingularityHandler" -> None. -- That one is my second favorite singularity handler! :) Commented Mar 26, 2020 at 12:59

When simplifying the functions and ComplexExpand the integrand, there are no problems with standard NIntegrate(besides converging slowly).

x1 = 1;
R[\[Tau]3_, \[Tau]4_] = (x1^2 + \[Tau]4^2)/(x1^2 + \[Tau]3^2);
S[\[Tau]3_, \[Tau]4_] = (\[Tau]3 - \[Tau]4)^2/(x1^2 + \[Tau]3^2);
a[\[Tau]3_, \[Tau]4_] =
1/4 Sqrt[4*R[\[Tau]3, \[Tau]4]*
S[\[Tau]3, \[Tau]4] - (1 - R[\[Tau]3, \[Tau]4] -
S[\[Tau]3, \[Tau]4])^2] //
FullSimplify[#, \[Tau]3 \[Element] Reals && \[Tau]4 \[Element]
Reals] &;

F[\[Tau]3_, \[Tau]4_] =
I Sqrt[-((1 - R[\[Tau]3, \[Tau]4] - S[\[Tau]3, \[Tau]4] -
4 I*a[\[Tau]3, \[Tau]4])/(1 - R[\[Tau]3, \[Tau]4] -
S[\[Tau]3, \[Tau]4] + 4 I*a[\[Tau]3, \[Tau]4]))] //
FullSimplify[#, \[Tau]3 \[Element] Reals && \[Tau]4 \[Element]
Reals] &;

Phi[\[Tau]3_, \[Tau]4_] =
1/a[\[Tau]3, \[Tau]4] Im[
PolyLog[2,
F[\[Tau]3, \[Tau]4] Sqrt[
R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]] +
Log[Sqrt[R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]]*
Log[1 - F[\[Tau]3, \[Tau]4] Sqrt[
R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]]] //
FullSimplify[#, \[Tau]3 \[Element] Reals && \[Tau]4 \[Element]
Reals] &;


.

Integrand[\[Tau]3_, \[Tau]4_] =
1/(x1^2 + \[Tau]3^2)^2 Phi[\[Tau]3, \[Tau]4] //
FullSimplify[#, \[Tau]3 \[Element] Reals && \[Tau]4 \[Element]
Reals] & // ComplexExpand[#, TargetFunctions -> {Re, Im}] & //
Simplify[#, \[Tau]3 \[Element] Reals && \[Tau]4 \[Element] Reals] &

(*   (1/((1 + \[Tau]3^2) Abs[\[Tau]3 - \[Tau]4]))(2 Im[
PolyLog[2,
I Sqrt[((1 + \[Tau]4^2) (-1 + (2 I)/(
I + ((\[Tau]3 - \[Tau]4) \[Tau]4)/
Abs[\[Tau]3 - \[Tau]4])))/(\[Tau]3 - \[Tau]4)^2]]] +
ArcTan[Abs[\[Tau]3 - \[Tau]4] +
Sqrt[1 + \[Tau]4^2]
Sin[1/2 ArcTan[-(-1 + \[Tau]4^2) Abs[\[Tau]3 - \[Tau]4],
2 (\[Tau]3 - \[Tau]4) \[Tau]4]], -Sqrt[1 + \[Tau]4^2] Cos[
1/2 ArcTan[-(-1 + \[Tau]4^2) Abs[\[Tau]3 - \[Tau]4],
2 (\[Tau]3 - \[Tau]4) \[Tau]4]]] Log[(
1 + \[Tau]4^2)/(\[Tau]3 - \[Tau]4)^2])   *)


Standard integration an with higher accuracy

NIntegrate[
Integrand[\[Tau]3, \[Tau]4], {\[Tau]3, -\[Infinity], \[Infinity]}, {\
\[Tau]4, -\[Infinity], \[Infinity]}]

(*   32.4697   *)

(nint = NIntegrate[
Integrand[\[Tau]3, \[Tau]4], {\[Tau]3, -\[Infinity], \
\[Infinity]}, {\[Tau]4, -\[Infinity], \[Infinity]},
WorkingPrecision -> 25, AccuracyGoal -> 6,
PrecisionGoal -> 6]) // Timing

(*   {63.125, 32.46969700779309434717063}   *)

• So it seems that the integral is finite, as the plots from Alexei and myself show, as well as your analysis, I will look again at my function and see if I have not made a mistake somewhere, as I expect a divergence there.
– Pxx
Commented Mar 26, 2020 at 12:58
• I am still not totally convinced that the integral converges, but I have made sure there is no mistake in my typing of the functions. I have opened a new post on the math stackexchange to clear the issue about the divergence before I work further on the numerical side.
– Pxx
Commented Mar 26, 2020 at 22:31