NIntegrate 2D highly oscillatory function

Question

I am trying to integrate a function, and the error I get is greater than the result.

So I need to calculate

NIntegrate[osc[n,L,u,x]*f1[u, x]/f2[u, x], {u, 0, 1}, {x, 1, inf}]

where

f1[u_,x_] := u(1-u) ( 1-Cos[x]-x Sin[x] )

f2[u_,x_] := (1-u)^2 x^2 + 2u(1-u) x Sin[x] + 2( 1-Cos[x] ) u^2

and

osc[n_,L_,u_,x_]:= Exp[-I u x(n+L)+2I n ArcTan[(u(1-Cos[x]))/(u Sin[x]+x(1-u))]]

Here, n and L are parameters: n is an integer (1, 2, 3, …), and L is real > 0

The best I got was by using Method -> "LevinRule":

For example, by using

nn=1; LL = 1;

I get

int[n_,L_,u_,x_] := osc[n,L,u,x]*f1[u, x]/f2[u, x];
inf = Infinity;
NIntegrate[ int[nn, LL, u, x], {u, 0, 1}, {x, 1, inf}, Method -> 
   {"LevinRule", "Kernel" -> osc[nn, LL, u, x]}]

(* -0.0741726 + 0.0851585 I *)

with an error of 0.18027814

Note that the error is greater than the result.

Answer 1

It seems to me that infinitely many, semi periodic singularities makes for a very difficult integral. If you could find a substitution that would tame all the singularities at once, you could probably succeed in truly taming this integral.

I agree with the OP that each singularity is integrable, but one has also to wonder whether the sum of the infinitely many singularities exists. My feeling is that it does. Even though the magnitude of the integrand looks like it falls off as $1/x$, the oscillatory nature of the function does not rule out convergence, and the numerical calculations I get stay around 0.08(6) + 0.08(1) I, not too far from the OP's result.

My thought was that since the singularities occur as $u \rightarrow 1^-$, $x \rightarrow 2 \pi n^-$, $n$ an integer, we might divide the infinite strip {u, 0, 1}, {x, 0, Infinity} into infinitely many rectangles {u, 0, 1}, {x, 2 Pi n, 2 Pi + 2 Pi n} and sum the integrals. (One has to subtract the integral over {u, 0, 1}, {x, 0, 1} to get the OP's integral, of course.) A singularity occurs only in one corner of each rectangle. This suggests that using Duffy's coordinates might be useful. The substitution transforms a singularity at a corner to one along an edge. If the singularity is of the right type, it can be quite effective. In this case, it helps, but the singularity persists to be a nuisance. However, the IMT singularity handler helps, at least to the extent that with a higher working precision, the integral can be computed with the desired precision.

Summing the integrals over the rectangles can be done with NSum. This means approximating the sum using approximate terms. I'm not real sure what the compounding of the error implies. Given that all terms were computed with a relative precision of MachinePrecision/2 digits and the real and imaginary parts of the terms all have the same sign (except the imaginary part of the first term), the sum should be nearly as good as NSum usually calculates, I suppose.

Another strategy I used was to isolate the singularity in the corner. Away from the corner NIntegrate can accurately compute the integral with the "LevinRule" method. Luckily the singularity can be predicted to lie in the rectangle {u, 1 - 1/Max[n, 1]/4, 1}, {x, 3 Pi/2 + 2 Pi n, 2 Pi + 2 Pi n}. We just have to figure out how to integrate this singularity.

For some reason Method -> {"DuffyCoordinates", "Corners" -> {{1, 1}}} was insufficient to compute the integrals with the desired precision. But the NIntegrate tutorial linked above includes code for manually performing the substitution. We can then call NIntegrate specifying the singularity handler "SingularityHandler" -> "IMT", which finally seems to conquer the problem.

So finally, here is the code that implements the above strategy. The function DuffyCoordinates is given in this section of the NItegrate tutorial. I would include it, but I'm not sure about copyright restrictions. Warning: It takes a long time, about a minute or so per rectangle; hence the results are memoized. Print statements are included for the impatient. The singular integral always gives a singularity warning NIntegrate::slwcon, so I supressed it with Quiet.

ClearAll[intpart7, f1, f2, osc, int];

(* OP's definitions *)
f1[u_, x_] := u (1 - u) (1 - Cos[x] - x Sin[x])
f2[u_, x_] := (1 - u)^2 x^2 + 2 u (1 - u) x Sin[x] + 2 (1 - Cos[x]) u^2
osc[n_, L_, u_, x_] := 
 Exp[-I u x (n + L) + 2 I n ArcTan[(u (1 - Cos[x]))/(u Sin[x] + x (1 - u))]]
nn = 1; LL = 1;
int[u_, x_] := osc[nn, LL, u, x]*f1[u, x]/f2[u, x]

(* integrates over the n-th rectangle *)
intpart7[n_Integer] := Block[{x, u},
   Print[n];
   intpart7[n] = (* memoize for reuse *)
    NIntegrate[            (* nonsingular integral *)
      Boole[0 <= u <= 1 - 0.25/Max[n, 1] || 2 Pi n <= x <= 2 Pi n + 3 Pi/2]*int[u, x],
      {u, 0, 1}, {x, 2 Pi n, 2 Pi + 2 Pi n},
      MaxRecursion -> 25, 
      Method -> {"GlobalAdaptive", 
        Method -> {"LevinRule", "Kernel" -> osc[nn, LL, u, x]}}] +
     Quiet[Print["."];
      NIntegrate @@ Join[  (* singular part *)
        DuffyCoordinates[
         int[u, x], {u, 1 - 1/Max[n, 1]/4, 1},
         {x, 3 Pi/2 + 2 Pi n, 2 Pi + 2 Pi n}, {1, 1}],
        {PrecisionGoal -> 8, WorkingPrecision -> 20, 
         MaxRecursion -> 25, 
         Method -> {"GlobalAdaptive", "SingularityHandler" -> "IMT"}}
        ], NIntegrate::slwcon]
   ];

Grand result

intzero = NIntegrate[int[u, x], {u, 0, 1}, {x, 0, 1}]
(sum = NSum[intpart7[n], {n, 0, Infinity}]) // AbsoluteTiming
sum - intzero
(*
  -0.0751088 + 0.0187646 I                (* integral 0 < x < 1 *)
  {1891.108293, -0.161177 + 0.0992955 I}  (* {time, sum} *)

  -0.0860682 + 0.0805309 I                (* final integral *)
*)

No error estimate warnings!

Remark: The computed terms may be inspected with ? intpart7 or DownValues[intpart7]. The indices n computed are

Cases[DownValues[intpart7], HoldPattern[intpart7[n_Integer]] :> n, Infinity]
(*
  {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
   19, 20, 21, 22, 23, 24, 25, 26, 30, 93305, 93306, 93307}
*)

Accuracy of results

We can see that using more terms in the estimation by NSum tends to increase the magnitude of the result, but by less and less. I let this run, while I was doing other stuff. It probably takes a total of two hours, including the time above. NSum does not seem to estimate the tail of the sum very accurately. It's possible that it diverges very slowly, but the OP seems to think it converges. I have not thought very much about it.

seq = Table[NSum[intpart7[n], {n, 0, Infinity}, NSumTerms -> n] - intzero, {n, 15, 100, 10}]
(*
{-0.0860682 + 0.0805309 I,
 -0.0861501 + 0.0805336 I,
 -0.0860646 + 0.0805342 I,
 -0.0861001 + 0.0805403 I,
 -0.0861384 + 0.080535 I,
 -0.0861612 + 0.0805338 I,
 -0.0861782 + 0.080533 I,
 -0.0861947 + 0.080534 I,
 -0.0862018 + 0.0805343 I}
*)

ListLinePlot[Abs@Differences@seq, DataRange -> {15, 95}]

Visualization of the Duffy transformation

The plots below show the real and imaginary parts of the integrand in the singular corner of the domain for n = 10 before and after the substitution of Duffy's coordinates. One can see that the singularity is moved to the edge as required by the "IMT" singularity handler.

Block[{n = 10, style, fn},
 {style[Re], style[Im]} = 
  Take["DefaultPlotStyle" /. (Method /. Charting`ResolvePlotTheme["Default", Plot3D]), 2];
 fn["int"] = {##} &;
 fn["Duffy"] = {First@DuffyCoordinates[##], {u, 0, 1/n}, {x, 1 - 2/n, 1}} &;
 GraphicsGrid@
  Table[Plot3D[part[#1], #2, #3, PlotRange -> 5, MaxRecursion -> 3, 
      PlotPoints -> 35, PlotStyle -> style[part], 
      AxesLabel -> Automatic, PlotLabel -> part[f]] & @@ 
    fn[f] @@ {int[u, x], {u, 1 - 1/Max[n, 1]/4, 1}, {x, 3 Pi/2 + 2 Pi n, 2 Pi + 2 Pi n}, {1, 1}},
   {part, {Re, Im}}, {f, {"int", "Duffy"}}
   ]
 ]

Answer 2

If you are happy with the findings so far, and just need help on dynamically creating the exclusions, you might want to proceed with the following index-generation code:

DeleteDuplicates@{#1,#2,Sequence@@Range[Ceiling[#2,2\[Pi]],Floor[#3,2\[Pi]],2\[Pi]],#3}&

which can be called with the variable and its lower and upper bounds like so, e.g.:

DeleteDuplicates@{#1,#2,Sequence@@Range[Ceiling[#2,2\[Pi]],Floor[#3,2\[Pi]],2\[Pi]],#3}&[x,1,10]

So, your NIntegrate becomes:

NIntegrate[int[nn, LL, u, x], {u, 0, 1},
  DeleteDuplicates@{#1,#2,
  Sequence@@Range[Ceiling[#2,2\[Pi]],Floor[#3,2\[Pi]],2\[Pi]],#3}&[x,xmin,xmax],
  Method->{"LevinRule", "Kernel"->osc[nn,LL,u,x]}]