Rigorous proof that $dx dy=r\ dr\ d\theta$ I get the graphic explanation, i.e. that the area $dA$ of the sector's increment can be looked upon as a polar "rectangle" as $dr$ and $d\theta$ are infinitesimal, but how do you prove this rigorously?

1. Geometrically, the exact area would be $$\frac{(r+dr)^2d\theta}{2} - \frac{r^2d\theta}{2}$$$$= (r + \frac{dr}{2}) dr d\theta  $$$$= r dr d\theta + \frac{dr^2 d\theta}{2}.$$
How do we get rid of $\frac{dr^2 d\theta}{2}$? Is it too insignificant in value compared to $r dr d\theta$ so that it vanishes? If so, is it the reason it gets ignored? And what's more, geometrically speaking, why would the two areas - $dx dy$ and $r dr d\theta$, be equal?
2. Approaching it algebraically, setting $$x = r\sin\theta$$$$y = r \cos\theta$$
gives $$\frac{dx}{d\theta} = r\cos\theta,  dx = r\cos{\theta} d\theta$$$$ \frac{dy}{dr} = \cos\theta,  dy = \cos\theta dr. $$
Multiplying both equations, side by side, gives $$dxdy = r\cos^2\theta dr d\theta.$$
Again I get an extra term, which is $\cos^2 \theta$. In both cases I am unable to derive that $dxdy = rdrd\theta$.
What do I do wrong in my reasoning? It all makes me think that I'm getting something essential terribly wrong.
 A: In the geometric approach, $dr^2=0$ as it is not only small but also symmetric (see here).
In the algebraic, more rigorous approach, you are deriving $x$ by $\theta$ and $y$ by $r$, but you are forgetting the cross terms! These are function of two variables each, so you should do all partial derivatives:
$$
dx=\dfrac{\partial x}{\partial r}dr+\dfrac{\partial x}{\partial \theta}d\theta
$$
$$
dy=\dfrac{\partial y}{\partial r}dr+\dfrac{\partial y}{\partial \theta}d\theta\;.
$$
Deriving (you did already two of them, here we do all four):
$$
dx=\sin\theta\,dr+r\,\cos\theta\,d\theta
$$
$$
dy=\cos\theta\,dr-r\,\sin\theta\,d\theta\;.
$$
The exterior product of these two vectors (see link above, or think of the cross product, or think of the area of the parallelogram) is:
$$
dx \, dy = (\sin\theta\,dr)(-r\,\sin\theta\,d\theta)-(\cos\theta\,dr)(\cos\theta\,d\theta)
$$
$$
= -r(\sin^2\theta +\cos^2\theta)\,dr\,d\theta=-r\,dr\,d\theta.
$$
Usually one takes instead $x=\cos\theta,y=\sin\theta$, so that the minus sign becomes a plus. But anyway, what matters is the absolute value of that thing.
I think that what you were missing is the link between area (and more generally volumes) and exterior product. Look it up starting from the link above, it's very interesting.
A: The 'right-way' to do this is to use differential forms:
$$ dr \wedge d \theta = (\frac{\partial r}{\partial x} dx + \frac{\partial r}{\partial y} dy) \wedge ( \frac{\partial \theta}{\partial x} dx + \frac{\partial \theta}{\partial y} dy )= ( \frac{\partial r}{\partial x} \frac{\partial \theta}{\partial y} -  \frac{\partial r}{\partial y} \frac{\partial \theta}{\partial x} ) dx \wedge dy$$
Now just plug in $r$ and $\theta$ as functions of $x$ and $y$
A: How do you prove this rigorously?
The equality $dx\;dy= r\;dr\;d\theta$ by itself is nonsense (here, we don't have ordinary products in each side of the equality). What we should prove is
$$\iint_Qf(x,y)\;dx\;dy=\iint_A f(r\cos \theta,r\sin \theta)\;r\;dr\;d\theta\tag{$*$}.$$
I will not give a rigorous proof, but I will tell you how a formal proof can be done according to the approach of C. H. Edwards in his book Advanced calculus of several variables. Quoting this book:

The student has undoubtedly seen change of variables formulas such as $(*)$ which result from changes from rectangular coordinates to polar coordinates. The appearance of the factor $r$ in the formula is sometimes "explained" by mythical pictures, such as figure in the OP's post, in which it is alleged that $dA$ is an "infinitesimal" rectangle with sides $dr$ and $r\;d\theta$, and therefore has area $r \;dr\; d\theta$. In this section we shall give a mathemtaically acceptable explanation of the origin of such factors in the transformation of multiple integrals from one coordinate system to another.

To start, let us "remember" three things about double integrals:
$\quad$ (i) From the geometric interpretation of the double integral, if $X$ is a very small region and $(x_0,y_0)$ is some point in $X$, then
$$\iint_Xf(x,y)\;dx\;dy\cong f(x_0,y_0) \Delta X,$$
where $\Delta X$ is the area of $X$.
$\quad$ (ii) The double integral of a function $g$ over a rectangle $A\subset\mathbb{R}^2 $ can be written as a limit of Riemann sums:
$$\iint_Ag(x,y)\;dx\;dy=\lim_{\mathcal{\|\mathcal{P}\|}\to0}\sum_{i=1}^k g(x_i,y_i)\Delta A_i,$$
where $\mathcal{P}=\{A_1,...,A_k\}$ is a partition of $A$, $(x_i,y_i)$ is a point in $A_i$ and $\Delta A_i$ is the area of $A_i$. Thus, for a sufficiently fine partition,
$$\sum_{i=1}^k g(x_i,y_i)\Delta A_i\cong\iint_Ag(x,y)\;dx\;dy.$$
$\quad$ (iii) In the integral notation, $(x,y)$ is a dummy variable and thus
$$\iint g(x,y)\;dx\;dy=\iint g(r,\theta)\;dr\;d\theta.$$
Assume that we want to calculate the double integral
$$\iint_Q f(x,y)\;dx\;dy,\tag{0}$$
where $Q$ is an annular sector in $\mathbb{R}^2$ (see picture below).

We want to "transform" this integral into an integral over a rectangle. Why? Because we expect that the "transformed integral" will be easier to compute (by virtue of Fubini's theorem).
Note that the rectangle $A=\{(r,\theta)\in\mathbb{R}^2\mid a\leq r\leq b\text{ and }\alpha\leq \theta\leq\beta \}$ is transformed into the annular sector $Q$ by the mapping $T:\mathbb{R}^2\to\mathbb{R}^2$ given by $$T(r,\theta)=(r\cos\theta,r\sin\theta).\tag{1}$$

So,
$$Q=T(A).\tag{2}$$
We shall see that the mapping $T$ can be used to transform $(0)$ into an integral over the rectangle $A$, as we want. For this, we have to answer the

Big question: Given a rectangle $R$, what is the relationship between the area $\Delta R$ of $R$ and the area $\Delta T(R)$ of $T(R)$?

Answer: Look at the picture above. If $(c,d)$ is the center-point of $A$, then
$$\text{area of } Q=\frac{b^2(\beta-\alpha)}{2}-\frac{a^2(\beta-\alpha)}{2}=\frac{1}{2}(b+a)(b-a)(\beta-\alpha)=c(\text{area of } A)$$
The same is true for $R$ and thus
$$\Delta T(R)=p\Delta R,\tag{3}$$
where $p$ is the abscissa of the center-point of $R$.
Now, let $\mathcal{P}=\{A_1,..., A_k\}$ be a partition of $A$ into very small subintervals. For each $i=1,...,k$, let $(x_i,y_i)$ be the center-point of $A_i$. Then, $T(x_i,y_i)$ is a point in $T(A_i)$ and
$$T(A)=T(A_1)\cup\cdots \cup T(A_k)\tag{4}$$
with disjoint union. It follows that
\begin{align}
\iint_Qf(x,y)\;dx\;dy&=\iint_{T(A)} f(x,y)\;dx\;dy\tag{by (2)}\\\\
&=\sum_{i=1}^k\iint_{T(A_i)} f(x,y)\;dx\;dy\tag{by (4)}\\\\
&\cong\sum_{i=1}^k f(T(x_i,y_i))\Delta T(A_i)\tag{by (i)}\\\\
&=\sum_{i=1}^kf(T(x_i,y_i))x_i\Delta A_i\tag{by (3)}\\\\
&\cong \iint_A f(T(x,y))\;x\;dx\;dy\tag{by (ii)}\\\\
&=\iint_A f(x\cos y,x\sin y)\;x\;dx\;dy\tag{by (1)}\\\\
&=\iint_A f(r\cos \theta,r\sin \theta)\;r\;dr\;d\theta\tag{by (iii)}
\end{align}
By now, working with one sufficiently fine partition, we have deduced that
$$\iint_Qf(x,y)\;dx\;dy\cong\iint_A f(r\cos \theta,r\sin \theta)\;r\;dr\;d\theta.$$
However, with an "epsilon argument" it is possible to prove that we can make this approximation as good as we want and thus the equality $(*)$ is indeed true (the details are made in the book). Quoting the said book again, the proof of $(*)$

is a matter of supplying enough epsilonics to convert the above discussion into a precise argument.

A: In fact, what you did wrong is in confusing partial and total differentials (copy below) :

$$dx = \sin{\theta}dr+r\cos{\theta} d\theta$$
$$ dy = \cos\theta dr- r\sin{\theta} d\theta $$
One have to compute the determinant of the Jacobian matrix :
$dxdy=
\begin{Vmatrix}
  \frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} \\
  \frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} \\ 
 \end{Vmatrix}dr\;d\theta = 
\begin{Vmatrix}
  \sin(\theta) & r \cos(\theta) \\
  \cos(\theta) & -r \sin(\theta) \\ 
 \end{Vmatrix}dr\;d\theta =
r\left|-\sin^2(\theta)-\cos^2(\theta)\right|\,dr\,d\theta  $
$$dx\,dy= r\,dr\,d\theta$$ 
A: Along with other commenters, I was also confused by the accepted answer of @geodude. I'd never heard of wedge or exterior products before. But now I understand and can explain why @geodude's answer confused me.
When @geodude wrote "$dx\,dy = \ldots$", I thought I could just plug into $dx\,dy$ the formulas for $dx$ and $dy$ from the two lines previous. No! That's not right. $dx\,dy$ was meant to be the exterior product $dx\wedge dy$.
Really, you should be looking for a relationship between an area in the $(r,\theta)$ space (i.e., the exterior product $dr\wedge d\theta$) and an area in the $(x,y)$ space (i.e., $dx\wedge dy$).
I was tempted and thought $dx\,dy$ was an area in the $(x,y)$ space because it looked like the area formula for a rectangle in the $(x,y)$ space. I was wrong. $dx\wedge dy$ is the area.
I'll write out the answer, combining that of @geodude and @Tryst.with.Freedom.
Let $T: \mathbb{R}^2 \rightarrow \mathbb{R}^2,$
$$\mathbf{T}\left(\begin{bmatrix}x \\ y \end{bmatrix}\right)=\begin{bmatrix}r\cos\theta \\ r\sin\theta \end{bmatrix}$$
then the total derivative of $\mathbf{T}$ is
$$ \mathbf{T}'\left(\begin{bmatrix}x \\ y \end{bmatrix}\right)=\begin{bmatrix}\cos\theta & -r\sin\theta\\ \sin\theta & r\cos\theta\end{bmatrix}.$$
In terms of differentials, that's
$$\begin{bmatrix}dx \\ dy \end{bmatrix}=\begin{bmatrix}\cos\theta & -r\sin\theta\\ \sin\theta & r\cos\theta\end{bmatrix}\begin{bmatrix}dr \\ d\theta \end{bmatrix}$$
or
$$\begin{align}dx &= \cos\theta dr -r\sin\theta d\theta \\ dy &= \sin\theta dr + r\cos\theta d\theta \end{align}.$$
BAD: This where I multiplied $dx$ and $dy$
$$\begin{align}
dx\,dy &= (\cos\theta dr -r\sin\theta d\theta)(sin\theta dr + r\cos\theta d\theta) \\
&= \cos\theta\sin\theta(dr)^2+r\cos^2\theta dr d\theta - r\sin^2\theta dr d\theta - r^2\cos\theta\sin\theta(d\theta)^2 \\
&=r\cos^2\theta dr d\theta - r\sin^2\theta dr d\theta\\
&=rdrd\theta(\cos^2\theta - \sin^2\theta) \\
&\neq rdrd\theta
\end{align}$$
I looked at all my trig identities. No dice. I couldn't get to $dx\,dy = rdrd\theta$. That's because it was never going to work!
GOOD: What I needed to do was the following:
$$\begin{align}
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^2\cos\theta\sin\theta)d\theta \wedge d\theta \\
&= (r\cos^2\theta) dr \wedge d\theta - (r\sin^2\theta) d\theta \wedge dr \\
&= (r\cos^2\theta) dr \wedge d\theta + (r\sin^2\theta) dr \wedge d\theta \\
&=  r(\cos^2\theta + \sin^2\theta)dr \wedge d\theta \\
&= r\,dr \wedge d\theta
\end{align}$$
I belive this allows you to heuristically interpret
$$ dxdy=rdrdθ $$
