# How to integrate functions of linearly interpolated data?

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?

• Your code doesn't work with the data you have posted. In fact First@Import.... is a header and the table is not formatted properly. It seems that the first 38 elements are headers and so is the last. Can you use Import[...,"Table"] and select only numbers Select[#, NumberQ]/@ then see if you can reproduce it?
– gpap
Commented Feb 1, 2013 at 13:33
• Plot[Derivative[-1][int][ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}, Evaluated -> True] looks okay? Commented Feb 1, 2013 at 13:57
• @ Oleksandr Yes, it is OK. Commented Feb 1, 2013 at 14:21
• @gpap I just checked again, it works (version 9.0.0)! ListLogPlot[data, PlotRange -> All, Joined -> True] gives me this. Commented Feb 1, 2013 at 14:30
• @gpap In versions <9 my Import code will not work since support of the jdx (JCAMP-DX) format is added in version 9. I'll try to find a workaround. Commented Feb 1, 2013 at 14:41

## Update:

After identifying and fixing a problem with data, we show that the setting the "MaxSubregions" option of "InterpolationPointsSubdivision" to the number intervals created by the interpolation points produces a highly accurate answer in a reasonably short amount of time.

### Verify the data!

Ahhh, this is maddening: there's problem with the data. One of the points (the 7th) is out of order and causes exactSum to be off by -9.28912*10^-13, which coincides the "error" of my first answer. If we fix this by sorting the data and recomputing the "exact" sum, the error is greatly diminished. The first answer (see below) turns out to be highly accurate and fast (for NIntegrate).

Here is the offending datum:

data[[4 ;; 8]]

(*  {{575.171, -2.41302*10^-7}, {575.231, -6.9385*10^-8},
{575.291, 8.54164*10^-7}, {575.29, 2.06936*10^-6}, (* <-- here *)
{575.35, 3.80004*10^-6}} *)


Here is the corrected exact sum, exactSumSorted, and the comparison with exactSum:

integrand[ww_, k_] := k/ww
dataF = {#1, integrand[#1, #2]} & @@@ data;
dataFSorted = {#1, integrand[#1, #2]} & @@@ Sort@data;
exactSum = MovingAverage[dataF[[All, 2]], 2].Differences[dataF[[All, 1]]]
exactSumSorted = MovingAverage[dataFSorted[[All, 2]], 2] .
Differences[dataFSorted[[All, 1]]]
(*              !
0.00003537589420603182
0.000035375893277120126
^            *)

exactSumSorted - exactSum
(*
-9.28912*10^-13
*)


Comparison of first answer with the sums:

ni - exactSum
ni - exactSumSorted
(*
-9.28912*10^-13
-1.76183*10^-19
*)


Now the relative error,

(ni - exactSumSorted)/exactSumSorted
(*
-4.98031*10^-15
*)


is much more in line with what is expected.

### The second integral

I omitted to deal with the second integral in my first answer, but it was what got me curious about what was going on. The graph of int2[ww] / ww should be piecewise concave up and the trapezoid rule should overestimate the integral, but the reverse was happening.

We can get the exact sum of the integrals over the subregions of the interpolating function in a way similar to how Alexey Popkov got the exact sum for the integral of the linear interpolation. First observe the following. If value of a function y = f[x] is linearly interpolated between the values y1, y2, as x ranges from x1 to x2, then we can compute the integral of f[x]/x as follows:

Integrate[((1 - t) y1 + t y2)/((1 - t) x1 + t x2) (x2 - x1), {t, 0, 1},
Assumptions -> {x2 > x1 > 0}]
(*
-y1 + y2 + ((-x2 y1 + x1 y2) Log[x2/x1])/(x1 - x2)
*)


[Update 2: Further explanation. First, it was easier for me to think abstractly in terms of x and y; the correspondence is x = ww and y = int2[ww]. Switching to TeX, a linear interpolation between any two values $A$, $B$ is given by $$(1-t)\,A + t\, B, \quad 0 \le t \le 1\,.$$ (Clearly it is linear in $t$ and ranges from $A$ to $B$ as $t$ ranges from $0$ to $1$. Applying this to a subregion $x_1 \le x \le x_2$ where $y$ is interpolated linearly between the values $y_1 = {\tt int2}(x_1)$ and $y_2 = {\tt int2}(x_2)$, we can paramterize the interpolation by $$y = (1-t)\,y_1 + t\,y_2,\quad x = (1-t)\,x_1 + t\,x_2 = x_1 + t\,(x_2 - x_1), \quad 0 \le t \le 1.$$ Substituting into the integral, this becomes $$\int_{ww_1}^{ww_2} {{\tt int2}(ww) \over ww}\;d(ww) = \int_{x_1}^{x_2} {y \over x}\;dx = \int_{0}^{1} {(1-t)\,y_1 + t\,y_2 \over (1-t)\,x_1 + t\,x_2}\;(x_2-x_1)\,dt\,,$$ where the factor $x_2-x_1 = dx/dt$. This corresponds to the Mathematica integral above. To get the integral over the whole data set, we just total up the values of this integral over all of the subregions, which we do below to get exactSum2.]

Thus the exact value (ignoring machine precision error) can be calculated with

With[{data = Sort@data},
exactSum2 = Total[
-First @ Differences[               (* -x2 y1 + x1 y2 *)
Reverse @ Transpose @ Most[data] Transpose @ Rest[data]] *
Log @ Ratios[data[[All, 1]]] /    (* Log[x2/x1] *)
Differences[data[[All, 1]]]       (* x1 - x2 -- note minus in front *)
] + data[[-1, 2]] - data[[1, 2]]    (* telescoping y2 - y1 *)
]
(*
0.00003537589326824584
*)


The integral may be computed, using the same "InterpolationPointsSubdivision" options as in the first answer, with

int2 = Interpolation[Sort@data, InterpolationOrder -> 1];

(ni2 = NIntegrate[int2[ww]/ww,
Evaluate[{ww}~Join~First@int2["Domain"]],
Method -> {"InterpolationPointsSubdivision",
"MaxSubregions" -> 1 + Length[First@int2["Coordinates"]]}
]) // AbsoluteTiming
(*
{29.845430, 0.00003537589326907005}
^         *)


Compared to Alexey Popkov's similar version, we see the above is a little faster and agrees with Alexey's up to the last digit.

(ni2AP = NIntegrate[int2[ww]/ww,
Evaluate[{ww}~Join~First@int2["Domain"]],
Method -> {"InterpolationPointsSubdivision",
"MaxSubregions" -> 10^9,
Method -> {"LocalAdaptive", MaxPoints -> Infinity, Method -> {"TrapezoidalRule"}}}
]) // AbsoluteTiming
(*
{31.126306, 0.000035375893269070053}
^     ^    *)


The error is about the best that could be expected:

{ni2, ni2AP} - exactSum2
(*
{8.24204*10^-16, 8.2421*10^-16}
*)


For one thing, as a very rough estimate, we can estimate that about four digits or more of precision might be lost just from differences of the first coordinate:

Max @ Abs[(Rest[#] - Most[#])/Most[#]] &@ Sort@data[[All, 1]]
(*
0.000119089
*)


[Update 2: I should have stressed that this is a crude estimate I used to give me an idea of what a reasonable error between NIntegrate and exactSum might be. The errors are probably randomly distributed and some probably cancel others. Also the greatest loss of precision in a single difference would be indicated by Min instead of Max (about 6 digits). So the actual error might vary widely from this estimate. The calculation above is based on this sort of example: Consider the two numbers below and their absolute difference:

0.00001234123
0.00001234567
-------------
0.00000000444 (abs. val. of difference)


We went from 7 digits to 3 digits of precision. The relative difference is

(0.00001234567 - 0.00001234123)/0.00001234567
(*
0.00035964
*)


The number of digits lost is indicated by the logarithm:

Log10[0.00035964]
(*
-3.44413
*)


End of update 2.]

To get a complete "InterpolationPointsSubdivision", a sufficient number of subregions needs to be allowed with the "MaxSubregions" option. This allows NIntegrate to handle the 50K+ subintervals in the interpolating function int.

data = First @
Import["http://webbook.nist.gov/cgi/cbook.cgi?JCAMP=C67561&Index=17&Type=IR"];
integrand[ww_, k_] := k/ww
dataF = {#1, integrand[#1, #2]} & @@@ data;
int = Interpolation[dataF, InterpolationOrder -> 1];

(ni =
NIntegrate[
int[ww],
Evaluate[{ww} ~Join~ First@int["Domain"]],
Method -> {"InterpolationPointsSubdivision",
"MaxSubregions" -> 1 + Length[First@int["Coordinates"]]}]
) // AbsoluteTiming

(* {0.129937, 0.0000353759} *)


I'm not sure why one more subregion is needed than coordinates/data points. I suppose it is because N coordinates divide the real numbers into N+1 regions, ignoring the interval of integration.

## Postscript: More of the story

Every now and then I tried NIntegrate with some other option, even increasing WorkingPrecision, to see if the integral would get closer to the exact value. It was the stubborn resistance of NIntegrate to get the relative error down to close to 10^-15, or even 10^-12, that made me suspect that NIntegrate was right and exactSum was wrong.

Here is something I found along the way that pleased me. Nothing earth shattering, but when a numerical method returns an error that is essentially zero, it's fun.

I used arbitrary precision numbers with 30 digits of precision. They can be quite a bit slower. Here is the exact value.

(With[{data = SetPrecision[Sort@data, 30]},
exactSum2p = Total[
-First @ Differences[               (* -x2 y1 + x1 y2 *)
Reverse @ Transpose @ Most[data] Transpose @ Rest[data]] *
Log @ Ratios[data[[All, 1]]] /    (* Log[x2/x1] *)
Differences[data[[All, 1]]]       (* x1 - x2 -- note minus in front *)
] + data[[-1, 2]] - data[[1, 2]]    (* telescoping y2 - y1 *)
]) // AbsoluteTiming
(*
{0.722450, 0.00003537589326907000957268}
*)


The machine precision calculation took 0.008799 sec. NIntegrate will be much slower, too.

int2p = Interpolation[SetPrecision[Sort@data, 30], InterpolationOrder -> 1];
(ni2p = NIntegrate[int2p[ww]/ww,
Evaluate[{ww} ~Join~ First@int2p["Domain"]],
Method -> {"InterpolationPointsSubdivision",
"MaxSubregions" -> 1 + Length[First@int2p["Coordinates"]]},
WorkingPrecision -> 30
]) // AbsoluteTiming
(*
{160.618752, 0.0000353758932690700095726756654469}
*)


Here is a comparison of the results:

Print[ni2p // InputForm]
Print[exactSum2p // InputForm]
(*
0.00003537589326907000957267566544694703035073570136128996521.891735884879235
0.000035375893269070009572675665446947030350749449535372155830.
^ 22nd            ^ 39th (1st difference)  *)

ni2p - exactSum2p // InputForm
(*
026.343028466547207  (* StandardForm: 0.*10^-27 *)
*)


The last output means the first 26+ digits after the decimal place are zero. This corresponds to exactSum2p having almost 22 digits of precision and its first nonzero digit being the 5th after the decimal point.

• (+1). Your solution gives amazing speedup in integration of pure InterpolatingFunction and significant speedup in integration of the int2[ww]/ww. But "MaxSubregions" suboption is completely undocumented: searching on all Wolfram websites gives no results. Where did you find it? Commented Dec 25, 2013 at 16:49
• @AlexeyPopkov I put in a wrong option, e.g. Method -> {"InterpolationPointsSubdivision", "Foo" -> 1}, and Mma told the valid ones. Commented Dec 25, 2013 at 16:59
• It is amazing! Such a crucially needed option can be found only esoterically! Even the form Method -> {"InterpolationPointsSubdivision"} is undocumented (I cannot find it in the Documentation). Commented Dec 25, 2013 at 17:03
• It is interesting that NIntegrate[int2[ww]/ww, {ww, data[[1, 1]], data[[-1, 1]]}, Method -> {"InterpolationPointsSubdivision", "MaxSubregions" -> 10^9}] gives even better performance but "MaxSubregions" -> 10^10 is forbidden. Commented Dec 25, 2013 at 17:24
• @AlexeyPopkov After ?**Subdivision*, one finds NIntegrateInterpolationPointsSubdivision. Thence Options[NIntegrateInterpolationPointsSubdivision] yields the default settings, too. This is good, too: ?NIntegrateStrategiesDump* -- lots of options and properties. Commented Dec 25, 2013 at 17:56

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

• (+1) It is strange but it seems that specifying the "InterpolationPointsSubdivision" preprocessor changes nothing in both cases: compare NIntegrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}, Method -> {"LocalAdaptive", "Partitioning" -> 10}, PrecisionGoal -> 9] with NIntegrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}, Method -> {"InterpolationPointsSubdivision", Method -> {"LocalAdaptive", "Partitioning" -> 10}}, PrecisionGoal -> 9]. The outputs are identical. Commented Feb 2, 2013 at 2:22
• I also tried to switch "InterpolationPointsSubdivision" off and got the same result: NIntegrate[int[ww], {ww, dataF[[1, 1]], dataF[[-1, 1]]}, Method -> {"SymbolicPreprocessing", "InterpolationPointsSubdivision" -> False, Method -> {"LocalAdaptive", "Partitioning" -> 10}}, PrecisionGoal -> 9]. AbsoluteTimings are also identical. It is just "LocalAdaptive" strategy makes good job. Commented Feb 2, 2013 at 2:23
• @AlexeyPopkov: I tried to make sure that I just played around with those options and it turns out that you are right: With just "InterpolationPointsSubdivision" the difference to the "exact sum" is not better than Automatic, it is the combination of "LocalAdaptive", "Partitioning" and PrecisionGoal that seems to make a difference. Certainly there is room for more investigation and optimization... Commented Feb 2, 2013 at 19:58

UPDATE 2 start

As Michael has shown, semi-documented preprocessor "InterpolationPointsSubdivision" with undocumented "MaxSubregions" -> 10^9 suboption gives significant speedup in both cases. So the complete Method option should be:

Method -> {"InterpolationPointsSubdivision", "MaxSubregions" -> 10^9,
Method -> {"LocalAdaptive", MaxPoints -> Infinity,
Method -> {"TrapezoidalRule"}}}


The rest of the answer does not include this update.

It is interesting that switching off "SymbolicProcessing" suboption gives a little better performance in this case:

Method -> {"InterpolationPointsSubdivision", "MaxSubregions" -> 10^9,
"SymbolicProcessing" -> False}


UPDATE 2 end

UPDATE start

In both cases (integration of int[ww] and int2[ww]/ww, i.e. direct integration of the InterpolatingFunction and integration of the function of the InterpolatingFunction) selection of the "TrapezoidalRule" integration rule significantly improves the performance without noticeable loss of precision. So complete Method option should look like this in both cases:

Method -> {"LocalAdaptive", Method -> {"TrapezoidalRule"}}


The rest of the answer does not include this update.

UPDATE end

## 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]],
ni4 = NIntegrate[int2[ww]/ww, Evaluate@Prepend[dataF[[All, 1]], ww],
MaxRecursion -> 12, PrecisionGoal -> 16,
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]:= << StatisticsDataSmoothing

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.

• Commented May 25, 2016 at 7:29

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 *)

• My point is not integrating the InterpolatingFunction (it can be done with Integrate) but integrating a smooth function of linearly interpolated data. I even do not insist on using NIntegate for this. I still have no good solution for the case when the closed form of the integral cannot be derived. Commented Feb 3, 2013 at 11:32
• @halirutan: do you know where the "Points" option for "TrapezoidalRule" is documented and what exactly it does? Commented Feb 3, 2013 at 16:39
• @AlbertRetey Here it is documented. Commented Dec 27, 2013 at 18:50
• @AlexeyPopkov: thanks, that's useful. Actually it is where one could expect to find it, I obviously just didn't try hard enough... Commented Dec 27, 2013 at 23:36
• @AlbertRetey Sorry for not answering to your comment. Somehow I must have missed it. I see that it is already from February.. better late than never. Commented Dec 27, 2013 at 23:43