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?

=================

Try Interpolation[data, InterpolationOrder -> All] – for unstructured grids it may help.

– Vitaliy Kaurov

Jun 24 ’13 at 13:31

Could you provide an example data set, or code to generate one?

– Simon Woods

Jun 24 ’13 at 15:45

Thank you but that does not help.

– devilswild

Jun 25 ’13 at 9:28

1

This is a question of sampling, not interpolation. Plot3D is simply using a coarser mesh than ListPlot3D. To get similar results you could tell Plot3D to use a similar mesh with the options PlotPoints -> 100, MaxRecursion -> 0

– Simon Woods

Jun 25 ’13 at 10:05

Oh my god……. and i was thinking about that problem for weeks now. Thank you so much,

– devilswild

Jun 25 ’13 at 10:09

=================

2 Answers

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(πx/2)2f(x)=\sin(\pi x / 2)^2 and interpolate it in the interval {−1,1}\{-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/√(1−x2)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:

Create an Interpolation function

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

Near a single point, Taylor series is clearly superior, but expanding in terms of Chebyshev polynomials gives a better global bound. So, the phrase “best approximation” is dependent on your goal.

– rcollyer

Jun 25 ’13 at 12:46

@rcollyer yes you’re absolutely right. i guess the best approximation is from the plot itself and i was tempted to write something about getting a plot and extract the interpolated data from there….but then the issue changed and simon gave an answer the OP was actually searching for. But I have to say, by means of defending myself. The original question, as I understood, was about how ListPlot3D does the interpolation…now it became a more practical one.

– Stefan

Jun 25 ’13 at 19:19

but maybe I just understood the question wrong, which wouldn’t be the first … sigh

– Stefan

Jun 25 ’13 at 19:21

not likely. I was just nit-picking on nomenclature.

– rcollyer

Jun 25 ’13 at 19:49

Understood 🙂 jokester

– Stefan

Jun 25 ’13 at 20:29

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 xyxy 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 xyxy coordinates:

ListPlot3D[points, Mesh -> All, NormalsFunction -> None,

ViewPoint -> {0, 0, Infinity}]

DelaunayMesh[Most /@ points]