You can look at the example worksheet on interpolation and smoothing ( http://www.maplesoft.com/support/help/Maple/view.aspx?path=examples%2fInterpolation_and_Smoothing )
At bottom, I've also included a simple method using inverse-distance and weighted averaging (which is fast, but which extrapolates poorly).
Sorry, for some reason I am unable to attach images of the plots right now.
I have used your third column as the dependent data, but it should be straightforward to change that. Just make sure that the dependent data is the 3rd column of the Matrix formed like, say, <P|T|M>.
restart:
P:=Vector([19.02,18.64,18.54,20.5,31.07,32.64,30.01,28.04,40.46,
40.34,39.36,39.18,50.17,51.03,51.39,52.65,51.14,60.66,
59.11,59.94,58.1,58.27,57.74],
datatype=float[8]):
T:=Vector([311.15,307.05,304.65,301.15,298.15,304.64,311.15,300.23,
298.5,304.65,308.6,310.81,300.75,302.63,304.55,311.2,
308.45,311.75,308.15,304.65,302.95,301.85,300.65],
datatype=float[8]):
M:=Vector([6.83,5.94,6.33,5.56,4.2588,4.4338,4.4338,4.2588,
0.684,1.1050,0.82132,0.88382,1.3934,0.57132,
0.92228,1.4588,0.97997,1.4878,1.3838,1.6723,
1.1765,1.3838,1.4588],datatype=float[8])*1E8:
First, one can construct the basic 3D point plot from the data. Some version of the point plot (with some shading scheme) may or may not be displayed alongside the various surfaces constructed below.
plots:-pointplot3d( <P|T|M>, symbolsize=20, axes=box );
Next, a simple triangular mesh can be used to interpolate linearly between the points.
Ptriang := plots:-surfdata( <P|T|M>, view=0..7e8, source=irregular ):
plots:-display( Ptriang,
plots:-pointplot3d(<P|T|M>,symbolsize=20,color=red),
axes=box );
Another approach is to smooth the data. This may not produce great results for such a small amount of data. There are several adjustable options for the method. It is not fast.
Ploess := Statistics:-ScatterPlot3D( <P|T|M>, view=0..7e8,
lowess, fitorder=2,
bandwidth=1/2, rule=1,
grid=[50,50], axes=box ):
plots:-display( Ploess,
plots:-pointplot3d(<P|T|M>,symbolsize=20,color=red) );
Another possible approach is to to use weighted average (weighted by distance or some metric). Below, a black-box procedure ff is created which takes any p-t point and returns the computed m-value. (A better scheme for this would be interpolation over Voronoi cells using so-called natural neighbors and radial basis functions, but I have only unfinished code for that at present.)
f := proc(x::float, y::float,
N::integer,
X::Vector(datatype=float[8]),
Y::Vector(datatype=float[8]),
Z::Vector(datatype=float[8]),
p::integer, R::float)
local i::integer,j::integer,res::float,innerres::float,dist::float;
innerres:=0.0;
for j from 1 to N do
dist:=sqrt((x-X[j])^2+(y-Y[j])^2);
innerres:=innerres+(max(0, abs(R-dist))/(R*dist))^p;
end do;
res:=0.0;
for i from 1 to N do
dist:=sqrt((x-X[i])^2+(y-Y[i])^2);
res := res + Z[i]*(max(0, abs(R-dist))/(R*dist))^p/innerres;
end do;
res;
end proc:
try
ff:=Compiler:-Compile(f);
catch:
ff:=proc(x,y,N,P,T,M,p,R) evalhf(f(x,y,N,P,T,M,p,R)); end proc;
end try:
We can query procedure ff at any (p,t) point.
ff(40.2, 304.7, 23, P, T, M, 3, 4.0);
which produces,
8
1.10586391484649837 10
And ff can now be used directly in the plot3d command,
Pinvdist := plot3d('ff'(p, t, 23, P, T, M, 3, 20.0),
p=min(P)..max(P), t=min(T)..max(T),
numpoints = 900, view=0..7e8):
plots:-display( Pinvdist,
plots:-pointplot3d(<P|T|M>,symbolsize=20,color=red),
axes=box );
