# Multidimensional NIntegral with singularity

I am dealing with an integral

$F(k)=-\int d^2p d^2q\frac1{\lvert p\rvert}\frac1{\lvert p-q\rvert}(\frac1{\lvert k-q\rvert}-\frac1{\lvert q\rvert})$

when I am evaluating it at some point, say k=30

NIntegrate[-p q/p/Sqrt[p^2 + q^2 - 2 p q Cos[\[Phi]p - \[Phi]q]](1/
Sqrt[30^2 + q^2 - 60 q Cos[\[Phi]q]]-1/q), {p, 0, 50}, {q, 0,
50}, {\[Phi]p, 0, 2 \[Pi]}, {\[Phi]q, 0, 2 \[Pi]}]


MMA gives 1397.63 with warning

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.

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
integration. NIntegrate obtained -3120.06+0.187845 I and 1266.8675450228204
for the integral and error estimates.


The integtal has singularity but should be converget. I tried to increase Maxrecursion Minrecursion or Working precision does not help.

Besides, if I split the integral in to 2 parts from bracket of the formula and integrate separately

NIntegrate[-p q/p/Sqrt[p^2 + q^2 - 2 p q Cos[\[Phi]p - \[Phi]q]]/
Sqrt[30^2 + q^2 - 60 q Cos[\[Phi]q]], {p, 0, 50}, {q, 0,
50}, {\[Phi]p, 0, 2 \[Pi]}, {\[Phi]q, 0, 2 \[Pi]}]


answer is -3120.06 + 0.187845 I

NIntegrate[-p q/p/Sqrt[p^2 + q^2 - 2 p q Cos[\[Phi]p - \[Phi]q]]/
q, {p, 0, 50}, {q, 0, 50}, {\[Phi]p, 0, 2 \[Pi]}, {\[Phi]q, 0,
2 \[Pi]}]


answer is -4638.95 + 1.66533*10^-16 I

with imaginary parts

I saw some other post with similar question, some answer suggest using

    Method -> "MonteCarlo", PrecisionGoal -> 1


It does work for me, gives result 1220.21 with no warning I don't know how much can I trust this result Because the answer changes slightly together with warning when I increase Precisiongoal

NIntegrate::maxp: The integral failed to converge after 50100 integrand
evaluations. NIntegrate obtained -3233.79 and 35.28825355386109 for the
integral and error estimates.


I need trustful answer for other part of the program, any suggestions?

• The p q/p reduces to q -- is that a typo? Jul 25 '17 at 4:11
• @Michael E2 not typo, it can be reduced to q
– p.s
Jul 25 '17 at 4:27
• @p.s The integrand in NIntegrate is different from that in $F(k)$. The integrand should read $\frac{q-\sqrt{k^2-2 k p \cos (x)+p^2}}{\sqrt{k^2-2 k p \cos (x)+p^2} \sqrt{p^2-2 p q \cos (x-y)+q^2}}$ Jul 25 '17 at 8:27
• corrected the typo
– p.s
Jul 25 '17 at 17:23

We can do one integral symbolically and cut the integration region in half by symmetry. We still need to raise working precision to get only a few digits of accuracy:

intDp = Integrate[
-p q/p/Sqrt[p^2 + q^2 - 2 p q Cos[ϕp - ϕq]]/Sqrt[30^2 + q^2 - 60 q Cos[ϕq]],
{p, 0, 50},
Assumptions -> 0 < p < 50 && 0 < q < 50 && 0 < ϕp < 2 Pi && 0 < ϕq < 2 Pi];

2 NIntegrate[
intDp, {q, 0, 50}, {ϕp, 0, 2 π}, {ϕq, 0, ϕp},
WorkingPrecision -> 50, PrecisionGoal -> 6] // AbsoluteTiming


NIntegrate::slwcon: Numerical integration converging too slowly...

(*  {94.9907, -3313.8596254170187313848506470289485935167540168596}  *)


The NIntegrate::slwcon is a warning but not necessarily an error.

EDIT #3 complete revision

This text is intended to exhibit the current status of the problem comprising the discussions.

Main issues are

1. transformation of p and q to cartesian coordinates
2. convergence of the integral in the origin (obvious after 1.)
3. introduction of explicit cut-off for large p and q
4. model case to judge the accuracy of NIntegrate

For the ease of reading I start again from the beginning.

Let us write down the twofold two-dimensional integral step by step.

Define length of vector x

abs[x_] := Sqrt[x.x]


Vectors pv, qv and kv (x and y are angles between 0 and 2 pi)

pv = p {Cos[x], Sin[x]};
qv = q {Cos[y], Sin[y]};
kv = {k, 0};


The volume elements are

d^2p = p dp dx
d^2 q = q dq dy


Hence the integrand of F(k) is

int = Simplify[
p q 1/(abs[pv] abs[pv - qv] ) (1/abs[kv - qv] - 1/q), {p > 0, q > 0}]

(* Out[1020]= (q - Sqrt[k^2 + q^2 - 2 k q Cos[y]])/(Sqrt[p^2 + q^2 - 2 p q Cos[x - y]] Sqrt[k^2 + q^2 - 2 k q Cos[y]]) *)


In Latex

$$int = \frac{q-\sqrt{k^2-2 k q \cos (y)+q^2}}{\sqrt{k^2-2 k q \cos (y)+q^2} \sqrt{p^2-2 p q \cos (x-y)+q^2}}$$

Typical behaviour of the integrand

r := RandomReal[];

qr = -Log[r]; pr = -Log[r]; xr = 2 \[Pi] r; yr = 2 \[Pi] r;
Plot[int /. {p -> pr, q -> qr, x -> xr, y -> yr}, {k, 0, 15},
PlotRange -> All,
PlotLabel -> "Monte-Carlo check of integrand\nas a function of k",
AxesLabel -> {"k", "int"}]


Polar coordinates in 2-dimensional p-q-space

It is convenient to consider p and q formally as cartesian coordinates and make a coordinate transformation to polar coordinates (notice the factor c from the volume element dpdq = dc du)

intc = Simplify[
c int /. {p -> c Cos[u], q -> c Sin[u]}, {c > 0, k > 0, 0 < x < 2 \[Pi],
0 < y < 2 \[Pi], 0 < u < \[Pi]/2}]

(* Out[1086]= (c Sin[u] - Sqrt[
k^2 - 2 c k Cos[y] Sin[u] +
c^2 Sin[u]^2])/Sqrt[(k^2 - 2 c k Cos[y] Sin[u] + c^2 Sin[u]^2) (Cos[u]^2 +
Sin[u]^2 - Cos[x - y] Sin[2 u])] *)


in Latex

$$intc = \frac{c \sin (u)-\sqrt{c^2 \sin ^2(u)-2 c k \sin (u) \cos (y)+k^2}}{\sqrt{(1-\sin (2 u) \cos (x-y)) \left(c^2 \sin ^2(u)-2 c k \sin (u) \cos (y)+k^2\right)}}$$

Monte-Carlo check of integrand

cr = -Log[r]; xr = 2 \[Pi] r; yr = 2 \[Pi] r; ur = r \[Pi]/2;
Plot[intc /. {c -> cr, x -> xr, y -> yr, u -> ur}, {k, 0, 30},
PlotRange -> All,
PlotLabel ->
"Monte-Carlo check of integrand\ntransformed to Cartesian coordinates\nas a \
function of k", AxesLabel -> {"k", "intc"}]


Sumarizing the (numerical) integral of the OP is given by the function

f[kk_, R_] :=
Integrate[intc /. k -> kk, {c, 0, R}, {u, 0, \[Pi]/2}, {x, 0, 2 \[Pi]}, {y,
0, 2 \[Pi]}]
fN[kk_, R_] :=
NIntegrate[
intc /. k -> kk, {c, 0, R}, {u, 0, \[Pi]/2}, {x, 0, 2 \[Pi]}, {y, 0,
2 \[Pi]}]


Here we have introduced a cutoff radis R in order to avoid the divergence at c=[Infinity].

A typical numerical result is

fN[30, 10]


During evaluation of In[1158]:= 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. >>

During evaluation of In[1158]:= 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 -745.356 and 27.317349155628584 for the integral and error estimates. >>

(* Out[1158]= -745.356 *)


Reliability

In order to get a feeling of the reliability of the numeral results we study the integral dxdydu over the integrand at c = 0

intc0 = -(1/Sqrt[1 - Cos[x - y] Sin[2 u]]);


The exact value of the integral over intc0 is (see proof below)

ic0 = -((2 \[Pi]^4)/(Gamma[5/8]^2 Gamma[7/8]^2));
% // N

(* Out[1139]= -79.7335 *)


This value can be used to ckeck the prcision of various numerical approaches

ic0N1 = NIntegrate[intc0, {u, 0, \[Pi]/2}, {x, 0, 2 \[Pi]}, {y, 0, 2 \[Pi]}]


During evaluation of In[1125]:= 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. >>

During evaluation of In[1125]:= 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 -79.4062 and 0.14460520038353433 for the integral and error estimates. >>

(* Out[1125]= -79.4062 *)

ic0N2 = NIntegrate[intc0, {u, 0, \[Pi]/2}, {x, 0, 2 \[Pi]}, {y, 0, 2 \[Pi]},
Method -> "MonteCarlo", PrecisionGoal -> 2, WorkingPrecision -> 10]

(* Out[1126]= -78.43086734 *)


We observe very good agreement using the first method and confirm the error estimate of Mathematica.

Calculation of the exact value of the integral over intc0

Brute force Mathematica gives

Integrate[intc0, {u, 0, \[Pi]/2}, {x, 0, 2 \[Pi]}, {y, 0,
2 \[Pi]}] (* Mathematica error *)

(* Out[1120]= 0 *)


Which is oviously wrong!

We can find the correct result expanding the square root into a series and integrating term by term

Sum[Binomial[-(1/2), n] (-1)^n Cos[x - y]^n Sin[2 u]^n, {n,
0, \[Infinity]}] // FullSimplify

(* Out[1133]= 1/Sqrt[1 - Cos[x - y] Sin[2 u]] *)

Integrate[Binomial[-(1/2), n] (-1)^n Cos[x - y]^n Sin[2 u]^n, {u, 0, \[Pi]/2},
Assumptions -> n > -1]

(* Out[1134]= ((-1)^n Sqrt[\[Pi]]
Binomial[-(1/2), n] Cos[x - y]^n Gamma[(1 + n)/2])/(2 Gamma[1 + n/2])


The x and y integral is

Simplify[2 \[Pi] Integrate[ Cos[x]^n , {x, 0, 2 \[Pi]}],
n \[Element] Integers]

(* Out[1113]= -(((-2)^(1 + n) (1 + (-1)^n) \[Pi]^3)/(Gamma[1/2 - n/2]^2 Gamma[1 + n])) *)


And the sum becomes

Sum[-(((-2)^(1 + n) (1 + (-1)^n) \[Pi]^3)/(
Gamma[1/2 - n/2]^2 Gamma[1 + n])) ((-1)^n Sqrt[\[Pi]]
Binomial[-(1/2), n] Gamma[(1 + n)/2])/(2 Gamma[1 + n/2]), {n,
0, \[Infinity]}]

(* Out[1136]= (2 \[Pi]^4)/(Gamma[5/8]^2 Gamma[7/8]^2) *)

% // N

(* Out[1137]= 79.7335 *)


Q.E.D.

My original post

This is not a solution but a clarification of my statement of the comment. The next step to be done is to show that the integral is divergent so that a numerical treatment is ruled out.

There is an inconsistency in the OP: the integrand in the definition of F(k) differs from the integrand used in NIntegrate. Here we calculate the correct Expression of the Integrand.

Let us write down the twofold two-dimensional integral step by step.

Define the length of vector x

abs[x_] := Sqrt[x.x]


Vectors pv, qv and kv (x and y are angles between 0 and 2 pi)

pv = p {Cos[x], Sin[x]};
qv = q {Cos[y], Sin[y]};
kv = {k, 0};


The volume elements of the integrals are

d^2 p = p dp dx
d^2 q = q dq dy


Hence the integrand of F(k) is

int = Simplify[
p q 1/(abs[pv] abs[pv - qv] ) (1/abs[kv - qv] - 1/q), {p > 0, q > 0}]

(* Out[133]= (q - Sqrt[k^2 + q^2 - 2 k q Cos[y]])/(Sqrt[p^2 + q^2 - 2 p q Cos[x - y]] Sqrt[ k^2 + q^2 - 2 k q Cos[y]]) *)


In Latex

$$int = \frac{q-\sqrt{k^2-2 k q \cos (y)+q^2}}{\sqrt{k^2-2 k q \cos (y)+q^2} \sqrt{p^2-2 p q \cos (x-y)+q^2}}$$

• Thanks for this statement. I copied the wrong line from MMA and it is corrected now. The expression you gave is correct and simplified. Sorry I did not check carefully. The integral is convergent at $p,q\to0$ and there is an up cutoff at high p q .
– p.s
Jul 25 '17 at 17:24
• Can you prove your statement: "The integral is convergent at p,q→0" ? Jul 25 '17 at 21:33
• by counting power $\frac{0^4}{0^3}$
– p.s
Jul 25 '17 at 21:45