Estimating Partials on discrete data
MAT329, Fall, 2015
Andy Long

Let’s revisit the swimming pool problem (#6, p. 1005), and think about the grid of data that we need to solve it:

data=-{{2,3,4,6,7,8,8},
{2,3,4,7,8,10,8},
{2,4,6,8,10,12,10},
{2,3,4,5,6,8,7},
{2,2,2,2,3,4,4}};
ListPlot3D[data,Mesh->None,InterpolationOrder->0,ColorFunction->"SouthwestColors"]
ListPlot3D[data,Mesh->5,InterpolationOrder->1,ColorFunction->"SouthwestColors"]

One strategy is to use the simple rectangular rules, of which four stand out (corresponding to the right- and left-rectangle rules). We can think of putting a “mask” over the data, where some values get a weight of 0 (those along part of the boundary) and some get a weight of 1: then we might average these estimates, to get a better estimate:

weights00={
{1,1,1,1,1,1,0},
{1,1,1,1,1,1,0},
{1,1,1,1,1,1,0},
{1,1,1,1,1,1,0},
{0,0,0,0,0,0,0}};
weights00*data
estimate00=Total[Total[weights00*data]]
weights01={
{0,0,0,0,0,0,0},
{1,1,1,1,1,1,0},
{1,1,1,1,1,1,0},
{1,1,1,1,1,1,0},
{1,1,1,1,1,1,0}};
weights01*data
estimate01=Total[Total[weights01*data]]
weights10={
{0,1,1,1,1,1,1},
{0,1,1,1,1,1,1},
{0,1,1,1,1,1,1},
{0,1,1,1,1,1,1},
{0,0,0,0,0,0,0}};
weights10*data
estimate10=Total[Total[weights10*data]]
weights11={
{0,0,0,0,0,0,0},
{0,1,1,1,1,1,1},
{0,1,1,1,1,1,1},
{0,1,1,1,1,1,1},
{0,1,1,1,1,1,1}};
weights11*data
estimate11=Total[Total[weights11*data]]

estimate=(estimate00+estimate01+estimate10+estimate11)/4
volume=Abs[estimate]*5*5
weights=(weights00+weights01+weights10+weights11)

Now I want to illustrate how we use a grid of 9 values to approximate all the partial derivatives needed to generate a quadratic approximation to a function; here are the partials.

DiscreteFxApprox[data_,m_,n_,h_]:=(data[[m+1]][[n]]-data[[m-1]][[n]])/(2*h);
DiscreteFyApprox[data_,m_,n_,k_]:=(data[[m]][[n+1]]-data[[m]][[n-1]])/(2*k);
DiscreteFxxApprox[data_,m_,n_,h_]:=(data[[m+1]][[n]]-2data[[m]][[n]]+data[[m-1]][[n]])/(2*h);
DiscreteFyyApprox[data_,m_,n_,k_]:=(data[[m+1]][[n]]-2data[[m]][[n]]+data[[m-1]][[n]])/(2*k);
DiscreteFxyApprox[data_,m_,n_,h_,k_]:=(data[[m+1]][[n+1]]-data[[m+1]][[n-1]]-data[[m-1]][[n+1]]+data[[m-1]][[n-1]])/(4*h*k);

Now I’ll use these on “panels” of nine points, and move them across the data to create quadratic approximations (but these won’t fit -- interpolate -- the data). 

One of our standard strategies in creating integration schemes is to fit data with simple functions, and then integrate them. That’s what we are doing here: fitting a quadratic to the data, and then we will integrate it.

ix=2
iy=4
p1=Plot3D[data[[ix]][[iy]]
+DiscreteFxApprox[data,ix,iy,1](x-ix)
+DiscreteFyApprox[data,ix,iy,1](y-iy)
+.5*DiscreteFxxApprox[data,ix,iy,1](x-ix)^2
+.5*DiscreteFyyApprox[data,ix,iy,1](y-iy)
+ DiscreteFxyApprox[data,ix,iy,1,1](x-ix)(y-iy),
{x,ix-1,ix+1},
{y,iy-1,iy+1}];
chooseData[ix_,iy_]:=data[[ix]][[iy]]
p2=ListPointPlot3D[{
{ix-1,iy-1,chooseData[ix-1,iy-1]},
{ix,iy-1,chooseData[ix,iy-1]},
{ix+1,iy-1,chooseData[ix+1,iy-1]},
{ix-1,iy,chooseData[ix-1,iy]},
{ix,iy,chooseData[ix,iy]},
{ix+1,iy,chooseData[ix+1,iy]},
{ix-1,iy+1,chooseData[ix-1,iy+1]},
{ix,iy+1,chooseData[ix,iy+1]},
{ix+1,iy+1,chooseData[ix+1,iy+1]}
},PlotStyle->PointSize[.1]];
Show[p1,p2,PlotRange->All]

plot[ix_,iy_]:=Plot3D[data[[ix]][[iy]]
+DiscreteFxApprox[data,ix,iy,1](x-ix)
+DiscreteFyApprox[data,ix,iy,1](y-iy)
+.5*DiscreteFxxApprox[data,ix,iy,1](x-ix)^2
+.5*DiscreteFyyApprox[data,ix,iy,1](y-iy)
+ DiscreteFxyApprox[data,ix,iy,1,1](x-ix)(y-iy),
{x,ix-1,ix+1},
{y,iy-1,iy+1}]
plot[2,3]

We can creat multiple, overlapping quadratics over the region, and then integrate each, and combine them (somehow).

plots=Table[plot[ix,iy],{ix,2,4,1},{iy,2,6,1}];
Show[plots,PlotRange->All]

plots=Table[plot[ix,iy],{ix,2,4,2},{iy,2,6,2}];
Show[plots,PlotRange->All]