Boundaries' bug in NIntegrate?

Question

This question is related to this answer of Antonov who discussed the customization of numerical integration strategies. Even though it's a inspiring method, there is a bug which make the result wrong.

One can found this bug without looking into the details of the codes. In that post, the correct half way integration is

Integrate[BesselJ[y, x^3], {x, 0, ∞}]

(*==> ConditionalExpression[Gamma[1/6 + y/2]/(3*2^(2/3)*Gamma[5/6 + y/2]), Re[y] > -(1/3)]*)

However, the printed result (HypergeometricPFQ involved) actually equals to

Integrate[BesselJ[y, x^3], {x, 0, 1}]

I've found that when the infinity boundary is involved, the returned value is {0,1} instead of {0,Infinity}.

Here I copy the codes (and modified a little bit) in that post:

Clear[SemiSymbolic];
Options[SemiSymbolic] = {"AnalyticalVariables" -> {}};
SemiSymbolicProperties = Options[SemiSymbolic][[All, 1]];
SemiSymbolic /: 
  NIntegrate`InitializeIntegrationStrategy[SemiSymbolic, nfs_, ranges_, 
   strOpts_, allOpts_] :=
  Module[{t, anVars},
   t = NIntegrate`GetMethodOptionValues[SemiSymbolic, SemiSymbolicProperties, 
     strOpts];
   If[t === $Failed, Return[$Failed]];
   {anVars} = t;
   SemiSymbolic[First /@ ranges, anVars]
   ];

SemiSymbolic[vars_, anVars_]["Algorithm"[regions_, opts___]] :=
  Module[{ranges, anRanges, funcs, t},

   (******************************************)
   (***The boundary of x is wrong, which is shown by the next Print code***)
   Print["ranges0=", Through[regions@"Boundaries"]];
   (*******************************************)

   ranges = Map[Flatten /@ Transpose[{vars, #@"Boundaries"}] &, regions];
   ranges = Map[Flatten, ranges, {-2}];
   anRanges = Map[Select[#, MemberQ[anVars, #[[1]]] &] &, ranges];
   ranges = Map[Select[#, ! MemberQ[anVars, #[[1]]] &] &, ranges];
   funcs = (#@"Integrand"[])@"FunctionExpression"[] & /@ regions;

   t = MapThread[
     Integrate[#1, Sequence @@ #2, 
       Assumptions -> (#[[2]] <= #[[1]] <= #[[3]] & /@ #3)] &, {funcs, 
      anRanges, ranges}];
   Print["SemiSymbolic::Integrate's result:", t];

   If[! FreeQ[t, Integrate], Return[$Failed]];

   Total[MapThread[
     NIntegrate[#1, Sequence @@ #2 // Evaluate, 
       Sequence @@ DeleteCases[opts, Method -> _] // Evaluate] &, {t, ranges}]]
   ];

One can see the bug by evaluating

NIntegrate[BesselJ[y, x^3], {x, 0, \[Infinity]}, {y, 0, 1}, Method -> {SemiSymbolic, "AnalyticalVariables" -> {x}}] 

The intermediate results printed are

ranges0={{{0,1},{0,1}}}
SemiSymbolic::Integrate's result:{(2^-y HypergeometricPFQ[{1/6+y/2},{7/6+y/2,1+y},-(1/4)])/((1+3 y) Gamma[1+y])}

(*==>0.371471*)

The integration is the same if one changes the upper limit of x from Infinity to 1.

In fact, Mathematica can directly give the numerical result not far away from the correct one even though it complains about the low convergence.

(*Direct numerical integration, not precise*)
NIntegrate[BesselJ[y, x^3], {x, 0, \[Infinity]}, {y, 0, 1}]
(*==>0.52435*)

(*This is the correct one*)
NIntegrate[Integrate[BesselJ[y, x^3], {x, 0, \[Infinity]}], {y, 0, 1}]
(*==>0.524448*)

After trying some codes a little bit, I've also found some strange results.

Reap[NIntegrate[Sin[x]/x, {x, 1, 10}, IntegrationMonitor :> (Sow[Through[#1@"Boundaries"]] &)]]
Reap[NIntegrate[Sin[x]/x, {x, 1, \[Infinity]}, IntegrationMonitor :> (Sow[Through[#1@"Boundaries"]] &)]]

(*==>{0.712265, {{{{{1, 10}}}}}}*)
(*==>{0.624713, {}}*)

The second output is empty. Why?

How to understand the extracted boundaries?

If it is a bug, how to fix it?

Answer 1

The second output is empty. Why?

Because of NIntegrate's automatic selection of algorithms for oscillatory integrals -- they (tend to) modify the integrand and the regions.

If it is a bug, how to fix it?

The bug was in the implementation of SemiSymbolic; see below.

As for NIntegrate's framework this, right now, is not a bug, it is a feature. One general way to deal with the transformations happening internally is to simply prevent them or do them before-hand using the options "SingularityHandler"->None and/or Method->{"UnitCubeRescaling","FunctionalRangesOnly"->False, _}. In this case that approach did not produce numerical results; see below.

The bug and the fix

This is a bug in my implementation of SemiSymbolic. I should have used the line:

ranges = Map[Flatten /@ Transpose[{vars, #@"OriginalBoundaries"}] &, 
   regions];

instead of the line:

ranges = Map[Flatten /@ Transpose[{vars, #@"Boundaries"}] &, regions];

I changed the code accordingly. With that new version we get:

In[39]:= NIntegrate[BesselJ[y, x^3], {x, 0, \[Infinity]}, {y, 0, 1}, 
 Method -> {SemiSymbolic, "AnalyticalVariables" -> {x}}]

During evaluation of In[39]:= ranges0={{{0,1},{0,1}}}

During evaluation of In[39]:= orig ranges0={{{0,\[Infinity]},{0,1}}}

During evaluation of In[39]:= SemiSymbolic::Integrate's result:{Gamma[1/6+y/2]/(3 2^(2/3) Gamma[5/6+y/2])}

Out[39]= 0.524448

A work-around that didn't work out

Attempting the work-around using "UnitCubeRescaling" mentioned above does not produce numerical results:

In[6]:= NIntegrate[BesselJ[y, x^3], {x, 0, \[Infinity]}, {y, 0, 1}, 
 Method -> {"UnitCubeRescaling", "FunctionalRangesOnly" -> False, 
   Method -> {SemiSymbolic, "AnalyticalVariables" -> {x}}}]

During evaluation of In[6]:= ranges0={{{0,1},{0,1}}}

During evaluation of In[6]:= orig ranges0={{{0,1},{0,1}}}

During evaluation of In[6]:= SemiSymbolic::Integrate's result:{Integrate[BesselJ[y,(-1+1/(1-x))^3]/(1-x)^2,{x,0,1},Assumptions->{0<=y<=1}]}

Out[6]= NIntegrate[BesselJ[y, x^3], {x, 0, \[Infinity]}, {y, 0, 1}, 
 Method -> {"UnitCubeRescaling", "FunctionalRangesOnly" -> False, 
   Method -> {SemiSymbolic, "AnalyticalVariables" -> {x}}}]

For successful application of this work-around see the explanations of AdaptiveNumericalLebesgueIntegration.m .