-1

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!

Michael E2
  • 235,386
  • 17
  • 334
  • 747
EZLearner
  • 295
  • 2
  • 6

3 Answers3

4

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}

Carl Woll
  • 130,679
  • 6
  • 243
  • 355
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}
Alex Trounev
  • 44,369
  • 3
  • 48
  • 106
  • 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) – EZLearner Oct 13 '18 at 12:30
  • 2
    There are no oscillations, I calculated in 1 second In[7]:= i2[10^-6, 1] Out[7]= {1.33688, -11.1869} – Alex Trounev Oct 13 '18 at 12:40
2

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

  • 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 –  Oct 13 '18 at 12:53
  • 1
    @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 03:08
  • 1
    Maybe the OP wanted a bad answer quickly, which is admittedly easy to do. ;) – J. M.'s missing motivation Oct 14 '18 at 04: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. –  Oct 14 '18 at 08:28