Numerical Integration of a multidimensional Integral of a singular function

I need to perform a numerical integration as part a computing matrix elements for a given block diagonal matrix. I've been trying several approaches with NIntegrate but none have been successful and are really slow. Here is the code:

l = 1;
f[u_, v_, k_] := E^((-I)*k*u) - E^((-I)*k*v);
g[u_, v_, k_] := E^((-I)*k*u) + E^((-I)*k*v) - 2*Cos[k*l];
W[u_, v_, x_, y_, ϵ_] := (-4*N[Pi])^(-1))*Log[(-0.0116681^2)*(u - x - I*ϵ)*(v - y - I*ϵ)];
F[u_, v_, x_, y_, ϵ_, m_, n_] = FullSimplify[f[u, v, m]*W[u, v, x, y, ϵ]*Conjugate[f[x, y, n]]];
G[u_, v_, x_, y_, ϵ_, m_, n_] = FullSimplify[g[u, v, m]*W[u, v, x, y, ϵ]*Conjugate[g[x, y, n]]];


Then I integrate:

mfi5 = Table[NIntegrate[F[u, v, x, y, 10^(-5), a, b]*Boole[(u - x)*(v - y) < 0], {u, -l, l}, {v,-l, l}, {x, -l, l}, {y, -l, l}, Method ->{"GlobalAdaptive", Method -> {"MultidimensionalRule"}}, MinRecursion -> 3, MaxPoints-> 100000, PrecisionGoal -> 3, AccuracyGoal -> 5, WorkingPrecision ->  8] + NIntegrate[F[u, v, x, y, 10^(-5), a, b]*Boole[(u - x)*(v - y) > 0], {u, -l, l}, {v, -l, l}, {x, -l, l}, {y, -l, l},Method -> {"GlobalAdaptive",Method -> {"MultidimensionalRule"}}, MinRecursion -> 3, MaxPoints-> 100000, PrecisionGoal -> 3, AccuracyGoal -> 5,WorkingPrecision -> 8], {a, k[[25]], k[[sk5[[2]]]], stk5},{b, k[[25]], k[[sk5[[2]]]], stk5}];
mgi5 = Table[NIntegrate[G[u, v, x, y, 10^(-5), a, b]*Boole[(u - x)*(v - y) < 0], {u, -l, l}, {v, -l, l}, {x, -l, l}, {y, -l, l}, Method -> {"GlobalAdaptive", Method -> {"MultidimensionalRule"}}, MinRecursion -> 3, MaxPoints -> 100000, PrecisionGoal -> 3, AccuracyGoal -> 5, WorkingPrecision -> 8] + NIntegrate[G[u, v, x, y, 10^(-5), a, b]*Boole[(u - x)*(v - y) > 0], {u, -l, l}, {v, -l, l}, {x, -l, l}, {y, -l, l}, Method -> {"GlobalAdaptive", Method -> {"MultidimensionalRule"}}, MinRecursion -> 3, MaxPoints -> 100000, PrecisionGoal -> 3, AccuracyGoal -> 5, WorkingPrecision -> 8], {a, κ[[25]], κ[[sκ5[[2]]]], stκ5},{b, κ[[25]],κ[[sκ5[[2]]]], stκ5}];


where the table indexes are given by:

κ1 = Flatten[Values[NSolve[Tan[x*l] == 2*x*l && 0 < x < 100, x, Reals]]];
κ2 = Table[((2*j - 1)/(2*l))*N[Pi], {j, 33, 800}];
κ = Join[κ1, κ2];
k = Table[i*(N[Pi]/l), {i, 1, 800}];
sk5 = Position[k, val_ /; val >= 500] // Flatten;
sk10 = Position[k, val_ /; val >= 1000] // Flatten;
sk15 = Position[k, val_ /; val >= 1500] // Flatten;
sk20 = Position[k, val_ /; val >= 2000] // Flatten;
sk25 = Position[k, val_ /; val >= 2500] // Flatten;
sκ5 = Position[κ, val_ /; val >= 500] // Flatten;
sκ10 = Position[κ, val_ /; val >= 1000] // Flatten;
sκ15 = Position[κ, val_ /; val >= 1500] // Flatten;
sκ20 = Position[κ, val_ /; val >= 2000] // Flatten;
sκ25 = Position[κ, val_ /; val >= 2500] // Flatten;
stk5 = Floor[k[[sk5[[1]]]]/4];
stk10 = Floor[k[[sk10[[1]]]]/4];
stk15 = Floor[k[[sk15[[1]]]]/4];
stk20 = Floor[k[[sk20[[1]]]]/4];
stk25 = Floor[k[[sk25[[1]]]]/4];
stκ5 = Floor[κ[[sκ5[[1]]]]/4];
stκ10 = Floor[κ[[sκ10[[1]]]]/4];
stκ15 = Floor[κ[[sκ15[[1]]]]/4];
stκ20 = Floor[κ[[sκ25[[1]]]]/4];
stκ25 = Floor[κ[[sκ25[[1]]]]/4];


However the numerical integration is not converging

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 GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 26.0048 +1.33834*10^-15\ I and 8.159908870291751 for the integral and error estimates."


And I know for a fact that the integration result is wrong since I know that the eigenvalues are all of absolute value less than three and come in pairs $\{\lambda_i, 1-\lambda_i\}$. The $W$ function is badly divergent but the $i\epsilon$ term added is used to regularize it, so I don't know where the convergence problem comes from.

I tryed modifying integration properties such as:"WorkingPrecison","MaxErrorIncreases"; and with "Montecarlo" and "AdaptiveMonteCarlo" as suggested in other posts. Nothing has fixed the problem yet.

Any suggestions on how to tackle the integral or which integration strategy to use and avoid the convergence problems will be of great help!

• In the definition of mfi5 the second term uses Integrate rather than NIntegrate. I suspect you meant Nintegrate. If so, kindly correct your question. – Jack LaVigne Jun 13 '17 at 0:02
• You go to a great deal of trouble defining various \[Kappa] values but they are not used in the tables mfi5 nor mfg5. Typo? – Jack LaVigne Jun 13 '17 at 0:06
• Yes typos, typos everywhere! I think I have corrected them and added some more options I've tweaked inside NIntegrate. – Matt Rest Jun 13 '17 at 4:29

This is a partial answer, for the integration of your F[u, v, x, y, 10^(-5), a, b]

The idea is, to do as much integration, as possible, analytically.

Doing the whole integration at once, I got kernel crash, because of to less memory. Therefore I expanded the F.

fse = FullSimplify[F[u, v, x, y, 10^(-5), a, b],
Assumptions -> {a, b, x, y, v, u} \[Element] Reals] // Expand


This gives 4 terms easy to integrate and 4 mor complicated terms (not shown here)

In order to help Mathematica integrate analytically I rationalized the numbers and divided the Boole[(u - x)*(v - y) < 0] into two integrals with integration limits from {u, -1, x}, {v, y, 1} and {u, x, 1}, {v, -1, y}.

Here you see this is the same

Manipulate[{RegionPlot[(u - x)*(v - y) < 0, {u, -1, 1}, {v, -1, 1},
MaxRecursion -> 5],
RegionPlot[-1 < u < x && y < v < 1 ||
x < u < 1 && -1 < v < y, {u, -1, 1}, {v, -1, 1},
MaxRecursion -> 5]}, {x, -1, 1}, {y, -1, 1}]


Now doing 16 times analytical integration in the two dimensions u and v, what takes quite a few minutes

Table[intuv[j][x_, y_, a_, b_] =
Assuming[{a, b, x, y, v, u} \[Element] Reals && -1 <= x <=
1 && -1 <= y <= 1 && 1/(1 + x) \[Element] Reals &&
1/a \[Element] Reals,
FullSimplify[Integrate[fse[[j]], {u, -1, x}, {v, y, 1}]]], {j, 1,
8}]

Table[intuv[j + 10][x_, y_, a_, b_] =
Assuming[{a, b, x, y, v, u} \[Element] Reals && -1 <= x <=
1 && -1 <= y <= 1 && 1/(1 + x) \[Element] Reals &&
1/a \[Element] Reals,
FullSimplify[Integrate[fse[[j + 10]], {u, -1, x}, {v, y, 1}]]], {j,
1, 8}]


The remaining x and y integration is done numericaly. Since some intuv were conditional expressions, I splitted it.

Total@Table[
NIntegrate[intuv[j][x, y, 1, 1], {x, -1, 1}, {y, -1, 1}], {j, 1,
4}] +
Total@Table[
NIntegrate[
intuv[j][x, y, 1, 1][[1]] Boole[
intuv[j][x, y, 1, 1][[2]]], {x, -1, 1}, {y, -1, 1}], {j, 5, 8}] +
Total@Table[
NIntegrate[intuv[j][x, y, 1, 1], {x, -1, 1}, {y, -1, 1}], {j, 11,
14}] +
Total@Table[
NIntegrate[
intuv[j][x, y, 1, 1][[1]] Boole[
intuv[j][x, y, 1, 1][[2]]], {x, -1, 1}, {y, -1, 1}], {j, 15,
18}]

(*   -1.54262 + 1.06587 I   *)
`
• Oh great thanks for your answer! Indeed, trying to compute the whole integral numerically is quite demanding for Mathematica. I will try to do as much of the integration analytically since I believe the integral can be handled analytically by Mathematica. However, I don't quite get why you partition the integration region as you did, Isn't it simpler just to take constant integration limits? And seeing the results you obtain wouldn't you recommend to do all the integration analytically rather than using analytical integration in the u and v variables and Numerical integration in a and y? – Matt Rest Jun 20 '17 at 18:10