Generalizing Univariate Integration Techniques
MAT329, Fall, 2015
Andy Long

I want to think about generalizing the Trapezoidal rule, to start. It can be considered a method if integrating by fitting points by secant lines -- that is, with linear functions. To what extent can we generalize that?

We begin by remarking that, given four points, it’s not generally possible to fit a plane through them, which would be the first thought, perhaps. Consider the following “panel” of points, which is four points on a rectangular grid (in this case, generated by a function -- but you might think of it as gridded data points generated in any way).

f[x_,y_]:=x^2-y^2
x0=0;y0=1;x1=1;y1=2;
p1=Plot3D[f[x,y],{x,x0,x1},{y,y0,y1}];
p2=ListPointPlot3D[{
{x0,y0,f[x0,y0]},
{x0,y1,f[x0,y1]},
{x1,y0,f[x1,y0]},
{x1,y1,f[x1,y1]}
},PlotStyle->PointSize[Large]]
Show[p1,p2]

The univariate midpoint method can be thought of as fitting a tangent line at the midpoint and integrating that. 

The following illustrates the use of the tangent plane at each of the four points in the panel. 

Which point should we use? Why choose any one? Use them all, and then average them!

Clear[xs,ys,x,y,i,j,xp,yp,x0,x1,y0,y1]
xs={0,1};
ys={1,2};
x0=xs[[1]];
x1=xs[[2]];
y0=ys[[1]];
y1=ys[[2]];
f[x_,y_]:=x^2-y^2;
f[x_,y_]:=x Sin[ Pi/2 x]^2-x y+y^2;
ftan[f_,ix_,iy_]:=
Module[{i=ix,j=iy,xp=xs[[ix]],yp=ys[[iy]]},
fx[x_,y_]=D[f[x,y],x];
fy[x_,y_]=D[f[x,y],y];
p1=Plot3D[{f[x,y],f[xp,yp]+fx[xp,yp](x-xp)+fy[xp,yp](y-yp)},{x,x0,x1},{y,y0,y1}];
p2=ListPointPlot3D[{
{x0,y0,f[x0,y0]},
{x0,y1,f[x0,y1]},
{x1,y0,f[x1,y0]},
{x1,y1,f[x1,y1]}
},PlotStyle->PointSize[Large]];
Show[p1,p2]
]
{ftan[f,1,1],ftan[f,1,2],ftan[f,2,1],ftan[f,2,2]}
(* ,f[xp,yp]+fx[xp,yp](x-xp)+fy[xp,yp](y-yp) *)

If you integrate each of those over the panel, you would get the average of the four points.

Clear[x0,y0,h,k]
x1=x0+h;
y1=y0+k;
mx0=(f10-f00)/h
mx1=(f11-f01)/h
my0=(f01-f00)/k
my1=(f11-f10)/k
I00=Integrate[f00+mx0(x-x0)+my0(y-y0),{x,x0,x1},{y,y0,y1}]
I01=Integrate[f01+mx1(x-x0)+my0(y-y1),{x,x0,x1},{y,y0,y1}]
I10=Integrate[f10+mx0(x-x1)+my1(y-y0),{x,x0,x1},{y,y0,y1}]
I11=Integrate[f11+mx1(x-x1)+my1(y-y1),{x,x0,x1},{y,y0,y1}]
Simplify[I00+I01+I10+I11]

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"]

Mathematica uses interpolation of constant order (Interpolation order 0)  to generate a fit to the data -- but it looks pretty crude! That’s related to the idea of using representative sample values for a Riemann Sum. Notice that they dropped off some values (there is a 4x6 grid of rectangles). 

Mathematica uses interpolation of linear order (Interpolation order 1)  to generate a continuous function that fits the data; one can then use that to integrate. But it’s not really linear. Those are actually little hyperbolic paraboloids, as you might be able to see. That means that they’re quadratic functions. What IS linear is the boundaries. Let’s recreate the Mathematica method:

DiscreteForwardFxApprox[data_,m_,n_,h_]:=(data[[m+1]][[n]]-data[[m]][[n]])/h;
DiscreteForwardFyApprox[data_,m_,n_,k_]:=(data[[m]][[n+1]]-data[[m]][[n]])/k;
plot[ix_,iy_]:=Plot3D[(DiscreteForwardFyApprox[data,ix,iy,1](y-iy)+data[[ix]][[iy]])((ix+1)-x)+(DiscreteForwardFyApprox[data,ix+1,iy,1](y-iy)+data[[ix+1]][[iy]])(x-ix),
{x,ix,ix+1},
{y,iy,iy+1}]

x0=2;
y0=3;
x1=x0+1;
y1=y0+1;
p1=plot[x0,y0];
p2=ListPointPlot3D[{
{x0,y0,data[[x0]][[y0]]},
{x0,y1,data[[x0]][[y1]]},
{x1,y0,data[[x1]][[y0]]},
{x1,y1,data[[x1]][[y1]]}
},PlotStyle->PointSize[.05]];
Show[p1,p2]

What do you get if you integrate a panel? Let’s do that....

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

Now if you have a function and you can’t evaluate it analytically, you may want to approximate it locally by functions which you can integrate (e.g. quadratics). So here is an example of how this might work, where you can think of the points (a,b) as the midpoint of a subinterval.

Clear[f]
f[x_,y_]:=x^2+y^2-x-y+3;
fx[x_,y_]=D[f[x,y],x]
fy[x_,y_]=D[f[x,y],y]
fxx[x_,y_]=D[f[x,y],x,x]
fxy[x_,y_]=D[f[x,y],y,x]
fyy[x_,y_]=D[f[x,y],y,y]

a=1
b=1
TangentBowl[x_,y_]:=f[a,b]+(x-a)(fx[a,b]+0.5*fxx[a,b](x-a))+ (y-b)(fy[a,b]+0.5*fyy[a,b](y-b) + fxy[a,b](x-a))
Plot3D[TangentBowl[x,y],{x,-2,2},{y,-2,2}]
Plot3D[{f[x,y],TangentBowl[x,y]},{x,-2,2},{y,-2,2}]