Two dimensional Numerical integration If I am numerically integrating my function $f(x,y)$ on a two dimensional cartesian grid, say $[0,1]\times[0,1]$ with $\Delta x=\Delta y$ using the values at the center of each cell, what is the order of accuracy, and how do I see this?
 A: Lets at first consider a single square of size $\Delta x \cdot \Delta y$. The area is then
\begin{align}
I = \int_0^{\Delta y} \int_0^{\Delta x} f(x,y)\,dxdy
\end{align}
or if we rewrite it to have a formulation that is centered
\begin{align}
I = \int_0^{\Delta y} \int_0^{\Delta x} f(x,y)\,dxdy =
\int_{-\Delta y/2}^{\Delta y/2} \int_{-\Delta x/2}^{\Delta x/2} f(\Delta x/2+x,\Delta y/2 + y)\,dxdy
\end{align}
Now we apply $2D$ Taylor expansion to $f(\Delta x/2+x,\Delta y/2 + y)$ at expansion point $(\Delta x/2,\Delta y/2 )$ which yields
\begin{align}
f(\Delta x/2+x,\Delta y/2 + y) &=f(\Delta x/2,\Delta y/2 )\\
&+x f_x(\Delta x/2,\Delta y/2 )\\&+yf_y(\Delta x/2,\Delta y/2 )\\
&+O(x^2+y^2)
\end{align}
Now plug this into the integral to get
\begin{align}
I & =
\int_{-\Delta y/2}^{\Delta y/2} \int_{-\Delta x/2}^{\Delta x/2} f(\Delta x/2+x,\Delta y/2 + y)\,dxdy\\
&=\int_{-\Delta y/2}^{\Delta y/2} \int_{-\Delta x/2}^{\Delta x/2} 
f(\Delta x/2,\Delta y/2 )
+x f_x(\Delta x/2,\Delta y/2 )+yf_y(\Delta x/2,\Delta y/2 )
+O(x^2+y^2)\, dx dy \\
&=\int_{-\Delta y/2}^{\Delta y/2} \int_{-\Delta x/2}^{\Delta x/2} 
f(\Delta x/2,\Delta y/2 )\, dx dy\\
&\quad+\int_{-\Delta y/2}^{\Delta y/2} \int_{-\Delta x/2}^{\Delta x/2} 
x f_x(\Delta x/2,\Delta y/2 )+yf_y(\Delta x/2,\Delta y/2 ) dx dy\\
&\quad +\int_{-\Delta y/2}^{\Delta y/2} \int_{-\Delta x/2}^{\Delta x/2}  O(x^2+y^2)\, dx dy \\
&=f(\Delta x/2,\Delta y/2 )\Delta x\Delta y+0+\Delta y\int_{-\Delta x/2}^{\Delta x/2}  O(x^2)\, dx +\Delta x \int_{-\Delta y/2}^{\Delta y/2}  O(y^2)\, dy \\
&=f(\Delta x/2,\Delta y/2 )\Delta x\Delta y +O(\Delta y \Delta x^3+\Delta x\Delta y^3)
\end{align}
Now lets call the point $(\Delta x/2, \Delta y/2)=c$ for center and since $\Delta x = \Delta y $ we have
\begin{align}
I = f(c)\Delta x^2+O(\Delta x^4)
\end{align}
So if you approximate a single cell by its center $f(c)$ you are left with an error of $O(\Delta x^4)$.
But you are interested in the error for the whole domain. This means, that your error is 
\begin{align}
\epsilon = N_x\cdot N_y \cdot O(\Delta x^4) = N_x^2\cdot O(\Delta x^4)
= \frac{1}{\Delta x^2} O(\Delta x^4) = O(\Delta x^2)
\end{align}
So your method is of second order accuracy.
I also wrote a simple Matlab script, to plot the order of convergence which verifies the theoretical analysis.

And here is the Code, if you are interested to play around a little bit.
function quad2d
f = @(x,y) sin(x+y).*exp(y); % <- use whatever function you want
Iexact = integral2(f,0,1,0,1,'AbsTol',1e-10);
M=10;
err =zeros(M,2);
for k = 1:M
    N = 2^k;
    Iapprox=numericalIntegration(f,N);
    err(k,2) =norm(Iexact-Iapprox);
    err(k,1) = N;
end
loglog(err(:,1),err(:,2))
hold on
order1 = err(:,1).^(-1)*err(1,2);
order2 = err(:,1).^(-2)*err(1,2);
loglog(err(:,1),order1,'--r');
loglog(err(:,1),order2,'--g');
grid on
legend('computed error', 'order 1', 'order 2');
end


function Iapprox=numericalIntegration(f,N)
dx = 1/N; dy = dx;
Iapprox = 0.0;
for i=0:N-1
    for j=0:N-1
        Iapprox = Iapprox + dx*dy*f(dx/2+i*dx,dy/2+j*dy);
    end
end
end

