3
$\begingroup$

Consider the following code:

DistrToy[thh_, EM_] = 
  Exp[-((EM*(1 - 0.5 Cos[thh]^4))/125)]*Cos[thh]^20; 
TableDataTemp1 = 
   Table[{thh, Eh, DistrToy[thh, Eh]}, {thh, Subdivide[0, Pi, 314]}, {Eh, 125., 
     3000., 5.}];
TableData := Flatten[TableDataTemp1, {2, 1}]
DistrToyInterpolated = 
 Interpolation[TableData, InterpolationOrder -> 1]
Timing[NIntegrate[DistrToy[thh, Eh], {thh, 0, Pi}, {Eh, 125, 3000}]]
Timing[NIntegrate[
  DistrToyInterpolated[thh, Eh], {thh, 0, Pi}, {Eh, 125, 3000}]]

Here I define a function DistrToy[thh,Eh], make a table with rows {{thh,Eh, DistrToy[thh,Eh]}} and then make an interpolation. In the end I compare the time which is needed to integrate the initial function and the interpolation over the domain of the interpolation.

I found that the times of the integrations are different: the interpolated function is integrated 40 times longer than the initial function (0.078 s vs 3.125 s)! This makes all related computations extremely slow.

How to fix this issue?

$\endgroup$
  • $\begingroup$ I modified your code so it runs. I also changed the data table so that it includes Pi, otherwise evaluating your interpolating function at Pi generates errors. $\endgroup$ – Carl Woll Jun 4 at 19:49
  • $\begingroup$ @CarlWoll : thank you! I apologize for such inconveniences. $\endgroup$ – John Taylor Jun 4 at 19:50
3
$\begingroup$

Interpolating functions have lots of weak singularities, which require dense sampling to get an accurate estimate or require the integration region to be broken up according to the pieces of the piecewise interpolation, which is what the "InterpolationPointsSubdivision" strategy does.

i1 = NIntegrate[
    DistrToyInterpolated[thh, Eh], {thh, 0, Pi}, {Eh, 125, 3000}, 
    Method -> {"InterpolationPointsSubdivision", 
      "MaxSubregions" -> 315*3000/5, (* slight overestimate from interpolation *)
      Method -> {"GaussKronrodRule", "Points" -> 2}}, 
    MaxRecursion -> 0]; // AbsoluteTiming
i1
(*
  {0.286635, Null}
  75.15398362972101`
*)

Integrate[] will antidifferentiate an InterpolatingFunction exactly:

int = Integrate[DistrToyInterpolated[thh, Eh], thh, Eh]; // AbsoluteTiming
i2 = int /. {thh -> Pi, Eh -> 3000}
i1 - i2   (* error between NIntegrate and Integrate is zero *)
(*
  {1.29825, Null}
  75.15398362972101`
  0.                 
*)

Because the interpolating function is of low order, the Gauss-Kronrod rule integrates exactly without any recursive subdivision. The slowness, such as it is, is due to the number of subdivisions from the interpolation.

$\endgroup$
  • $\begingroup$ Thank you! I have few questions about this approach. 1) Consider the integration using AdaptiveMonteCarlo method. This integration, although less accurate, provides the answer for the integral over the interpolated function faster (although not enough) than the approach using InterpolationPointsSubdivision and GaussKronkodRule. However, by replacing InterpolationPointsSubdivision by AdaptiveMonteCarlo I do not change everything. Is there a way to use AdaptiveMonteCarlo method and speedup the integration? 2) Does this approach work for non-rectangular array? $\endgroup$ – John Taylor Jun 5 at 8:16
  • 1
    $\begingroup$ @JohnTaylor (1) I don't know much about the Monte Carlo methods. They're sort of a method of last resort in my mind, trading off accuracy for speed. If speed is a priority and only one or two digits of precision is acceptable, then it's probably the way to go. (2) Non-rectangular regions don't seem to be supported by InterpolationPointsSubdivision. That's a shame. $\endgroup$ – Michael E2 Jun 5 at 18:11
2
$\begingroup$

You can use NDSolveValue to integrate:

NDSolveValue[
    {
    Derivative[1,1][f][x,y] == DistrToyInterpolated[x,y],
    f[0,y] == DistrToyInterpolated[0,y],
    f[x,125] == DistrToyInterpolated[x,125]
    },
    f[Pi,3000],
    {x,0,Pi},
    {y,125,3000}
] //AbsoluteTiming

{0.072736, 75.0375}

$\endgroup$
  • $\begingroup$ Thank you! How do you think, should this approach using NDSolveValue work with more complicated functions? $\endgroup$ – John Taylor Jun 4 at 19:53
  • $\begingroup$ Also, should it work with complicated domains of the integration (where the domain for one variable depends on another variable)? $\endgroup$ – John Taylor Jun 4 at 19:55
  • 1
    $\begingroup$ @JohnTaylor Hard to say without examples. Non-square domains might require some non-trivial work. $\endgroup$ – Carl Woll Jun 4 at 19:59

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.