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?

Improve this question
$\endgroup$
2
  • $\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
    Commented Jun 4, 2019 at 19:49
  • $\begingroup$ @CarlWoll : thank you! I apologize for such inconveniences. $\endgroup$
    – John Taylor
    Commented Jun 4, 2019 at 19:50

2 Answers 2

Reset to default
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.

Improve this answer
$\endgroup$
2
  • $\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
    Commented Jun 5, 2019 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
    Commented Jun 5, 2019 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}

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thank you! How do you think, should this approach using NDSolveValue work with more complicated functions? $\endgroup$
    – John Taylor
    Commented Jun 4, 2019 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
    Commented Jun 4, 2019 at 19:55
  • 1
    $\begingroup$ @JohnTaylor Hard to say without examples. Non-square domains might require some non-trivial work. $\endgroup$
    – Carl Woll
    Commented Jun 4, 2019 at 19:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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