2
$\begingroup$

I am using Interpolation to construct an InterpolatingFunction from several points. I do not need a higher order InterpolationOrder than 1.

I wonder, in this case, is the contructed InterpolatingFunction simply a piecewise linear function? My goal is to Integrate that function, but due to performance reasons, I need to get the primitive function (evaluate the non-definite integral first). But I am wondering, whether this can be done and whether the result would be piecewise quadratic function so that the operation $F(x_i)$, where $F(x)=\int f(x)\,{\rm d} x$ and $f(x)$ is the InterpolatingFunction is as fast as plugging into a piecewise defined $\tilde{F}(x) = a_ix^2+b_ix +c_i$.

So in Mathematica, the code should look like:

f=Interpolation[data, InterpolationOrder -> 1];
F[x_]:=Integral[f[x],x];
result = Exp[F[#]]&/@biglistofx;

Are there any other performance issues that I should be aware when using Interpolation and its integration?

If the Interpolation does not produce a piecewise linear function, what would be my best approach to make the above calculation as fast as possible (note that Length@biglistofx $\sim 10^7$)?

$\endgroup$
4
$\begingroup$

Integrate will compute the antiderivative of an InterpolatingFunction:

data = Table[Sin[x], {x, 1., 10.}];
f = Interpolation[data, InterpolationOrder -> 1];
F[x_] = Integrate[f[x], x] (* with =, not := *)

Mathematica graphics

Plot[{f[x], F[x]}, {x, 1, 10}]

Mathematica graphics

Timing (InterpolatingFunction and Exp automatically threads over lists, so there's no need for Map; it's faster this way, too):

result = Exp[F[#]] &@ RandomReal[{1, 10}, 10^7]; // AbsoluteTiming
(*  {18.2411, Null}  *)

Length@result
(*  10000000  *)
$\endgroup$
  • 1
    $\begingroup$ You can also use Derivative[-1][f] to generate the pure function (i.e. InterpolatingFunction without the argument x) antiderivative. $\endgroup$ – Carl Woll Mar 14 '17 at 17:06
  • $\begingroup$ But are the InterpolatingFunction and its integral simple functions like piece wise defined linear/quadratic functions or are they more complex Mathematica internal expressions (and thus requiring more computation time)? Or in other words, can they be displayed as a piecewise definition? $\endgroup$ – leosenko Mar 14 '17 at 18:19
  • 1
    $\begingroup$ @leosenko The interpolating function is computed via an algorithm, no doubt optimized for its task. It's about 100 times slower than a simple math function like Exp and about 6 times faster than the correspond Piecewise function. I'd be surprised (and interested) if the interpolating function can be beat. $\endgroup$ – Michael E2 Mar 14 '17 at 18:51
1
$\begingroup$

Apart from Michael's and Carl's suggestions, you can also use NDSolve[] to generate the interpolant corresponding to the integral of the original InterpolatingFunction[]:

data = Table[Sin[x], {x, 1., 10.}];
f = Interpolation[data, InterpolationOrder -> 1];
F = NDSolveValue[{ff'[x] == f[x], ff[1] == 0}, ff, {x, 1, 10}];

Plot[{f[x], F[x]}, {x, 1, 10}]

plot of interpolant and its integral

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