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?
$\endgroup$ 11 Answer
$\begingroup$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
