we have
$$
u_x=f_rr_x, u_{xx}=f_rr_{xx}+r_x(f_{rr}r_x)\text{ similarly for } y,
$$
$$
u_{xx}+u_{yy}=f_rr_{xx}+r_x(f_{rr}r_x)+f_rr_{yy}+r_y(f_{rr}r_y)
=f_{rr}(r_x^2+r_y^2)+f_r(r_{xx}+r_{yy}).
$$
with
$$
r_x=\frac{x}{\sqrt{x^2+y^2}}, r_{xx}=\frac{y^2}{(x^2+y^2)^{3/2}} \text{ similarly for } y
$$
we have
$$
u_{xx}+u_{yy}=f_{rr}+\frac{1}{r}f_r.
$$
of course if $f$ depends on $\theta$ it gets more complicated
Here is my proof:
$$ \frac{\partial}{\partial x}=\frac{\partial}{\partial r} \frac{\partial r}{\partial x} +\frac{\partial}{\partial \theta}\frac{\partial \theta}{\partial x}$$ so by applying the product rule
$$ \frac{\partial^2}{\partial x^2}=\frac{\partial}{\partial x}\left(\frac{\partial}{\partial r} \frac{\partial r}{\partial x}+\frac{\partial}{\partial \theta}\frac{\partial \theta}{\partial x}\right)=
\frac{\partial^2}{\partial r^2}\left(\frac{\partial r}{\partial x}\right)^2+
\frac{\partial}{\partial r}\frac{\partial^2 r}{\partial x^2} +
\frac{\partial^2}{\partial \theta^2}\left(\frac{\partial \theta}{\partial x}\right)^2+
\frac{\partial}{\partial \theta}\frac{\partial^2 \theta}{\partial x^2}$$
$$ +\frac{\partial^2 }{\partial r \partial \theta}\frac{\partial r}{\partial x} \frac{\partial \theta}{\partial x}$$
You get the same relation with $y$. Now you just have to calculate the derivatives of $r,\theta$ with respect to $x,y$. You have
$$ r=\sqrt{x^2+y^2},\ \theta=\arctan \frac{y}{x}$$.
What follows is a simple calculus exercise on derivatives. You just need to prove that
$$\left(\frac{\partial r}{\partial x}\right)^2+\left(\frac{\partial r}{\partial y}\right)^2=1, $$
etc.(the relations you need so that when you sum $\frac{\partial^2}{\partial x^2}$ and $\frac{\partial^2}{\partial y^2}$ you'll get the polar form of the laplacian).
Going on the lines you started, you shouldn't change $\sin \theta,\cos \theta$ in terms of $x,y$. Just calculate the next derivative with respect to $r,\theta$, using appropriately the formula of the partial derivative of composition of functions.
Best Answer
The meaning of $dA=r dr d\theta$ is that when you integrate over a region with respect to $dA$, you get the same result as you do by parametrizing it in polar coordinates and integrating with respect to the differential $r dr d \theta$.
Thus your question (because $d\theta$ is not a problem) reduces to showing that "integrating" against a differential like "$(dr)^2$" gives zero. What this actually means is that sums of the form $S=\sum_i f(r_i) (\Delta r)_i^2$ should all go to zero as you refine the partition.
For bounded $f$ for example $|S| \leq ML\delta$ where $M$ is the bound on $|f|$, $L$ is the total length of the interval, and $\delta$ is the longest of the $(\Delta r)_i$'s. This indeed goes to zero as $\delta \to 0$.
The point can be understood without inspecting all the details by just counting how many terms there are (call that $n$) and how small the individual terms are relative to that (which is on the order of $1/n^2$).