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"}}
]
]

WorkingPrecision
might help. $\endgroup$ – Jinxed Feb 20 '15 at 22:25int[]
andinf
? $\endgroup$ – Dr. belisarius Feb 20 '15 at 22:39