# How does ListInterpolation work?

ListInterpolation takes an array of data and interpolates between the entries. Let us interpolate between the elements of the identity matrix. For reference, this is what the identity matrix looks like:

MatrixPlot[IdentityMatrix[5]]


Now let's interpolate, and request first order interpolation:

if = ListInterpolation[IdentityMatrix[5], {{0, 1}, {0, 1}}, InterpolationOrder -> 1];
Plot[if[x, x], {x, 0, 1}]


The result is a piecewise function as expected, but not a piecewise linear function! Why? What does first order interpolation mean here?

The kind of linear interpolation I am familiar with is what ListDensityPlot does:

ListDensityPlot[IdentityMatrix[5], DataRange -> {{0, 1}, {0, 1}}, Mesh -> All]


Construct a Delaunay triangulation of the data and interpolate on each triangle.

This is not what ListInterpolation does (or even Interpolation when the data is strictly on a square grid). What does ListInterpolation do, then?

Inspired by this question:

• Well, it is doing bilinear interpolation in your first example, I would think. You definitely will get something with quadratic pieces if you sample across diagonals. Jun 15, 2016 at 14:00
• @J.M. Well, we all have to learn at some point :-) Yes, you are right, and this is due to my ignorance of math. But I won't delete it. If you convert that to a one-line answer, I'll accept it. Jun 15, 2016 at 14:01
• Docs say the only methods available to ListInterpolation are Spline and Hermite Jun 15, 2016 at 14:02
• Currently my comment is not that pedagogically useful; I'll try to come up with a sufficiently illuminating example first, unless somebody beats me to it. :) Jun 15, 2016 at 14:04
• Szabolcs, please consider leaving the question up for a little longer. I was not aware of that behavior either, and I would find it instructive to see an interesting example, if @J.M. or others have the time / inclination. Jun 15, 2016 at 14:43

I'll just treat the simplest case of interpolating across four noncoplanar points. The extension to general grids is straightforward (one just joins multiple pieces appropriately).

Consider the following bivariate interpolating polynomial:

ip[x_, y_] = InterpolatingPolynomial[{{{0, 0}, 1}, {{0, 1}, -1/2},
{{2, 0}, 1/2}, {{2, 1}, 0}}, {x, y}];


and the following InterpolatingFunction[]:

if = ListInterpolation[{{1, -1/2}, {1/2, 0}}, {{0, 2}, {0, 1}}, InterpolationOrder -> 1];


They are the same:

Plot3D[ip[x, y] - if[x, y] // Chop, {x, 0, 2}, {y, 0, 1}]


Let's have a look at what ip[] looks like:

Expand[ip[x, y]]
1 - x/4 - (3 y)/2 + (x y)/2


Huh. That last term certainly isn't linear. What gives?

In a bilinear interpolation such as this one, although the procedure is something along the lines of "linearly interpolate for each grid line in one direction, and then linearly interpolate those results", the result is inevitably quadratic. This should come as no surprise to people familiar with the hyperbolic paraboloid, which is one of the simplest examples of a ruled surface, or a surface formed by sweeping out a line in space. In fact, z == 1 - x/4 - (3 y)/2 + (x y)/2 is indeed a hyperbolic paraboloid.

{Plot3D[ip[x, y], {x, 0, 2}, {y, 0, 1}, MeshFunctions -> {#1 - #2 &}],
Plot[ip[x, x], {x, 0, 1}]} // GraphicsRow