Numerical Integration of a Gaussian Distribution in Polar Coordinates

I want to numerically evaluate a 2D-integral of a specific probability distribution over some given area (I use MATLAB so all the code below is MATLAB code). I broke down the problem so that it appears in the following integration:
\begin{equation}
\int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} \frac{1}{2 \pi \sigma^2} e^{-\frac{x^2+y^2}{2\sigma^2}} dx dy = 1
\end{equation}

Doing this numerically in MATLAB with cartesian coordinates looks something like this:

Dx = 0.05;
Dy = 0.05;

x = -1:Dx:1;
y = -1:Dy:1;

[XX YY] = meshgrid(x,y); 
pdf1 = 1/(2*pi*sigma^2) * exp(- (XX.^2 + YY.^2)/(2*sigma^2));


area1 = sum(sum(pdf1)) * Dx * Dy % numerical integration
error1 = abs(1-area1) % almost no error: 1.1102e-15

However, when I try the same thing in a straightforward manner with polar coordinates, I get the following:

rsteps = 50;
tsteps = 20;


Dr = 1/rsteps; 
Dt = 2*pi/tsteps;

r = [0:1:rsteps]*Dr;
t = [0:1:tsteps-1]*Dt;

[RR TT] = meshgrid(r,t);
pdf2 = RR/(2*pi*sigma^2) .* exp(- (RR.^2)/(2*sigma^2));

area2 = sum(sum(pdf2)) * Dr * Dt % area in radius-phase-plane

error2 = abs(1-area2) % substantial part is missing: 0.0033

At first I thought that I did some error with the conversion to polar coordinates, but letting the value rsteps go to very large number, the error becomes almost zero, but the convergence is pretty bad (e.g., for 50000 steps in $r$-direction, the error is of the order $10^{-9}$ which is way worse than in cartesian coordinates with a much cruder grid).

My intuition says that my straightforward numerical integration fails because I cannot properly evaluate the distribution at the origin ($x=0$, $y=0$) with the polar coordinates, so there is an area element missing (something that would look like a circle around the origin in cartesian coordinates). Is that correct or is there another problem? How can I fix this with a reasonable grid size anyway. Ideally, I would like to use the same grid, and achieve the same error as in cartesian coordinates.

Thanks for your help!

Answer

This is to be expected; in fact we can calculate from the error you report that your $\sigma$ value must have been around $0.1$.

From the way you set up your $r$ grid, you're effectively using the trapezoidal rule. You've got $N+1$ points and you're normalizing by $N$. The two outer points should be weighted by $1/2$, but since the values there are $0$ that doesn't matter. The error of the trapezoidal rule can be estimated by expanding the function at the centre of the trapezoid and integrating the missing quadratic term, which yields

$$\int_{-h/2}^{h/2}\frac{f''(x_0)}2(x-x_0)^2\mathrm dx=\frac{f''(x_0)}{12}h^3\;,$$

where $h$ is the interval length. Approximating the sum over all intervals by an integral then yields

$$
\int_a^b\frac{f''(x_0)}{12}h^3\frac1h\mathrm dx=\frac{h^2}{12}\left(f'(b)-f'(a)\right)=\frac{(b-a)^2}{12N^2}\left(f'(b)-f'(a)\right)\;,
$$

where $N$ is the number of intervals. In your case $b-a=1$, $N=50$, $f'(b)\approx0$ and $f'(a)=1/\sigma^2$ (after cancelling $2\pi$ with the angular integration), so the error would be expected to be about $-(12\cdot50^2\sigma^2)^{-1}$, which is $-0.0033$ for $\sigma\approx0.1$.

The reason you didn't have this problem in the first integration is that you had $f'(a)\approx f'(b)\approx0$ in that case.

I guess the moral of the story is to use quadrature rules with higher error orders.