# 4D NIntegrate with singularities

I need to integrate a function in a 4D region (x1,y1,x2,y2), which explodes whenever x1=x2&&y1=y2. The code is as following:

    {Lx, Ly} = {1/50, 1/50};

f1[x_, y_, n_, m_] := Cos[n*Pi*(Lx + x)/(2*Lx)]*Sin[m*Pi*(Ly + y)/(2*Ly)];
f2[x_, y_, n_, m_] := Sin[n*Pi*(Lx + x)/(2*Lx)]*Cos[m*Pi*(Ly + y)/(2*Ly)]*m*Lx/(n*Ly);
dist[x1_, y1_, x2_, y2_] := Sqrt[(x1 - x2)^2 + (y1 - y2)^2];

int[f_, n_, m_, i_, j_] := int[f, n, m, i, j] =
NIntegrate[f[x1, y1, n, m]*f[x2, y2, i, j]/dist[x1, y1, x2, y2],
{y1, -Lx, Lx}, {x1, -Lx, Lx}, {y2, -Ly, Ly}, {x2, -Lx, Lx}]


Then, I want to evaluate int using either f1 or f2, for positive integer (n,m,i,j), for instance using AbsoluteTiming[int[f1,1,2,1,2]]. I may need to compute these integrals thousands of times, for the different possible combinations of (n,m,i,j) - in the most complex cases, each of them may reach up to 14, for instance. So I must be able to calculate the integral as fast as possible.

However, the method is simply not working so well, for some. When I compute, for instance: int[f1,4,1,1,1], the computation simply lasts for ever and I receive some scary error messages:

NIntegrate::errprec: Catastrophic loss of precision in the global error estimate due to insufficient WorkingPrecision or divergent integral.


I have tried multiple integration strategies and methods, but the thing that seems to work best is to simply define the integration boundaries using {y1,-Ly,Ly},{x1,-Lx,Lx},{y2,-Ly,y1,Ly},{x2,-Lx,x1,Lx}, taking advantage of the integration order to define the singularity locations.

The integration method "DuffyCoordinates" performed quite fast, but I didn't know whether I was correctly using the "Corners" option, e.g.:

 int2[f_,n_,m_,i_,j_] := NIntegrate[f[x1, y1, n, m]*f[x2, y2, i, j]/dist[x1, y1, x2, y2],
{y1, -Lx, Lx}, {x1, -Lx, Lx}, {y2, -Ly, Ly}, {x2, -Lx, Lx},Method -> {"DuffyCoordinates","Corners" -> {0, 0, 1, 1}}]


... But it also gets stuck for (n,m,i,j)=(4,1,1,1). So I don't know what else to do:

• Am I doing something wrong with this couple of integration strategies, or missing something obvious?
• I suspect much of the problems are due to the integral actually evaluating to 0, maybe because the functions are somehow odd (whatever that means in 4D?), and I'm integrating them in a symmetric interval. But I don't know how I could define my function int in a way that it would check whether the integrand is odd in this 4D symmetrical integration volume.

Edit 15/08/2018 corrected the expression of the function dist[]

Cross-posted: http://community.wolfram.com/groups/-/m/t/1404998

• You can probably add some small value inside Sqrt to avoid singularities. Have you tried Monte-Carlo based approaches? – David Baghdasaryan Aug 15 '18 at 8:03
• But wouldn't that just become a different function? I tried Method->"MonteCarlo", and it was super fast, but setting different random seeds the results were very inconsistent. I might be doing something wrong, though, since I don't have experience working with MonteCarlo methods, but that really makes me doubt the method – miguel Aug 15 '18 at 9:24
• probably not much different if that value enough small. You can specify PrecisionGoal to a precision you need. See the documentation for that option. – David Baghdasaryan Aug 15 '18 at 9:54
• I've corrected the definition for the function dist[], that sqrt should be in the denominator of the integrand, I'm sorry for wasting your time. Thanks for the tips. These approaches seem to have some convergence issues, and the estimate for the error seems to be almost as large as the estimate for the value of the integral itself. So, for instance running int[f1,1,1,1,4], the integral and error estimates are 6.4062*10^-8 and 8.40064*10^-8. Also, running int[f1,1,1,1,1], the value obtained by the automatic method is 0.0000194422, whereas "AdaptiveQuasiMonteCarlo" gets 0.0000167858 – miguel Aug 15 '18 at 12:57
• Please include links to crossposts of the same question. (General policy is discussed here: mathematica.meta.stackexchange.com/questions/367/…) – Michael E2 Aug 15 '18 at 19:03

I am not sure how helpful this answer is. I have not analyzed the integrand and integral much, just applied couple of tricks.

## Plan

1. Identify which corners have singularities.

2. Split the integral to have one singularity only in the integral terms.

3. Specify the singularity corners to DuffyCoordinates.

## Definitions

Clear[int4]
Options[int4] = Options[NIntegrate];
int4[f_, n_, m_, i_, j_, opts : OptionsPattern[]] :=
NIntegrate[
f[x1, y1, n, m]*f[x2, y2, i, j]/dist[x1, y1, x2, y2], {y1, 0,
Lx}, {x1, 0, Lx}, {y2, 0, Ly}, {x2, 0, Lx}, opts]

Clear[int5]
Options[int5] = Options[NIntegrate];
int5[f_, n_, m_, i_, j_, opts : OptionsPattern[]] :=
NIntegrate[
f[x1, y1, n, m]*f[x2, y2, i, j]/dist[x1, y1, x2, y2], {y1, -Lx,
0}, {x1, -Lx, 0}, {y2, -Ly, Ly}, {x2, -Lx, 0}, opts]


## Identify singularity corners

The identification is done by using the trace of the sampling points of "AdaptiveMonteCarlo". We see that the corners {0,0,0,0} and {0,1,0,1} along the axes 2 and 4 have (most likely) singularities.

res =
Reap[int3[f1, 4, 1, 1, 1, MinRecursion -> 3, MaxRecursion -> 20,
Method -> "AdaptiveMonteCarlo", PrecisionGoal -> 2,
EvaluationMonitor :> Sow[{y1, x1, y2, x2}]]];

inds = Flatten[Table[{i, j}, {i, 1, 4}, {j, i + 1, 4}], 1];

ListPlot[RandomSample[res[[2, 1]][[All, #]], 12000],
AspectRatio -> Automatic, PlotLabel -> #] & /@ inds


## Split integral and specify corners

AbsoluteTiming[
res1 = int4[f1, 4, 1, 1, 1, MaxRecursion -> 20,
Method -> {"DuffyCoordinates", "Corners" -> {0, 1, 0, 1},
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 20000}},
PrecisionGoal -> 4]
]

(* {63.621, 5.30172*10^-7} *)

AbsoluteTiming[
res2 = int5[f1, 4, 1, 1, 1, MaxRecursion -> 20,
Method -> {"DuffyCoordinates", "Corners" -> {0, 0, 0, 0},
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> 20000}},
PrecisionGoal -> 4]
]

(* {35.0088, -6.6503*10^-7} *)

res1 + res2

(* -1.34857*10^-7 *)
`

These give error message but the results might be good enough.