How does ListPlot3D interpolate data?

Question

I have data of the form {{x1, y1, f(x1, y1}, {x2, y2, f(x2, y2)}, ...}. I need to interpolate the data to get an interpolated f(x, y).

When I just do Interpolation[data], I get a function, but the result is not very smooth, independent on the interpolation order (there is a difference between the orders, but it does not get better).

When I plot the data with ListPlot3D[data], I get a really nice result. It is really smooth. When I change the interpolation order in ListPlot3D, the result does seem to not change, except when going from 0th order to 1st order.

So I would like to get a function, which looks like the result of ListPlot3D. How does ListPlot3D interpolate the data? Or can I get a function of the ListPlot3D result?

Edit

First of all, thanks a lot for the detailed answers. I just tried the InterpolationPolynomial, but it turned out that it is too demanding for our machine. Maybe, if I post two pictures showing what I'm seeing, one of you will show me how to interpolate my data in the proper way.

First the plot made with ListPlot3D, the nice one:

And the interpolated function plotted:

My data grid has the resolution of 0.2, so I have 101 x 101 data points; that is, the grid is based on {x, -10, 10, .2}, {y, -10, 10, .2}).

As one can clearly see, the second plot shows some edges and is not as smooth as the first one. Maybe someone has an idea?

Answer 1

I will answer the question posed in the title:

How does ListPlot3D interpolate data?

It uses a Delaunay triangulation of the projection of the points on the $xy$ plane. This is essentially the same as using first-order Interpolation (which at present is the only interpolation scheme available for unstructured data). Observe:

SeedRandom[0];
points = RandomReal[1, {100, 3}];

ListPlot3D[points, NormalsFunction -> None]

f = Interpolation[points, InterpolationOrder -> 1];
Plot3D[f[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> {0, 1}, 
 PlotPoints -> 25, MaxRecursion -> 5]

Compare the top view of the mesh used by ListPlot3D to the Delaunay triangulation of just the $xy$ coordinates:

ListPlot3D[points, Mesh -> All, NormalsFunction -> None, 
 ViewPoint -> {0, 0, Infinity}]

DelaunayMesh[Most /@ points]

Answer 2

Quality of Interpolation

I don't have any insight into the interpolation scheme of Mathematica, but I can describe how I would interpolate toward the best approximation.

Mathematica has a whole battery of numerical interpolation functionality. But how to achieve the same quality, as the specific plot do?

One interpolation function that comes to my mind here is InterpolatingPolynomial.

Let's take a function $f(x)=\sin(\pi x / 2)^2$ and interpolate it in the interval $\{-1, 1\}$.

n = 50;
f[x_] := Sin[Pi x/2]^2;
data = Table[{x, f[x]}, {x, -1, 1, 2/n}];

Now let's do some interpolation with data:

ipo[x_] = InterpolatingPolynomial[data // N, x];

And compare the original function with the interpolated polynomial:

Show[GraphicsGrid[{{
    Plot[#, {x, -1, 1}, Frame -> True, Axes -> False,
        DisplayFunction -> Identity]& /@ {f[x], ipo[x]}
}}]]

In order to qualify the interpolation there is a quantification step, to determine the error as a function of n. With this it is possible to determine the best approximation in some interval.

The above function is a smooth function on the specific interval, but what if we have to interpolate a function with a specific singularity near some interval?

If this is the case we have to deal with the so called Runge phenomenon.

Let's take for instance the function as defined in the above article on Wikipedia:

f[x_] := 1/(1 + 25 x^2);

If we re-evaluate everything we get this plot:

Here the interpolation oscillates toward the end of the interval. If something like this happens it is best to use nonequidistant x-values. On that specific case the best asymptotic density is given by the formula:

$1/\sqrt{(1 - x^2)}$

These nodes are known as Chebyshev nodes

So what we see here is, that the specific plot functions work roughly in two steps:

1. Create an Interpolation function
2. Qualify the Interpolation function through quantification

I have to leave now, but try to come back to your interesting question again and get more into the specifics; how I think Mathematica does the interpolation or how I'd do it if I were Mathematica (that's not the case, which is good luck for Mathematica).