Interpolation and Smoothing

All Products Maple MapleSim

Home : Support : Online Help : Applications and Example Worksheets : Mathematical Visualization : Interpolation and Smoothing

Interpolated Plotting and Smoothing:

Example Worksheet

The topic of this worksheet is interpolation and smoothing of given 2-dimensional and 3-dimensional data. The organization is by dimension, task, and regularity of the data.

One common task is to generate only a plot of the data, as a curve or surface which passes through or approximates the data points.

Another task is to generate a procedure or piecewise spline expression which approximates the data and which can be queried for a value at any individual point nearby the original data points.

Another choice is whether to produce a curve or surface which passes through the given data points (and so interpolates them) or which simply approximates them (by smoothing). Is the data known to contain noise or measurement error? If the dependent data is expected to contain error or noise then it is reasonable to fit a smoothed surface or curve to the data, where the smoothed fit does not necessarily pass through or match the dependent data at the given data points. On the other hand, if the data is expected to be fully correct then it is reasonable to produce an interpolated curve or surface which passes through all the given data points.

A curve or surface may also be plotted directly from the data, or indirectly by first creating a procedure or (in the 2D case) a piecewise spline, and then subsequently supplying that to the plot or plot3d command. The constructed procedure would accept individual 1D or 2D independent data points and compute a dependent scalar value for that input.

Another relevant distinction for 3D data is whether the independent portion of the data lies on a regular grid in the 2D plane or whether it is comprised of points which are irregularly spaced in both directions.

This worksheet elaborates and compares these approaches for 2D and 3D plotting.

>	restart; with(LinearAlgebra): with(CurveFitting): with(Interpolation): with(plots): with(Statistics): randomize():

1-dimensional case

The x-data points are taken uniformly in the range, as the integers from 1 to 15. But they could just as well be unequally spaced, as long as they are ordered.

Here, the data2D is a 15x2 Matrix with the first column being the ordered x-data points and the second column being their y-values.

>	xdata := <($1..15)>: ydata := RandomVector(15, generator = 1.0 .. 5.0): data2D:=< xdata \| ydata >:

Here is the plot command producing a linear interpolation. The goal is to produce a smooth curve instead.

>	display( plot(data2D, color=black), pointplot(data2D, symbolsize=10, color=red) );

We produce an interpolating function using splines, and plot this function.

>	f := Interpolate(xdata, ydata): display(plot(f, 1..15, color=black), pointplot(data2D, symbolsize=10, color=red) );

The above method of evaluating the interpolating spline is simplest. It uses the Interpolate command from the Interpolation package. Behind the scenes, it uses the ArrayInterpolation command from the CurveFitting package, which is called in an efficient manner for all the interpolating points that the plotting procedure wants to evaluate. Using options, one can select between different types of interpolants.

>	g := Interpolate(xdata, ydata, method=cubic): h := Interpolate(xdata, ydata, method=radialbasisfunction, thinplatespline): display(plot([f, g, h], 1..15, color=[black, blue, green]), pointplot(data2D, symbolsize=10, color=red) );

Another use of an interpolant generated by the Interpolate command is to do numeric integration. The following commands determine the area under the three curves in the previous plot.

>	int(f, 1 .. 15, numeric); int(g, 1 .. 15, numeric); int(h, 1 .. 15, numeric);

$40.5747583655129$

$40.6203688066937$

$40.62177489$

(1.1)

Another option is local optimization. We have chosen the data randomly so far, which is not very suitable for demonstrating this functionality. Let's use an example that is more suitable. We assume the following data come from an experiment where we want to find the (locally) maximal value. We see that the interpolated curve has a maximum close to the maximum of the underlying ground truth (which would be the maximum 1 at $x = π$ ). This example uses the command Optimization[Maximize].

ground_truth := x -> 1/(1 + (x - Pi)^2):
xsampled := [seq(i, i = 0. .. 5, 0.7)]:
ysampled := map(ground_truth, xsampled):

interpolated := Interpolate(xsampled, ysampled):
optimum := Optimization[Maximize](unapply(interpolated(x), x, numeric), 0 .. 5);

plots:-display(plot([ground_truth, interpolated], 0 .. 5, linestyle=["dash", "solid"]),
pointplot(<xsampled, ysampled>, symbolsize=10, color=blue),
pointplot([optimum[2][1], optimum[1]], symbolsize=25, color=green, symbol=soliddiamond));

$[0.9635817101770758, [\begin{array}{c} 3.1369851865331473 \end{array}]]$

Efficiency

The object returned by the Interpolate command uses the command CurveFitting[ArrayInterpolation] behind the scenes, for the methods that are listed on the Interpolate help page as accepting 1-dimensional data only. There is some overhead in setting up each individual call to that command; the Interpolate object minimizes this overhead, but there is significant overhead inherent in the model where there is a call to ArrayInterpolation for each value computed. Conversely, if you can arrange to compute many values with a single call, this is more efficient. We will illustrate this below, reusing the first example set of data, and using the command CodeTools:-Usage to print information about memory and time resources used.

A relatively large number of 5000 interpolation points will be used in order to more clearly demonstrate some performance differences, and each example will be run 10 times. This large number of points produces some jaggedness in the plot, but we're not primarily concerned with the appearance of the resulting plot. We create this plot in three different ways:

1.	With a single direct call to the ArrayInterpolation command, on all 5000 points at once.

2.	Using the result from Interpolate, f, to compute all 5000 points at once; this uses a calling sequence we haven't seen yet in this document, where you submit a Vector of data points to f.

3.	Letting the plot command call f, with a separate call for each point.

>	interpolationpoints := Vector(5000, i -> i*15/5000, datatype=float[8]):

>	CodeTools:-Usage(plot(<interpolationpoints \| ArrayInterpolation(data2D, interpolationpoints, method=spline)>), iterations=10);

memory used=3.26MiB, alloc change=32.00MiB, cpu time=28.80ms, real time=29.20ms, gc time=5.60ms

>	CodeTools:-Usage(plot(<interpolationpoints \| f(interpolationpoints)>), iterations=10);

memory used=3.23MiB, alloc change=0 bytes, cpu time=22.40ms, real time=22.70ms, gc time=0ns

>	CodeTools:-Usage(plot(f, 1 .. 15, sample=convert(interpolationpoints, list), adaptive=false), iterations=10);

memory used=89.99MiB, alloc change=8.00MiB, cpu time=869.60ms, real time=811.40ms, gc time=118.80ms

We see that the first two options are about equally fast, and the third is 30-40 times slower. So: it is very advantageous for efficiency if one can pick a large number of points to interpolate at, at once.

In order to minimize overhead, the object created by Interpolate does the following:

- Upon construction of the object, it preprocesses the input data, by converting it to a floating-point rtable and sorting it. This means that this happens only once, rather than with each call.

- The object pre-allocates a single element rtable with the correct data type for the result. This prevents any allocation from taking place for calls that seek to compute a single value.

Obtaining a closed-form expression

The above methods can be used for any application where you need to evaluate the interpolant numerically; it cannot, however, produce a closed-form (typically piecewise) expression for the interpolating function. Such closed-form expressions can be obtained from some commands in the CurveFitting package.

>	Spline(xdata, ydata, v);

$\{\begin{array}{c} 6.40737377947708 - 3.29095019314989 v + 1.38269059814989 {(v - 1)}^{3} & v < 2 \\ - 0.506079211218933 + 0.857121601299780 v + 4.14807179444967 {(v - 2)}^{2} - 2.86512326974945 {(v - 2)}^{3} & v < 3 \\ 1.67454797453381 + 0.557895380950769 v - 4.44729801479868 {(v - 3)}^{2} + 1.92150000984791 {(v - 3)}^{3} & v < 4 \\ 11.6691339690577 - 2.57220061910286 v + 1.31720201474505 {(v - 4)}^{2} + 0.988945535357804 {(v - 4)}^{3} & v < 5 \\ - 14.0309216587970 + 3.02904001646066 v + 4.28403862081847 {(v - 5)}^{2} - 3.81310060227913 {(v - 5)}^{3} & v < 6 \\ 3.66736375129320 + 0.157815451260210 v - 7.15526318601892 {(v - 6)}^{2} + 4.56321490775871 {(v - 6)}^{3} & v < 7 \\ 5.42148701435834 - 0.463066197501500 v + 6.53438153725721 {(v - 7)}^{2} - 3.98741633475571 {(v - 7)}^{3} & v < 8 \\ - 0.883660344840631 + 0.643447872745792 v - 5.42786746700992 {(v - 8)}^{2} + 2.95974476326412 {(v - 8)}^{3} & v < 9 \\ 14.4367227496527 - 1.33305277148167 v + 3.45136682278246 {(v - 9)}^{2} - 1.79965285430079 {(v - 9)}^{3} & v < 10 \\ 1.05068589165074 + 0.170722311180876 v - 1.94759174011991 {(v - 10)}^{2} + 1.21867068593904 {(v - 10)}^{3} & v < 11 \\ 2.95265048389397 - 0.0684491112418367 v + 1.70842031769720 {(v - 11)}^{2} - 0.717063285455365 {(v - 11)}^{3} & v < 12 \\ - 11.2438018327704 + 1.19720166778647 v - 0.442769538668893 {(v - 12)}^{2} + 0.313852219882421 {(v - 12)}^{3} & v < 13 \\ - 12.1009477207848 + 1.25321925009595 v + 0.498787120978371 {(v - 13)}^{2} - 1.00162767607432 {(v - 13)}^{3} & v < 14 \\ 15.4985347316823 - 0.754089536170272 v - 2.50609590724459 {(v - 14)}^{2} + 0.835365302414864 {(v - 14)}^{3} & otherwise \end{array}$

(1.2.1)

2- and higher-dimensional case

Below we'll distinguish between data interpolation, which matches the given data points exactly, and data smoothing, which approximates noisy data.

Smoothing

In this section, we are assuming that we have a collection of 2- or higher-dimensional independent data for which the 1-dimensional dependent data has some error component or is known to be noisy. The presence of noise implies that whatever surface this data represents does not necessarily pass through the given dependent data points.

In this case, a surface is approximated by numerically smoothing the data using the lowess algorithm. For any given point to be plotted a window of close enough, surrounding data points will be used to compute a local, weighted, low order fit.

Scatterplot3D

The command ScatterPlot3D from the Statistics package provides a smoothed plot of the 2D noisy data. The independent data does not need to be regularly spaced, and is supplied as an n-by-3 Array or Matrix. Each of the n rows represents an individual point. The columns are the x-, y-, and z-values.

Here is an example. First, we'll look down upon the projection of the data points onto the x-y plane, and thus visualize the layout and spacing of the 2D dependent data values.

>	X := Sample(Uniform(-50,50),175):

>	Y := Sample(Uniform(-50,50),175):

>	Zerror := Sample(Normal(0,100),175):

>	Z := Array(1..175,(i)->-(sin(Y[i]/20)*(X[i]-6)^2+(Y[i]-7)^2+Zerror[i])):

>	XYZ := Matrix([[X],[Y],[Z]],datatype=float[8])^%T;

>	ScatterPlot3D(XYZ, axes=box, orientation=[20,0,0]);

>	ScatterPlot3D(XYZ, lowess, grid=[25,25], axes=box, orientation=[20,70,0]);

Here is another example, which reads in the n-by-3 data from a file.

>	M:=ImportMatrix(cat(kernelopts(mapledir),"/data/plotting/irregulardata3D.csv"), source=csv,datatype=float[8]):

>	Ppt := pointplot3d(M,axes=box,style=patchnogrid,symbolsize=10,color=red, labels=[x,y,z]):

>	display(Ppt,orientation=[90,0,0]);

The default behavior for the ScatterPlot3D command is to use an order 2 quadratic form, and to use a window of width 1/3 of the original domain in each of the x- and y-directions.

On the right below is the smoothed surface, while on the left is a slightly transparent view of the same surface overlaid with a point-plot of the original data points.

>	Ploess := ScatterPlot3D( M, lowess ):

>	display(Array([ display([Ppt,Ploess], transparency=0.2), display(Ploess,style=patch,color=RGB(0.0,0.4,0.6)) ]),view=700..860,axes=box,orientation=[54,78,-15], transparency=0.0);

Lowess

If you need to do more than just plot the surface, you can use the CurveFitting[Lowess] command. It returns a procedure that runs the lowess algorithm. This procedure can deal with arbitrary-dimensional input. Of course, if the input points are 2-dimensional, one can use this procedure for a plot, as well.

>	p := Lowess(M):

p(0, 0);

$841.914000008457$

(2.1.2.1)

>	optimum := Optimization:-Maximize(p(x, y), x=-6000 .. 6000, y=-6000 .. 6000);

$optimum ≔ [849.801754007939735, [x = 1865.27968390193, y = −1327.90514681541]]$

(2.1.2.2)

>	plots:-display( plot3d(p, -6000 .. 6000, -6000 .. 6000), pointplot3d([op(eval([x, y], optimum[2])), optimum[1]], symbolsize=25, color=red, orientation=[54,78,-15]));

Below, we plot the (numerically estimated) derivative of the smoothed surface in the x- and y-direction, in red and blue, respectively.

>	plot3d([fdiff(p, [1], [x, y]), fdiff(p, [2], [x, y])], x = -6000 .. 6000, y = -6000 .. 6000, color=[red, blue]);

Interpolation

In this section we are assuming that we have a collection of 2- or higher-dimensional independent data for which the 1-dimensional dependent data is accepted as correct. The goal is to produce a surface which must pass through all of the data points. That is, at every x-y point of the independent data the height of the plotted surface must match the corresponding value of the dependent (z) data.

In this case, a surface is computed numerically by interpolating the data.

This falls into one of two distinct situations. The first case has all the independent x-y data forming a regular grid, and the second case consists of irregularly spaced independent data points in the x-y plane. In both cases, the surface may be plotted using the surfdata command.

Data points in a grid

In the case of a uniform 2-dimensional grid of data a surface plot can be generated using the surfdata command.

The data is comprised of a grid of x- and y-points. For the examples of this section, the data points are taken uniformly in each direction. The x-values are taken as the integers from 1 to 7, and the y-values are taken as the integers from 1 to 9. (Non-uniform grid data can be processed in the same way, but for non-grid data, see the later sections.)

View of data points in the x-y plane

Here is a view of the input data points, which in this example is a full grid, uniformly spaced in both the x and y directions.

>	xpts:=<($1..7)>: ypts:=<($1..9)>:

>	display( seq(plottools:-line([xpts[i],1],[xpts[i],9]), i=1..7), seq(plottools:-line([1,ypts[j]],[7,ypts[j]]), j=1..9), pointplot([seq(seq([xpts[i],ypts[j]],j=1..9),i=1..7)], color=red) );

The given data values corresponding to each data point are stored in data which is a 7x9 Array. A random collection of values is generated for this example.

>	data3D:=Array( LinearAlgebra:-RandomMatrix(7, 9, generator=0.2 .. 0.8), datatype=float[8]):

Computations for such data can be done using some commands in the Interpolation package: the commands listed on the help page for SplineInterpolation. These commands call CurveFitting[ArrayInterpolation] behind the scenes. They accept arbitrary-dimensional grid data; in order to be able to show plots, we use 2-dimensional input points.

>	f := SplineInterpolation([xpts, ypts], data3D);

$Interpolation :- SplineInterpolation ([[\begin{array}{c} 1.0 \\ 2.0 \\ 3.0 \\ 4.0 \\ 5.0 \\ 6.0 \\ 7.0 \end{array}], [\begin{array}{c} 1.0 \\ 2.0 \\ 3.0 \\ 4.0 \\ 5.0 \\ 6.0 \\ 7.0 \\ 8.0 \\ 9.0 \end{array}]], [\begin{array}{c} 0.4394237129718905 & 0.2334864716054194 & 0.3530279668184577 & 0.5811382828416922 & 0.329132906753384 & 0.4257115852204389 & 0.7407807077973092 & 0.6045309713542488 & 0.5830824563358177 \\ 0.6537581716442946 & 0.3283500920129456 & 0.5803027267434254 & 0.49909705039389723 & 0.5344330142446085 & 0.36088329441371786 & 0.4782985287500852 & 0.739674339722364 & 0.4606682853258782 \\ 0.44754345383522176 & 0.37997366187623427 & 0.6282167186778296 & 0.4567297883426742 & 0.39488827113152647 & 0.38565812448548487 & 0.40260018867213376 & 0.43445161412115535 & 0.493520873105826 \\ 0.26813123992450677 & 0.7684096138696237 & 0.2725921057892862 & 0.5356245165007512 & 0.2854508694100452 & 0.3507484652973021 & 0.5119334605079597 & 0.2318917785390107 & 0.7145829159779697 \\ 0.2622267959819816 & 0.740189323829678 & 0.4680381928798836 & 0.2872738022939716 & 0.6146244293662912 & 0.5322672878666747 & 0.7086220044857923 & 0.22497458657142758 & 0.4467450022571856 \\ 0.6487037917023593 & 0.3293460551570467 & 0.28331834999781264 & 0.42989787998858314 & 0.5682950163278608 & 0.6582734443248712 & 0.24373262362755188 & 0.4052123669097832 & 0.511058941604045 \\ 0.3682588366844777 & 0.5216256708793106 & 0.6899911099144254 & 0.4393607107008185 & 0.26521738844863446 & 0.5812471055229301 & 0.39596458852053373 & 0.6142710156337192 & 0.775137710060507 \end{array}], verify = false)$

(2.2.1.1)

>	g := LinearInterpolation([xpts, ypts], data3D);

$Interpolation :- LinearInterpolation ([[\begin{array}{c} 1.0 \\ 2.0 \\ 3.0 \\ 4.0 \\ 5.0 \\ 6.0 \\ 7.0 \end{array}], [\begin{array}{c} 1.0 \\ 2.0 \\ 3.0 \\ 4.0 \\ 5.0 \\ 6.0 \\ 7.0 \\ 8.0 \\ 9.0 \end{array}]], [\begin{array}{c} 0.4394237129718905 & 0.2334864716054194 & 0.3530279668184577 & 0.5811382828416922 & 0.329132906753384 & 0.4257115852204389 & 0.7407807077973092 & 0.6045309713542488 & 0.5830824563358177 \\ 0.6537581716442946 & 0.3283500920129456 & 0.5803027267434254 & 0.49909705039389723 & 0.5344330142446085 & 0.36088329441371786 & 0.4782985287500852 & 0.739674339722364 & 0.4606682853258782 \\ 0.44754345383522176 & 0.37997366187623427 & 0.6282167186778296 & 0.4567297883426742 & 0.39488827113152647 & 0.38565812448548487 & 0.40260018867213376 & 0.43445161412115535 & 0.493520873105826 \\ 0.26813123992450677 & 0.7684096138696237 & 0.2725921057892862 & 0.5356245165007512 & 0.2854508694100452 & 0.3507484652973021 & 0.5119334605079597 & 0.2318917785390107 & 0.7145829159779697 \\ 0.2622267959819816 & 0.740189323829678 & 0.4680381928798836 & 0.2872738022939716 & 0.6146244293662912 & 0.5322672878666747 & 0.7086220044857923 & 0.22497458657142758 & 0.4467450022571856 \\ 0.6487037917023593 & 0.3293460551570467 & 0.28331834999781264 & 0.42989787998858314 & 0.5682950163278608 & 0.6582734443248712 & 0.24373262362755188 & 0.4052123669097832 & 0.511058941604045 \\ 0.3682588366844777 & 0.5216256708793106 & 0.6899911099144254 & 0.4393607107008185 & 0.26521738844863446 & 0.5812471055229301 & 0.39596458852053373 & 0.6142710156337192 & 0.775137710060507 \end{array}], verify = false)$

(2.2.1.2)

>	plot3d([f(x,y), g(x,y) + 0.2], x=1..7, y=1..9, color=[red, blue]);

>	Optimization:-Minimize(f(x, y), x = 1..7, y=1..9);

$[0.236597474643406969, [x = 3.83202160187721, y = 5.29959188777011]]$

(2.2.1.3)

>	fsolve({(a,b)->a-b, (a,b)->f(a,b)-0.5}, [1..7, 1..9]);

$[2.4966859680690522, 2.4966859680690522]$

(2.2.1.4)

>	f(op(%)) - 0.5;

$−7.10542735760100 10^{−15}$

(2.2.1.5)

The plot above uses the normal plot3d command, which does not know to select its sample points at the exact locations where our data points are, because it does not receive the actual data, only the interpolating object. The surfdata command does know to do this, because it is given the actual data.

The first 3D plot below is the piecewise-planar surface produced by default. The goal is to produce a smoother surface instead, and this will be accomplished by using the gridsize and interpolation options of the surfdata command.

The plotted surface is being overlaid here with both the 3D point-plot as well as the (linear) patch-surface, so that the computed surface may be visually demonstrated to be authentic with respect to the original discrete data.

>	ptsplot:=PLOT3D( GRID(1..7, 1..9, data3D, STYLE(POINT) ) ): gridplot:=PLOT3D( GRID( 1..7, 1..9, data3D, STYLE(WIREFRAME)) ):

The following plot3d options control the general look of these plots, and are reused in this subsection. They are assigned to a single, reusable name so that the differences between methods is more clearly illustrated.

>	lookandfeel := axes=boxed, style=surface, glossiness=1.0, lightmodel=light1, color=RGB(204/255, 84/255, 0/255), view=[1..7, 1..9, 0..1], orientation = [60, 51, -6]:

>	display( ptsplot,gridplot, surfdata(data3D, 1..7, 1..9, lookandfeel) );

On the left below is the default surface produced by the surfdata command, without any additional interpolation.

On the right below is the surface produced by the surfdata command using the new interpolating optional keyword parameters gridsize and interpolation.

The gridsize parameter is an equation of the form gridsize=list, where the the list contains a pair of positive integers denoting the size of the uniform grid of points at which interpolation of the original data will be performed.

The interpolation parameter is an equation of the form interpolation=list, where the list contains options accepted by the CurveFitting:-ArrayInterpolation command.

>	display( Array([ surfdata( data3D, 1..7, 1..9, lookandfeel ), surfdata( data3D, 1..7, 1..9, gridsize=[28,36], interpolation=[method=spline], lookandfeel) ]) );

display( Array([

     display( ptsplot, gridplot,
              surfdata( data3D, 1..7, 1..9, lookandfeel ) ),

     display( ptsplot, gridplot,
              surfdata( data3D, 1..7, 1..9,
                        gridsize=[28,36], interpolation=[method=spline],
                        lookandfeel) )

                 ]) );

Note that the 3D plot renderer does its own small amount smoothing of the surface. Hence, even when using the purely linear method of the computational interpolation scheme, the plot on the right below shows a modest level of surface smoothing.

display( Array([

     display( ptsplot, gridplot,
              surfdata( data3D, 1..7, 1..9, lookandfeel ) ),

     display( ptsplot, gridplot,
              surfdata( data3D, 1..7, 1..9,
                        gridsize=[25,33], interpolation=[method=linear],
                        lookandfeel) )

                 ]) );

The refined grid on which values are interpolated does not need to contain a whole integer multiple (in either direction) of the number points of the original grid of data points.

display( Array([

     display( ptsplot, gridplot,
              surfdata( data3D, 1..7, 1..9, lookandfeel ) ),

     display( ptsplot, gridplot,
              surfdata( data3D, 1..7, 1..9,
                        gridsize=[33,33], interpolation=[method=cubic],
                        lookandfeel) )

                 ]) );

Non-uniform data grid

This is the case of irregularly spaced data. In the case of two-dimensional independent data and one dependent dimension, this is typically supplied as an n-by-3 Array or Matrix. Each row represents a point, which the columns interpreted as x-, y-, and z-value.

Here below is an example. First, we'll look down upon the projection of the data points onto the x-y plane and thus visualize the layout and spacing of the 2D dependent data values.

This is the same data that was used for the ScatterPlot3D example in the Smoothing section. But in contrast to the smoothing example we are now supposing that the data points are to be interpolated, which is to say that the surface must pass directly through the data points.

>	M:=ImportMatrix(cat(kernelopts(mapledir),"/data/plotting/irregulardata3D.csv"), source=csv,datatype=float[8]): points := M[.., 1..2]: values := M[.., 3]:

>	Ppt := pointplot3d(M,axes=box,style=patchnogrid,symbolsize=15,color=red, labels=[x,y,z]):

>	display(Ppt,orientation=[90,0,0]);

There are many different approaches to working with such data; we will examine the five that are implemented in Maple. The first two methods, natural neighbor interpolation and linear triangular interpolation, determine a so-called Delaunay triangulation of the convex hull of the input points, and use this triangulation to determine which points contribute to the interpolated value. For visualizing such a Delaunay triangulation, one can use the surfdata command.

The option `source=irregular` is supplied to denote that this data is not to be interpreted as forming a regularly spaced grid with only three values along one of the two independent data dimensions.

>	sd := surfdata(M,source=irregular,axes=box, thickness=0): plots:-display(Ppt, sd);

If you want to do computations with the interpolant, rather than just generating a plot, you can use the linear triangular method of interpolation. The interpolant is exactly the surface you see in the plot. You can obtain this interpolant from the Interpolation[Interpolate] command by giving it the option $method = lineartriangular$ , or equivalently by calling the LinearTriangularInterpolation command in the same package. The interpolant is undefined at points outside the convex hull of the input points.

>	f_lt := Interpolate(points, values, method = lineartriangular);

$f_lt ≔ (Linear Triangular interpolation object with 176 sample points)$

(2.2.2.1)

>	f_lt(-4000, 2000);

$790.016559118906$

(2.2.2.2)

Natural neighbor interpolation is the default method employed by Interpolation[Interpolate] for two-dimensional input points (so the problem is three-dimensional if you include the values). You can get the same object by calling the NaturalNeighborInterpolation command in the same package. This interpolation method is based on the Delaunay triangulation shown above, but leads to a smoother interpolant than the linear triangular method. As a consequence of its dependence on the Delaunay triangulation, the interpolated value is undefined if the sample point is not within the convex hull of the input points. This, combined with the sampling grid of the $plot3d$ command, causes the jagged edges in the plot below.

>	f_nn := Interpolate(points, values);

$f_nn ≔ (Natural Neighbor interpolation object with 176 sample points)$

(2.2.2.3)

>	f_nn(-4000, 2000); plots:-display(Ppt, plot3d(f_nn, -6100 .. 6000, -6600 .. 6600, grid=[149,149]));

$790.114010448548$

The previous two implementation methods are implemented in Maple only for two-dimensional input points. The radial basis function, inverse distance weighted, and Kriging interpolation methods are implemented for any dimension - but here we will show them with two-dimensional input points.

The radial basis function interpolant is a linear combination of radial basis functions centered at each input point, with the coefficients chosen so that the result is an interpolant. You can use this interpolation method by calling the Interpolation[Interpolate] command with the option $method = radialbasisfunction$ , or equivalently by calling the RadialBasisFunctionInterpolation command in the same package. Several different radial basis functions are available; they are listed on the RadialBasisFunctionInterpolation help page and they can be selected by passing the name of the radial basis function to either of these commands (in some cases with an optional numeric shape parameter).

>	f_rb := Interpolate(points, values, method = radialbasisfunction);

$f_rb ≔ (\begin{array}{c} Ra&DifferentialD;ⅈal Basⅈs Functⅈon ⅈnt&ExponentialE;rpolatⅈon obȷ&ExponentialE;ct wⅈth 176 sampl&ExponentialE; poⅈnts \\ Ra&DifferentialD;ⅈal Basⅈs Functⅈon: multⅈqua&DifferentialD;rⅈc \end{array})$

(2.2.2.4)

>	f_rb_g := Interpolate(points, values, method = radialbasisfunction, inversequadratic, .0008);

$f_rb_g ≔ (\begin{array}{c} Ra&DifferentialD;ⅈal Basⅈs Functⅈon ⅈnt&ExponentialE;rpolatⅈon obȷ&ExponentialE;ct wⅈth 176 sampl&ExponentialE; poⅈnts \\ Ra&DifferentialD;ⅈal Basⅈs Functⅈon: ⅈnv&ExponentialE;rs&ExponentialE;qua&DifferentialD;ratⅈc \end{array})$

(2.2.2.5)

>	f_rb(-4000, 2000); f_rb_g(-4000, 2000); plots:-display(Array(1..2, map(f -> plots:-display(Ppt, plot3d(f, -6100 .. 6000, -6600 .. 6600), view = [DEFAULT, DEFAULT, 700..900]), [f_rb, f_rb_g])));

$791.031389266952488$

$787.067250992000595$

The inverse distance weighted interpolant is obtained by weighting the values associated with all points with the distance to the query point. An extra parameter can be supplied, which is the distance beyond which input points should be disregarded completely. You can use this interpolation method by calling the Interpolation[Interpolate] command with the option $method = inversedistanceweighted$ , or equivalently by calling the InverseDistanceWeightedInterpolation command in the same package

>	f_id := Interpolate(points, values, method=inversedistanceweighted);

$f_id ≔ (\begin{array}{c} Inv&ExponentialE;rs&ExponentialE; Dⅈstanc&ExponentialE; W&ExponentialE;ⅈght&ExponentialE;&DifferentialD; ⅈnt&ExponentialE;rpolatⅈon obȷ&ExponentialE;ct wⅈth 176 sampl&ExponentialE; poⅈnts \\ Ra&DifferentialD;ⅈus of ⅈnflu&ExponentialE;nc&ExponentialE;: ⅈnfⅈnⅈty \end{array})$

(2.2.2.6)

>	f_id_2000 := Interpolate(points, values, method=inversedistanceweighted, 2000);

$f_id_2000 ≔ (\begin{array}{c} Inv&ExponentialE;rs&ExponentialE; Dⅈstanc&ExponentialE; W&ExponentialE;ⅈght&ExponentialE;&DifferentialD; ⅈnt&ExponentialE;rpolatⅈon obȷ&ExponentialE;ct wⅈth 176 sampl&ExponentialE; poⅈnts \\ Ra&DifferentialD;ⅈus of ⅈnflu&ExponentialE;nc&ExponentialE;: 2000. \end{array})$

(2.2.2.7)

>	f_id(-4000, 2000); f_id_2000(-4000, 2000); plots:-display(Array(1..2, map(f -> plots:-display(Ppt, plot3d(f, -6100 .. 6000, -6600 .. 6600, grid=[149, 149])), [f_id, f_id_2000])));

$792.621960368492296$

$786.057624146955391$

Finally, Kriging is an interpolation method based on a statistical model of how distance between two points influences the correlation between the values at those two points. This model is called a variogram. Several parameterized variograms are implemented in Maple, as listed on the Interpolation[Kriging][SetVariogram] page. Variogram parameters can be inferred automatically, but the best results are obtained when the variogram is chosen manually. Sometimes theoretical considerations about the problem, for example from literature research, can give reasons to choose a particular variogram; at other times, you can estimate a variogram from the input points and values. That is what we do in this case. We start by examining the empirical variogram, that is, the graph of distances between pairs of points and absolute differences in value between them. This can be done with the Interpolation[Kriging][DisplayVariogram] command. We can call this command once we have created a Kriging object, using the command Interpolation[Interpolate] with the option $method = kriging$ , or the equivalent command Interpolation[Kriging].

>	f_kr := Interpolate(points, values, method=kriging);

$f_kr ≔ (\begin{array}{c} Krⅈgⅈng ⅈnt&ExponentialE;rpolatⅈon obȷ&ExponentialE;ct wⅈth 176 sampl&ExponentialE; poⅈnts \\ Varⅈogram: Sph&ExponentialE;rⅈcal(81.,1444.,5634.358496) \end{array})$

(2.2.2.8)

>	DisplayVariogram(f_kr, output=groupedpoints);

This shows that the values are generally closer when points are closer, but it is hard to see what the relationship might look like. We tweak the display a bit by grouping distance-difference pairs together and averaging them. That has happened already for the graph above, but we need to do so more aggressively.

>	binned := DisplayVariogram(f_kr, output=bins, bounds=25): binmeans := map(convert, map(Statistics:-Mean, binned), list): binmeans_t := ListTools:-Transpose(binmeans):

>	pointplot(binmeans, view=[0 .. max(binmeans_t[1]), 0 .. max(binmeans_t[2])]);

This suggests that the minimum expected difference is about 100, the maximum is about 1600, and that maximum is reached when the distance between the points is around 8000. This would be reasonably well represented by the variogram $Circular (100, 1600, 8000)$ .

>	SetVariogram(f_kr, Circular(100, 1600, 8000));

$(\begin{array}{c} Krⅈgⅈng ⅈnt&ExponentialE;rpolatⅈon obȷ&ExponentialE;ct wⅈth 176 sampl&ExponentialE; poⅈnts \\ Varⅈogram: Cⅈrcular(100,1600,8000) \end{array})$

(2.2.2.9)

Now that we have set the variogram, we can evaluate the Kriging interpolator at any point, and create a plot of the interpolator.

>	f_kr(-4000, 2000); plots:-display(Ppt, plot3d(f_kr, -6100 .. 6000, -6600 .. 6600, grid=[99,99]));

$790.783008163623322$

The statistics behind the method give us another output: the expected variance associated with a predicted value, which gives a measure for our confidence in that value. We can ask for that for one point, or show it as the surface color in a 3D plot.

>	f_kr(-4000, 2000, output=variance);

$227.482481165330682$

(2.2.2.10)

>	plots:-display(Ppt, plot3d(f_kr(x, y), x=-6100 .. 6000, y=-6600 .. 6600, color = [f_kr(x, y, output=variance)/800, 0.8, 0.8, colortype=HSV]));

We see that areas close to many points give higher confidence than areas that are further away from most points.

	See Also
	For additional details on the commands used in the worksheet, including additional options to control the interpolating and smoothing methods, see the following help pages: Interpolation[Interpolate] Statistics[Lowess] Spline ArrayInterpolation plot plot,options plot3d surfdata ScatterPlot3D display RandomMatrix

Return to Index for Example Worksheets

Download Help Document

Maple

Maple 附加模块

MapleSim

MapleSim 附加模块

系统工程

项目服务

Maple T.A. and Möbius

教育

行业

汽车与航空航天

机器人

机器设计和工业自动化

其他

应用领域

产品价格

购买

院系/全校正版授权

Maplesoft精英维护升级计划

支持

产品培训

产品在线帮助

研讨会与活动

出版

资源中心

示例和应用

社区

关于 Maplesoft

媒体中心

用户社区

联系我们

Online Help

All Products Maple MapleSim

Maple

数学软件

Maple 附加模块

MapleSim

多学科系统级建模仿真

MapleSim 附加模块

系统工程

项目服务

Maple T.A. and Möbius

教育

行业

汽车与航空航天

机器人

机器设计和工业自动化

其他

应用领域

产品价格

购买

院系/全校正版授权

Maplesoft精英维护升级计划

支持

产品培训

产品在线帮助

研讨会与活动

出版

资源中心

示例和应用

社区

关于 Maplesoft

媒体中心

用户社区

联系我们

Online Help

All Products Maple MapleSim