2
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ 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])/(1 + x) + (Cos[1] Sin[x])/(1 + x), the integral of each of which it estimates at one go. The IntegrationMonitor is not called at all. $\endgroup$ – Michael E2 Nov 29 '16 at 18:46
  • $\begingroup$ @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. $\endgroup$ – luyuwuli Nov 30 '16 at 2:15
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Nov 30 '16 at 2:49
  • $\begingroup$ @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. $\endgroup$ – luyuwuli Nov 30 '16 at 3:07
  • $\begingroup$ Anton answers every week or so. I suppose he's busy. He may yet reply. $\endgroup$ – Michael E2 Nov 30 '16 at 3:21
3
$\begingroup$

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 .

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.