How to integrate functions of linearly interpolated data?

Question

At first, consider integration of pure InterpolatingFunction.

Importing some data (works in v.9, for earlier versions one can use this link to download zipped M file which can be imported by data=Import["data.zip", "data.m"]):

data = First@
   Import["http://webbook.nist.gov/cgi/cbook.cgi?JCAMP=C67561&Index=17&Type=IR"];

Transforming the data according to an integrand:

integrand[ww_, k_] := k/ww
dataF = {#1, integrand[#1, #2]} & @@@ data; 

Now interpolation and integration:

In[13]:= int = Interpolation[dataF, InterpolationOrder -> 1];
ni = NIntegrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}]

During evaluation of In[13]:= NIntegrate::ncvb: NIntegrate failed to
converge to prescribed accuracy after 9 recursive bisections in ww
near {ww} = {1053.057940138439}. NIntegrate obtained 0.000035786622823277624
and 1.6106181332828143*^-6 for the integral and error estimates. >>

Out[14]= 0.0000357866

I am confused by unexpectedly high absolute error (see next) and strange error Message generated by NIntegrate when integrating function of linearly interpolated data. Note that NIntegrate fails to converge on linear (!) piece of the polyline:

In[15]:= Nearest[dataF[[All, 1]], 1053.057940138439`, 2]

Out[15]= {1053.0705299372542, 1053.010264968627}

Integrate is able to integrate pure InterpolatingFunctions:

i = Integrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}]

0.00003537589327711996

We can check it by direct integration:

exactSum = 
 MovingAverage[dataF[[All, 2]], 2].Differences[dataF[[All, 1]]]

0.00003537589420603186

The differences:

{ni, i} - exactSum

{4.10729*10^-7, -9.28912*10^-13}

One can see the large error for NIntegate output. Increasing PrecisionGoal and MaxRecursion does not give much:

ni = NIntegrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}, 
   MaxRecursion -> 25, PrecisionGoal -> 16] // Quiet
exactSum - ni

0.000035342690008474926

3.32042*10^-8

At this point my questions are:

1. Why NIntegrate fails to converge when interpolating the polyline? Why it has biggest trouble not even on a node but on a linear piece of the polyline?
2. Is it possible to achieve for NIntegrate at least the same precision as it is for Integrate? Or, better, as it is for direct summation?

Now I wish to integrate a function of linearly interpolated data:

int2 = Interpolation[data, InterpolationOrder -> 1];
ni2 = NIntegrate[int2[ww]/ww, {ww, data[[1, 1]], data[[-1, 1]]}, 
   MaxRecursion -> 22, PrecisionGoal -> 19] // Quiet

0.000035342690032342044

Integrate cannot work with functions of interpolated data at all. Direct summation gives:

area[{w1_, k1_}, {w2_, k2_}] = 
  Integrate[(a + b ww)/ww, {ww, w1, w2}, PrincipalValue -> True] /. 
    First@Solve[a + b w1 == k1 && a + b w2 == k2, {a, b}] // Simplify;
exactSum = Total[area @@@ Partition[data, 2, 1]]

exactSum - ni2

0.000035375894186915904

3.32042*10^-8

Again, I get very large error. It is not always possible to get closed-form formula for direct summation. So the final question is: what is the best way to integrate smooth functions of linearly interpolated data with high precision?

Answer 1

NIntegrate has many advanced options that let you control which algorithms and strategies it will use. I'm quite sure that you can find a set of options that will make NIntegrate work well enough for the desired task, but of course these numerical algorithms will never be quite as fast and precise as an exact solution, which your sums and Integrates results basically are.

A relatively obvious but poorly documented and probably still not optimal choice for integrands containing interpolating functions is the preprocessor method "InterpolationPointsSubdivision". With it and the strategy "LocalAdaptive", I already get a substantial improvement by just "blindly" playing with the available options, e.g. for your first example:

ni = 
 NIntegrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}, 
  Method -> {"InterpolationPointsSubdivision", 
    Method -> {"LocalAdaptive", "Partitioning" -> 10}},
  PrecisionGoal -> 9
  ]
exactSum - ni

(*
==> 0.0000353759
*)

(*
==> 6.65151*10^-13
*)

and also for your second example:

ni3 = 
 NIntegrate[int2[ww]/ww, {ww, data[[1, 1]], data[[-1, 1]]}, 
  Method -> {"InterpolationPointsSubdivision", 
    Method -> {"LocalAdaptive", "Partitioning" -> 15}},
  PrecisionGoal -> 9]
ni3 - exactSum

(*
==> 0.0000353759
*)

(*
==> 1.99279*10^-12
*)

With the given options NIntegrate doesn't throw any errors/warnings but is relatively slow. I'm not considering myself an expert in numeric integration and it looks like your own background would probably let you choose something even better. I would suggest to look at the tutorial advanced numeric integration and see how far you get. I see there some options which I don't fully understand but could well imagine to be of help for this case. You also might try the WRI support, if such a question is forwarded to the developers you might get a very good answer (if not, answers are not always very impressive). If my experience with NDSolve can be extrapolated to NIntegrate I would expect there is a set of options which should provide a good result in decent time, but that kind of problem (56000 data points with several sharp peaks and quite some noice or oscillations) is clearly something that needs an approach which gives NIntegrate some hints about what the characteristics of that function are...

Edit: as Alexey has found the improvement in precision is actually independent of the method setting "InterpolationPointsSubdivision", it's rather the combination of "LocalAdaptive", "Partitioning" and PrecisionGoal that makes a difference. Thus for this example these somewhat simpler settings seem to give results with equal quality:

 NIntegrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}, 
  Method -> {"LocalAdaptive", "Partitioning" -> 10},
  PrecisionGoal -> 9
 ]

I don't understand the precise effect of each of these settings in every detail and just stopped playing when I found a combination that seemed to work well. For other problems it might make sense to play with any of these or -- even better -- to study the problem and documentation in more depth...

Answer 2

You can try to force NIntegrate into doing a similar approach than you did by calculating the sum directly. This reduces your error quite a bit, although I don't understand the reasons behind this exercise:

ni = NIntegrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}, 
  Method -> {"TrapezoidalRule", 
    "Points" -> Length[dataF]}, MaxRecursion -> 0]
(* difference is -1.97858*10^-11 *)

Answer 3

Solution using NIntegrate of Mathematica 6+ (9.0.1)

One can achieve for integrating InterpolatingFunction by NIntegrate the same precision as with Integrate by specifying explicit set of Exclusions and MaxRecursion->0 (although the timings are disappointing):

integrand[ww_, k_] := k/ww
dataF = {#1, integrand[#1, #2]} & @@@ data;
int = Interpolation[dataF, InterpolationOrder -> 1];
ni1 = NIntegrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}, 
    MaxRecursion -> 0, PrecisionGoal -> 16, 
    Exclusions -> dataF[[2 ;; -2, 1]]]; // AbsoluteTiming
ni2 = NIntegrate[int[ww], Evaluate@Prepend[dataF[[All, 1]], ww], 
    MaxRecursion -> 0, PrecisionGoal -> 16]; // AbsoluteTiming
i1 = Integrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}]; // AbsoluteTiming
exactSum = 
   MovingAverage[dataF[[All, 2]], 2].Differences[
     dataF[[All, 1]]]; // AbsoluteTiming
{ni1, ni2} - i1
{ni1, ni2, i1} - exactSum
{ni1, ni2, i1, exactSum} // InputForm

During evaluation of In[1]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 0 recursive bisections in ww near {ww} = {1034.51}. NIntegrate obtained 0.00003537589327712014` and 5.3740022820734247`*^-20 for the integral and error estimates. >>

Out[5]= {171.515625, Null}

During evaluation of In[1]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 0 recursive bisections in ww near {ww} = {575.35}. NIntegrate obtained 0.000035375893282044965` and 3.621199488108469`*^-12 for the integral and error estimates. >>

Out[6]= {150.656250, Null}

Out[7]= {0.250000, Null}

Out[8]= {0.031250, Null}

Out[9]= {1.76183*10^-19, 4.925*10^-15}

Out[10]= {-9.28912*10^-13, -9.23987*10^-13, -9.28912*10^-13}

Out[11]//InputForm=
{0.00003537589327712014, 0.000035375893282044965, 0.00003537589327711996, 
 0.00003537589420603182}

Using the "LocalAdaptive" strategy and non-zero MaxRecursion in the same manner one can increase precision of integration of smooth function of InterpolatingFunction:

int2 = Interpolation[data, InterpolationOrder -> 1];
ni3 = NIntegrate[int2[ww]/ww, {ww, data[[1, 1]], data[[-1, 1]]}, 
    MaxRecursion -> 12, PrecisionGoal -> 16, Exclusions -> data[[2 ;; -2, 1]],
     Method -> "LocalAdaptive"]; // AbsoluteTiming
ni4 = NIntegrate[int2[ww]/ww, Evaluate@Prepend[dataF[[All, 1]], ww], 
    MaxRecursion -> 12, PrecisionGoal -> 16, 
    Method -> "LocalAdaptive"]; // AbsoluteTiming
area[{w1_, k1_}, {w2_, k2_}] = 
  Integrate[(a + b ww)/ww, {ww, w1, w2}, PrincipalValue -> True] /. 
    First@Solve[a + b w1 == k1 && a + b w2 == k2, {a, b}] // Simplify;
exactSum2 = Total[area @@@ Partition[data, 2, 1]]; // AbsoluteTiming
{ni1, ni2} - exactSum
{ni1, ni2, exactSum} // InputForm

Out[13]= {197.421875, Null}

Out[14]= {169.859375, Null}

Out[16]= {1.015625, Null}

Out[17]= {-9.28912*10^-13, -9.23987*10^-13}

Out[18]//InputForm=
{0.00003537589327712014, 0.000035375893282044965, 0.00003537589420603182}

---

NIntegrate of Mathematica 5.2

I should also point out that in Mathematica 5.2 one can get the same precision using only MaxRecursion and PrecisionGoal options and much faster, addition of exclusions makes it even faster:

In[4]:= << Statistics`DataSmoothing`

data = << data.m;
integrand[ww_, k_] := k/ww
dataF = {#1, integrand[#1, #2]} & @@@ data;
int = Interpolation[dataF, InterpolationOrder -> 1];
ni1 = NIntegrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}, 
    MaxRecursion -> 25, PrecisionGoal -> 16]; // AbsoluteTiming
ni2 = NIntegrate[int[ww], Evaluate@Prepend[dataF[[All, 1]], ww], 
    MaxRecursion -> 0, PrecisionGoal -> 16, 
    SingularityDepth -> 0]; // AbsoluteTiming
i = Integrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}]; // AbsoluteTiming
exactSum = 
   MovingAverage[dataF[[All, 2]], 
     2].(Rest@dataF[[All, 1]] - Most@dataF[[All, 1]]); // AbsoluteTiming
{ni1, ni2} - i
{ni1, ni2, i} - exactSum
{ni1, ni2, i, exactSum} // InputForm

From In[4]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration being 0, oscillatory integrand, or insufficient WorkingPrecision. If your integrand is oscillatory try using the option Method->Oscillatory in NIntegrate. More\[Ellipsis]

From In[4]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 26 recursive bisections in ww near ww = 1033.9707945203033`. More\[Ellipsis]

Out[9]= {79.5937500 Second, Null}

From In[4]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 2 recursive bisections in ww near ww = 575.2930004161132`. More\[Ellipsis]

Out[10]= {37.2656250 Second, Null}

Out[11]= {0.2656250 Second, Null}

Out[12]= {0.0468750 Second, Null}

Out[13]= {-5.34626*10^-14, 6.48336*10^-13}

Out[14]= {-9.82374*10^-13, -2.80576*10^-13, -9.28912*10^-13}

Out[15]//InputForm=
{0.000035375893223657405, 0.00003537589392545566, 0.00003537589327711996, 
 0.000035375894206031784}

In[16]:= int2 = Interpolation[data, InterpolationOrder -> 1];
ni3 = NIntegrate[int2[ww]/ww, {ww, data[[1, 1]], data[[-1, 1]]}, 
    MaxRecursion -> 25, PrecisionGoal -> 16]; // AbsoluteTiming
ni4 = NIntegrate[int2[ww]/ww, Evaluate@Prepend[data[[All, 1]], ww], 
    MaxRecursion -> 12, PrecisionGoal -> 16]; // AbsoluteTiming
area[{w1_, k1_}, {w2_, k2_}] = 
  Integrate[(a + b ww)/ww, {ww, w1, w2}, PrincipalValue -> True] /. 
    First@Solve[a + b w1 == k1 && a + b w2 == k2, {a, b}] // Simplify;
exactSum2 = Total[area @@@ Partition[data, 2, 1]]; // AbsoluteTiming

{ni3, ni4} - exactSum2
{ni3, ni4, exactSum2} // InputForm

From In[16]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration being 0, oscillatory integrand, or insufficient WorkingPrecision. If your integrand is oscillatory try using the option Method->Oscillatory in NIntegrate. More\[Ellipsis]

From In[16]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 26 recursive bisections in ww near ww = 1033.9707945203033`. More\[Ellipsis]

Out[17]= {84.7031250 Second, Null}

From In[16]:= NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration being 0, oscillatory integrand, or insufficient WorkingPrecision. If your integrand is oscillatory try using the option Method->Oscillatory in NIntegrate. More\[Ellipsis]

From In[16]:= NIntegrate::tmap: NIntegrate is unable to achieve the tolerances specified by the PrecisionGoal and AccuracyGoal options because the working precision is insufficient.  Try increasing the setting of the WorkingPrecision option.

Out[18]= {39.8593750 Second, Null}

Out[20]= {1.8281250 Second, Null}

Out[21]= {-9.83531*10^-13, -9.3084*10^-13}

Out[22]//InputForm=
{0.00003537589321637961, 0.000035375893269070094, 0.00003537589419991056}

---

All the above timings I have got on the same machine.