Simple proof of integration in polar coordinates?

Of course, if you break $\mathbb{R}^2$ into a polar grid

$\hspace{3.5cm}$

the small slightly curved rectangles have area $r\,\mathrm{d}\theta\,\mathrm{d}r$.

However, it seems that you are interested in looking at $$ \begin{align} \mathrm{d}y\,\mathrm{d}x &=(\sin(\theta)\,\mathrm{d}r+r\cos(\theta)\,\mathrm{d}\theta)(\cos(\theta)\,\mathrm{d}r-r\sin(\theta)\,\mathrm{d}\theta)\\ &=r\,\mathrm{d}\theta\,\mathrm{d}r \end{align} $$ and why the $\mathrm{d}r^2$ and $\mathrm{d}\theta^2$ terms disappear and the $\mathrm{d}r\,\mathrm{d}\theta$ and $\mathrm{d}\theta\,\mathrm{d}r$ have different signs.

Let's start with $$ \begin{align} \mathrm{d}x&=\cos(\theta)\,\mathrm{d}r-r\sin(\theta)\,\mathrm{d}\theta\\ \mathrm{d}y&=\sin(\theta)\,\mathrm{d}r+r\cos(\theta)\,\mathrm{d}\theta \end{align} $$ rewritten as $$ \begin{bmatrix}\mathrm{d}x\\\mathrm{d}y\end{bmatrix} =\begin{bmatrix}\cos(\theta)\\\sin(\theta)\end{bmatrix}\mathrm{d}r +\begin{bmatrix}-r\sin(\theta)\\r\cos(\theta)\end{bmatrix}\mathrm{d}\theta $$ Therefore, the displacements $\color{green}{\mathrm{d}r}$ and $\color{red}{\mathrm{d}\theta}$ get mapped to $\color{green}{\begin{bmatrix}\cos(\theta)\\\sin(\theta)\end{bmatrix}\mathrm{d}r}$ and $\color{red}{\begin{bmatrix}-r\sin(\theta)\\r\cos(\theta)\end{bmatrix}\mathrm{d}\theta}$ in $\mathbb{R}^2$:

$\hspace{3cm}$

where the area in gray is given by $\color{green}{\begin{bmatrix}\cos(\theta)\\\sin(\theta)\end{bmatrix}\mathrm{d}r}\times\color{red}{\begin{bmatrix}-r\sin(\theta)\\r\cos(\theta)\end{bmatrix}\mathrm{d}\theta}=r\,\mathrm{d}r\,\mathrm{d}\theta$.

The fact that the cross product is involved is the reason that the $\mathrm{d}r^2$ and $\mathrm{d}\theta^2$ terms disappear and the $\mathrm{d}r\,\mathrm{d}\theta$ and $\mathrm{d}\theta\,\mathrm{d}r$ have different signs. This, and its $n$-dimensional analogs, are why we use wedge products and differential forms when changing variables.

The area element needs to be computed carefully.

This is very informal, but perhaps you should think of the volume element as a pair, rather than just a product, as in: $$\binom{dx}{dy} = \begin{bmatrix} \cos \theta & - r \sin \theta \\ \sin \theta & r \cos \theta \end{bmatrix}\binom{dr}{d\theta}$$ Multiplying a set by a matrix $A$ corresponds to changing the volume by a factor $\det A$. In this case, $\det A = r$, so the volume element computation becomes, informally, $dxdy = r dr d\theta$.

Addendum: The two volume elements are, informally, $[r,r+dr]\times[\theta,\theta+d\theta]$, and $[x,x+dx]\times[y,y+dy]$.

The short answer is that you need to consider the wedge product between the differentials, not a symmetric product as you have written. The reason for this is addressed in Robjohn's excellent answer.

The volume element, $dV$, is formally given by the wedge product of $dx_1, \dots, dx_n$. This means that in $\mathbb{R}^2$ the volume element in Cartesian coordinates should technically be written as $$ dV = dx \wedge dy. $$ Note that the wedge product is antisymmetric, which means that in particular $dx \wedge dy = - dy \wedge dx.$ Taking this into consideration, if you perfrom this wedge product between the $dx$ and $dy$ you calculated, we have the following $$\begin{align*} dV &= dx \wedge dy \\ &= (\cos\theta ~dr - r\sin\theta ~d\theta) \wedge (\sin\theta ~dr + r\cos\theta~ d\theta) \\ &= \cos\theta \sin\theta ~dr\wedge dr + r \cos^2\theta ~dr \wedge d\theta - r\sin^2\theta~ d\theta \wedge dr -r \sin\theta\cos\theta~ d\theta \wedge d\theta\\ &= r(\cos^2\theta + \sin^2\theta) dr \wedge d\theta\\ &= r~ dr \wedge d\theta, \end{align*} $$ where we have also used the fact that $dr \wedge dr = d\theta \wedge d\theta = 0.$ If you calculate the determinant of the Jacobian you'll find $$\det\left(\dfrac{\partial(x,y)}{\partial{(r,\theta)}}\right) = r,$$ which concurs with the above calculation.