Numerical Integration of a multidimensional Integral of a singular function

Question

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!

Answer 1

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   *)