Problem with NIntegrate over a highly-oscillatory integrand

Question

I'm trying to numerically evaluate the integral $$\int_{a}^{b}\mathop{\mathrm{d}x}\int_{x}^{b}\frac{\sin(x-y)}{xy}\mathop{\mathrm{d}y}$$ using Mathematica. To do that, I the function

Si2[a_, b_] := NIntegrate[Sin[x - y]/(x y), {x, a, b}, {y, x, b},
  AccuracyGoal -> 25, PrecisionGoal -> 25, WorkingPrecision -> 40,
  MaxRecursion -> 1000000, Method -> "InterpolationPointsSubdivision"];

However, running Si2[.1,1] gives the error

NIntegrate::inumr: The integrand Sin[x-y]/(x y) has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1},{0,1}}.

However, I'm not integrating over $x=y=0$ (which is an obvious singularity).

Two notes:

1. This is an example of usage. In practice, I need to evaluate this function for parameters much closer to $a=0$ (e.g. $\log_{10}(a)\sim-6$).
2. Note that I use the InterpolationPointsSubdivision method because I saw in various answers that it is a good method to evaluate numerically a highly-oscillatory integrand. I tried to use few other methods, but got the same error.

Any advice? Thanks!

Answer 1

Something like that?

Si2[a_, b_] := 
 NIntegrate[Sin[x - y]/(x y), {x, a, b}, {y, x, b}, 
   AccuracyGoal -> 25, PrecisionGoal -> 25, WorkingPrecision -> 40, 
   MaxRecursion -> 100000000, {Method -> "QuasiMonteCarlo"}] // Quiet
Si2[0.1, 1]

The result is

-0.7030662536921237781417087351943610040114

and it's not very slow

In[65]:= % // AbsoluteTiming

Out[65]= {0.007766, -0.7030662536921237781417087351943610040114}

so you can try to increase maxrecursion, etc etc

Answer 2

Here is @JM's answer. First integrate the interior integral:

i1[a_, b_] = Integrate[Sin[x-y]/y, {y, x, b}, Assumptions -> 0<x<b]

(CosIntegral[b] - CosIntegral[x]) Sin[x] + Cos[x] (-SinIntegral[b] + SinIntegral[x])

Then, use this integral:

int[a_, b_, opts:OptionsPattern[NIntegrate]] := NIntegrate[i1[a, b]/x, {x, a, b}, opts]

Reproducing previous results:

int[.1, 1, WorkingPrecision->40] //AbsoluteTiming
int[10^-6, 1, WorkingPrecision->40] //AbsoluteTiming

{0.413422, -0.7031515701355189970289164059287652664106}

{0.882268, -11.18694435428611911727667338594134335609}

Answer 3

Why such a high accuracy? You can simply calculate the integral without options

i2[a_, b_] := 
     NIntegrate[NIntegrate[Sin[x - y]/(x y), {y, x, b}], {x, a, b}] // 
       Quiet // AbsoluteTiming


 i2[1/10, 1]

Out[]= {0.697864, -0.703152}