How to understand the result of integrating this `InterpolatingFunction`

Question

I have this function

f[x_] = Integrate[Interpolation[{{1, 2}, {2, 4}}, x, InterpolationOrder -> 0], x];

that I plot along with its integral, which clearly has the wrong slope (the yellow curve should have 0 slope at x=1). Differentiating the integral (which seems to run correctly) doesn't yield the original function at all:

Plot[Evaluate[{f[x], Integrate[f[x], x], D[Integrate[f[x], x], x]}], {x, 0, 3}]

I don't think the discontinuity of the integral is much of a deal, because the integral is continuous in the range where InterpolatingFunction doesn't warn about extrapolation.

What did I miss?

Answer 1

If "I think it should be considered a bug" is answer, then this is an answer. Otherwise, it's just an explanation of what gets constructed by the OP's code.

The computed integral of f is equivalent to

zz = Interpolation[{{{1}, 0, 0, 1}, {{2}, 2, 4, 2}}, InterpolationOrder -> 0]

and has the same sort of erroneous interpolation, e.g., zz[1] returns -1 instead of 0.

yy = Head@Integrate[f[x], x];
yy === zz
(*  True  *)

The problem seems to be that the integral of f[x] is not promoted from an order 0 to an order 1 interpolation. I believe that's a bug. Normally Integrate raises the interpolation order by 1.

However, the bugs aren't done yet. I tried changing the InterpolationOrder of the OP's example to 1, and it didn't fix everything. The integral was okay, but the derivative of the integral, which is quadratic, was wrong, that is, it seems to be cubic (or higher)!

Plot[
 Evaluate@{f[x], Integrate[f[x], x], D[Integrate[f[x], x], x]},
 {x, 1, 2}]

Note that Integrate does a symbolic antiderivative of an InterpolatingFunction. Since an InterpolatingFunction is a piecewise polynomial, this shouldn't be surprising. Like with other results of Integrate, you perhaps should not rely on the constant of integration that might be added, but it seems always to start from zero at the initial end point of the interval. Integrate constructs a new InterpolatingFunction by specifying function data at each grid point as follows: If at $x_j$, we have given $$f(x_j), f'(x_j), \dots, f^{(n)}(x_j)\,,$$ then in the antiderivative $F$ constructed by Integrate there will be $$F(x_j), F'(x_j), \dots, F^{(n)}(x_j), F^{(n+1)}(x_j)\,.$$ It also (except in the OP's order 0 case) adds one to the interpolation order. (One might wonder that since it adds two data points, i.e., the values of the next derivative at each end point of a subinterval, the order could go up by two in Hermite interpolation. But since we're integrating polynomials, we know the degree only goes up by one.)

Derivatives of an InterpolatingFunction are computed in a very different way. They are computed by InterpolatingFunction from the interpolation data for the original function IF. All that happens in D[IF[x], x] or IF' is that an integer is incremented, an integer that indicates which derivative of the interpolating data to compute on a numeric input for x. It seems that this algorithm doesn't compute the derivative correctly in the above example. In fact, the derivatives appear to be the derivatives of a degree-5 polynomial, corresponding to the six data points seen in zz:

InterpolatingPolynomial[{{{1}, 0, 0, 1}, {{2}, 2, 4, 2}}, x]

See What's inside InterpolatingFunction[{{1., 4.}}, <>]? for information about the innards of interpolating functions.