# Problem with NIntegrate over a highly-oscillatory integrand

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!

• InterpolationPointsSubdivision is only intended for integrands that involve an InterpolatingFunction object, so no surprise that it failed here. – J. M. is away Oct 13 '18 at 12:11
• @J.M.issomewhatokay. Your'e correct but it's not an issue - specifying Method->"OscillatorySelection", for example, returns the same error. – EZSlaver Oct 13 '18 at 12:28
• The inner integral should be expressible in terms of the sine and cosine integrals, so that you are left with a 1D integral. Have you tried it? – J. M. is away Oct 13 '18 at 12:30
• No, but I'll try. – EZSlaver Oct 13 '18 at 12:31
• I only have gedanken Mathematica, so please try With[{a = 1*^-6, b = 1}, NIntegrate[(CosIntegral[b] - CosIntegral[x]) Sinc[x] - Cos[x] (SinIntegral[b] - SinIntegral[x])/x, {x, a, b}]] – J. M. is away Oct 13 '18 at 12:34

## 4 Answers

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:= % // AbsoluteTiming

Out= {0.007766, -0.7030662536921237781417087351943610040114}


so you can try to increase maxrecursion, etc etc

• It is even faster than the automatic method. – Alex Trounev Oct 13 '18 at 12:50
• It gives an error, but yes it is and this is why I posted it. Every method I tried gives a similar result, but this is the fastest I found and one can go crazy with accuracy etc etc – Konstantinos Oct 13 '18 at 12:53
• @EZSlaver The answer is fast, but the accuracy is very bad, only about 4 digits of precision are correct. Strange that this is the accepted answer when the OP requested 25 digits of precision. – Carl Woll Oct 14 '18 at 3:08
• Maybe the OP wanted a bad answer quickly, which is admittedly easy to do. ;) – J. M. is away Oct 14 '18 at 4:38
• I never said that my answer is the precise one, I just made a suggestion as some sort of a starting point and maybe one could build from there. And that was before I saw @J.M.iscomputer-less comment. – Konstantinos Oct 14 '18 at 8:28

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}

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}

• I need high accuracy because the integrand is very oscillatory, and as I've said, I need it down to very values of a (and of b as well) – EZSlaver Oct 13 '18 at 12:30
• There are no oscillations, I calculated in 1 second In:= i2[10^-6, 1] Out= {1.33688, -11.1869} – Alex Trounev Oct 13 '18 at 12:40

It works with the Gauss-Kronrod Rule, too:

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

Si2[1/10, 1] // AbsoluteTiming
(*  {0.370963, -0.7031515701355190323189234140040463855837}  *)

Si2[1/10^6, 1] // AbsoluteTiming
(*  {14.9981, -11.18694435428611911727667338594134335609}  *)


Side note: Si2[0.1, 1] generates a NIntegrate::slwcon warning and doesn't finish after several minutes. Ii seems machine precision creeps into the NumericalFunction used for the integrand. The integrand's argument is rescaled so that when it runs from 0 to 1, x runs from 0.1 to 1. The rescaling is computed with the machine precision 0.1. It seems this causes the error effectively to be bounded below by machine-precision rounding error. It can be observed that the error plateaus by incrementing MaxRecursion one step at a time or by using IntegrationMonitor.