¿Cómo se demuestra esto rigurosamente?
La igualdad $dx\;dy= r\;dr\;d\theta$ por sí sola no tiene sentido (aquí, no tenemos productos ordinarios en cada lado de la igualdad). Lo que debemos probar es $$\iint_Qf(x,y)\;dx\;dy=\iint_A f(r\cos \theta,r\sin \theta)\;r\;dr\;d\theta\tag{$*$}.$$
No daré una prueba rigurosa, pero te diré cómo se puede hacer una prueba formal según el enfoque de C. H. Edwards en su libro Cálculo avanzado de varias variables. Citando este libro:
Sin duda, el estudiante ha visto fórmulas de cambio de variables como $(*)$ que resultan de cambios de coordenadas rectangulares a coordenadas polares. La aparición del factor $r$ en la fórmula a veces se "explica" con imágenes míticas, como la figura en el post original, en la que se afirma que $dA$ es un rectángulo "infinitesimal" con lados $dr$ y $r\;d\theta$, y por lo tanto tiene un área de $r \;dr\; d\theta$. En esta sección daremos una explicación matemáticamente aceptable del origen de tales factores en la transformación de integrales múltiples de un sistema de coordenadas a otro.
Para empezar, "recordemos" tres cosas sobre integrales dobles:
$\quad$ (i) Desde la interpretación geométrica de la integral doble, si $X$ es una región muy pequeña y $(x_0,y_0)$ es algún punto en $X$, entonces $$\iint_Xf(x,y)\;dx\;dy\cong f(x_0,y_0) \Delta X,$$ donde $\Delta X$ es el área de $X$.
$\quad$ (ii) La integral doble de una función $g$ sobre un rectángulo $A\subset\mathbb{R}^2 $ puede escribirse como límite de sumas de Riemann: $$\iint_Ag(x,y)\;dx\;dy=\lim_{\mathcal{\|\mathcal{P}\|}\to0}\sum_{i=1}^k g(x_i,y_i)\Delta A_i,$$ donde $\mathcal{P}=\{A_1,...,A_k\}$ es una partición de $A$, $(x_i,y_i)$ es un punto en $A_i$ y $\Delta A_i$ es el área de $A_i$. Así, para una partición suficientemente fina, $$\sum_{i=1}^k g(x_i,y_i)\Delta A_i\cong\iint_Ag(x,y)\;dx\;dy.$$
$\quad$ (iii) En la notación de la integral, $(x,y)$ es una variable dummy y por lo tanto $$\iint g(x,y)\;dx\;dy=\iint g(r,\theta)\;dr\;d\theta.$$
Supongamos que queremos calcular la integral doble $$\iint_Q f(x,y)\;dx\;dy,\tag{0}$$ donde $Q$ es un sector anular en $\mathbb{R}^2$ (ver la imagen abajo).
Queremos "transformar" esta integral en una integral sobre un rectángulo. ¿Por qué? Porque esperamos que la "integral transformada" sea más fácil de calcular (por virtud del teorema de Fubini).
Observa que el rectángulo $A=\{(r,\theta)\in\mathbb{R}^2\mid a\leq r\leq b\text{ y }\alpha\leq \theta\leq\beta \}$ se transforma en el sector anular $Q$ por la función $T:\mathbb{R}^2\to\mathbb{R}^2$ dada por $$T(r,\theta)=(r\cos\theta,r\sin\theta).\tag{1}$$
Entonces, $$Q=T(A).\tag{2}$$
Veremos que la función $T$ se puede usar para transformar $(0)$ en una integral sobre el rectángulo $A$, como queremos. Para esto, tenemos que responder la
Gran pregunta: Dado un rectángulo $R$, ¿cuál es la relación entre el área $\Delta R$ de $R$ y el área $\Delta T(R)$ de $T(R)$?
Respuesta: Mira la imagen anterior. Si $(c,d)$ es el punto central de $A$, entonces $$\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)$$ Lo mismo es cierto para $R$ y por lo tanto $$\Delta T(R)=p\Delta R,\tag{3}$$ donde $p$ es la abscisa del punto central de $R$.
Ahora, sea $\mathcal{P}=\{A_1,..., A_k\}$ una partición de $A$ en subintervalos muy pequeños. Para cada $i=1,...,k$, sea $(x_i,y_i)$ el punto central de $A_i$. Entonces, $T(x_i,y_i)$ es un punto en $T(A_i)$ y $$T(A)=T(A_1)\cup\cdots \cup T(A_k)\tag{4}$$ con unión disjunta. Se sigue que \begin{align} \iint_Qf(x,y)\;dx\;dy&=\iint_{T(A)} f(x,y)\;dx\;dy\tag{por (2)}\\\\ &=\sum_{i=1}^k\iint_{T(A_i)} f(x,y)\;dx\;dy\tag{por (4)}\\\\ &\cong\sum_{i=1}^k f(T(x_i,y_i))\Delta T(A_i)\tag{por (i)}\\\\ &=\sum_{i=1}^kf(T(x_i,y_i))x_i\Delta A_i\tag{por (3)}\\\\ &\cong \iint_A f(T(x,y))\;x\;dx\;dy\tag{por (ii)}\\\\ &=\iint_A f(x\cos y,x\sin y)\;x\;dx\;dy\tag{por (1)}\\\\ &=\iint_A f(r\cos \theta,r\sin \theta)\;r\;dr\;d\theta\tag{por (iii)} \end{align}
En este punto, trabajando con una partición suficientemente fina, hemos deducido que $$\iint_Qf(x,y)\;dx\;dy\cong\iint_A f(r\cos \theta,r\sin \theta)\;r\;dr\;d\theta.$$ Sin embargo, con un "argumento epsilon" es posible probar que podemos hacer esta aproximación tan buena como queramos y así la igualdad $(*)$ es verdadera (los detalles se encuentran en el libro). Citando nuevamente el libro mencionado, la prueba de $(*)$
es cuestión de suministrar suficientes argumentos epsilon para convertir la discusión anterior en un argumento preciso.