# Boundaries' bug in NIntegrate?

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 /:
NIntegrateInitializeIntegrationStrategy[SemiSymbolic, nfs_, ranges_,
strOpts_, allOpts_] :=
Module[{t, anVars},
t = NIntegrateGetMethodOptionValues[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, #[]] &] &, ranges];
ranges = Map[Select[#, ! MemberQ[anVars, #[]] &] &, ranges];
funcs = (#@"Integrand"[])@"FunctionExpression"[] & /@ regions;

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

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

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?

• The double exponential oscillatory rule doesn't subdivide the interval for NIntegrate[Sin[x]/x, {x, 1, Infinity}]. It translates the interval to {0, Infinity} and splits the integrand into (Cos[x] Sin)/(1 + x) + (Cos Sin[x])/(1 + x), the integral of each of which it estimates at one go. The IntegrationMonitor is not called at all. – Michael E2 Nov 29 '16 at 18:46
• @MichaelE2 In which cases will IntegrationMonitor not be called? You said that the integration "estimates at one go", I'm puzzled about this. I think NIntegrate[x,{x,0,1}] also estimates at one go, however, it still return the correct boundaries when IntegrationMonitor is called. – luyuwuli Nov 30 '16 at 2:15
• I don't know the full internals, but you can read up a little on the strategy in this tutorial. The DEO strategy approximates the integral with a truncated series, and perhaps does not use the "GlobalAdaptive" strategy of recursive subdivision. -- As for which strategies "Boundaries" works, I'm not sure.  Method -> "Trapezoidal" seems to always give {-Infinity, Infinity}, even for finite intervals. It might work reliably only for the global adaptive strategy. – Michael E2 Nov 30 '16 at 2:49
• @MichaelE2 Thanks, I'll spend some time studying that. However, I really wish the docs could be more clear about these ambiguities especially about the undocumented options. I've also comment below Antonov's answer, but he doesn't reply; that's why I make this post. – luyuwuli Nov 30 '16 at 3:07
• Anton answers every week or so. I suppose he's busy. He may yet reply. – Michael E2 Nov 30 '16 at 3:21

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];


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


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

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

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

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

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

Out= 0.524448


# A work-around that didn't work out

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

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

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

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

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

Out= 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 .