The surface integral over surface $S$ (which is given by $z=f(x,y)$, where $(x,y)$ is a point from the region $D$ in the $xy$-plane) is:
$$ \iint\limits_{S} g(x,y,z)\ dS = \iint\limits_{D} g(x,y,f(x,y))\ {{\sqrt {\,{{\left[ {{f_x}} \right]}^2} + {{\left[ {{f_y}} \right]}^2} + 1} \,dA}}$$
Does there exist a rigorous proof of this formula in real analysis.
On
There are various intuitive arguments showing that your formula, or the standard formula for the computation of surface area, is correct, and gives the expected results. But a full proof necessitates exact definitions of "surface" and "area", and the invariance of area under orthogonal mappings, etcetera.
If ${\rm d}S$ is a tiny piece of your surface $S$, located at the point $p\in S$, and this piece is projected orthogonally to the $(x,y)$-plane, then we obtain there a tiny "area element" ${\rm d}A$. The area scaling factor between these two pieces is $\cos\phi$, whereby $\phi$ is the angle between the tangent plane of $S$ at $p$ and the $(x,y)$-plane: $${\rm d}S={1\over\cos\phi}\,{\rm d}A\ .$$ It is easy to verify that $${1\over\cos\phi}=\sqrt{1+f_x^2+f_y^2}$$ when $S$ is given in the form $z=f(x,y)$.
On
Let $S:\mathbb{R}^2\rightarrow \mathbb{R}^3$ be a surface mapping. This is a more general form for a surface mapping which includes the specific type $S(u,v) = (u, v, f(u,v))$. What is the area of the surface when the $(u,v)$ domain is the rectangle $\mathcal{R}$ with corners $(0,0)$ and $(a,b)$? We can approximate the area as a sum of areas of parallelograms, each with three vertices on $S$. Create an $N$-by-$N$ rectangular grid on $\mathcal{R}$ (with each subrectangle having a width of $\Delta u = ah$ and height of $\Delta v = bh$, where $h=1/N$). Parallelogram $(i,j)$ has corners $P_0=S(u_i,v_j)$, $P_1=S(u_i+\Delta u, v_j)$, $P_2=S(u_i,v_j+\Delta v) + S(u_i+\Delta u,v_j) - S(u_i,v_j)$, and $P_3=S(u_i,v_j+\Delta v)$. The area of the parallelogram is $A_{ij} = |(P_1-P_0) \times (P_3-P_0)|$. An alternate definition of the partial derivative $S_u$ at the point $(u,v)$ is $$ S_u(u,v)h = S(u+h,v) - S(u,v) - R(h)h $$ where $R(h)$ is some continuous vector remainder function with $R(h)=0$. Using this definition, $$ \begin{align} A_{ij} &= \left|(S_u(u_i,v_j)\Delta u + R_1(h)h) \times (S_v(u_i,v_j)\Delta v + R_2(h)h)\right| \\ &= \left| (S_u \times S_v)\Delta u \Delta v + R_{ij}(h) h^2\right|. \end{align} $$ Using the triangle inequality, $$ |S_u \times S_v| \Delta u \Delta v - |R_{ij}(h)|h^2 \le A_{ij} \le |S_u\times S_v|\Delta u \Delta v + |R_{ij}(h)|h^2. $$ Using the squeeze theorem, we find that $$ \begin{align} \lim_{N\rightarrow\infty} \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} A_{ij} &= \lim_{N\rightarrow\infty} \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} |S_u(u_i,v_j)\times S_v(u_i,v_j)|\Delta u \Delta v \\ &= \iint_{\mathcal{R}} |S_u(u,v) \times S_v(u,v)| du dv. \end{align} $$ The sum of $|R_{ij}(h)|h^2$ is smaller than or equal to $R_{\textrm{max}}(h)$ where $R_{\textrm{max}}(x)$ is the maximum of all $|R_{ij}(x)|$. Since the limit of $R_{\textrm{max}}(x)$ as $x\rightarrow 0$ is zero, the limit of the sum of $|R_{ij}(h)|h^2$ is zero. This result can be extended to a general subset $\mathcal{R}$ of $\mathbb{R}^2$ since a general domain for integration in $\mathbb{R}^2$ can be approximated arbitrarily closely with a set of disjoint rectangles.
To strengthen the case that $$ \iint_{\mathcal{R}} |S_u(u,v) \times S_v(u,v)| du dv \tag{1}\label{surface_area_integral} $$ is the actual surface area, consider what happens when we move the fourth point of the parallelogram onto the surface $S$. This means replacing point $P_2$ with $P_2' = S(u_i+\Delta u, v_j+\Delta v)$. Now instead of a parallelogram, we have two triangles $P_0 P_1 P_3$ and $P_1 P_2' P_3$. The difference between the area of the parallelogram and the sum of the areas of these two triangles is the difference between the area of triangle $P_1 P_2 P_3$ (denote this area by $T$) and triangle $P_1 P_2' P_3$ (denote this area by $T'$). First consider the vector from $P_2$ to $P_2'$: $$ \begin{align} P_2' - P_2 &= (S(u_i+\Delta u,v_j+\Delta v) - S(u_i,v_j+\Delta v)) - (S(u_i+\Delta u,v_j) - S(u_i,v_j)) \nonumber \\ &= S_u(u_i,v_j+\Delta v)\Delta u - S_u(u_i,v_j)\Delta u + R_1(h)h \nonumber \\ &= (S_u(u_i,v_j+\Delta v) - S_u(u_i,v_j))\Delta u + R_1(h)h \nonumber \\ &= R(h)h \tag{2}\label{P2p_minus_P2} \end{align} $$ where $R_1(h)$ and $R(h)$ are some continuous vector functions with $R_1(0)=R(0)=0$. The last equality in \eqref{P2p_minus_P2} is true because we are assuming $S_u$ is a continuous function, meaning $S_u(u_i,v_j+\Delta v) - S_u(u_i,v_j)$ is a continuous function of $h$ equal to 0 at $h=0$. By the triangle inequality, $P_1 P_2$ differs from $P_1 P_2'$ by at most $|R(h)|h$. Similarly, $P_3 P_2$ differs from $P_3 P_2'$ by at most $|R(h)|h$. We are also assuming that $S_u \ne 0$ and $S_v \ne 0$ (except possibly at a finite number of points). So $c=P_1P_2 = \tilde{R}_1(h)h$ and $d=P_3P_2 = \tilde{R}_2(h)$ for continuous functions $\tilde{R}_1$ and $\tilde{R}_2$ with $\tilde{R}_1(0) \ne 0$ and $\tilde{R}_2(0)\ne 0$. Using the formula $2T=cd \sin\gamma$ (where $\gamma$ is the angle between edge $P_1P_2$ and edge $P_2 P_3$), we can say that $T=\tilde{R}_3(h) h^2$ (where $\tilde{R}_3$ is continuous with $\tilde{R}_3(0) \ne 0$) if we also assume that $\gamma$ is always greater than 0 and less than $\pi$ (which will be true if the angle between $S_u$ and $S_v$ is always greater than 0 and less than $\pi$). Using these facts about $c$, $d$, and $T$ along with Heron's formula for the area of a triangle shows that $T'^2 - T^2 = R_2(h)h^4$ for some continuous function $R_2$ with $R_2(0)=0$. Since $T' - T = (T'^2 - T^2) / (T' + T)$, $T'-T = R_3(h)h^2$ for some continuous function $R_3$ with $R_3(0)=0$. So moving $P_2$ onto the surface simply changes $A_{ij}$ by $R_3(h)h^2$. So the limit of the sum of areas of the $2N^2$ watertight triangles with all vertices on the surface is \eqref{surface_area_integral}.
On
I think an old school proof, and how I proved, is the following:
Imaging you make a grid in your surface, where from the above or below point of view, it seems like a normal homogeneous square grid, but when you see it for another point of view, it seems deformed.
That deformation is because you are proyecting that grid to a surface which is not parallel everywhere to that grid; you could try to graph any surface in geogebra to see what i mean.
That deformed grid has an horizontal component parallel to the $x$-axis and an horizontal component parallel to the $y$-axis. The horizontal components can be draw as vectors parallel to each axis. Now, What are thoose vector? Imaging you choose a parallelogram that aproximate to a component of the grid. To get the area $\Delta A$ of that parallelogram, we can use the fact that the norm of the cross product of two vectors is equal to the area of the parallelogram formed by these to vectors; so,the vectors need to be tangent to that surface, because they aproximate to the parallelogram of the grid, named: $\vec {r_x}$ and $\vec {r_y}$. The formula that serves to calculate the aproximated area of each parallelogram is:
$$\Delta A=\left\Vert\vec{r_x}\times\vec{r_y}\right\Vert$$
To get $\vec {r_x}$, think that the coordinates point in the x-axis, because is parallel to that, but also point in the z-axis, such that the vector is tangent to the surface. If take the difference of to vectors $p$ and $q$, where $p$ is the point where we are compute $\vec {r_x}$, so:
$$q - p = \vec {r_x}$$ Define $q$ as: $$q=(x+\Delta x,y,f(x+\Delta x,y))$$ and $p$ as: $$p=(x,y,f(x,y))$$ Then taking: $$\lim_{\Delta x\to 0}{q-p}=\lim_{\Delta x\to 0}{(x+\Delta x,y,f(x+\Delta x,y))-(x,y,f(x,y))}$$
$$=\lim_{\Delta x\to 0}{(\Delta x,0,f(x+\Delta x, y)-f(x,y))}$$ Sadly, that limit converges to $0$. To make it useful, we need to multiply by ${\Delta x}\over {\Delta x}$ as follow:
$$=\lim_{\Delta x\to 0}{(\Delta x,0,f(x+\Delta x, y)-f(x,y)){\Delta x}\over {\Delta x}}$$
$$=\lim_{\Delta x\to 0}{\left(1,0,{f(x+\Delta x, y)-f(x,y)\over \Delta x}\right)}\Delta x$$
The term ${f(x+\Delta x, y)-f(x,y)\over \Delta x}$ is the definition of the partial derivative of $f$ with respect to $x$, si the above limit result as:
$$\lim_{\Delta x\to 0}{\vec {r_x}}=\left(1,0,{\partial f\over\partial x}\right)dx$$
The same argument but for $\vec {r_y}$ can be done. We now have:
$$\vec{r_x}=\left(1,0,{\partial f\over\partial x}\right)dx;\ \vec{r_y}=\left(0,1,{\partial f\over\partial y}\right)dy$$
where the differentials $dx$ and $dy$ comes from multiplying by $\Delta x\over\Delta x$ and $\Delta y\over\Delta y$ as showed above and taking the limit: $\lim_{(\Delta x,\Delta y)\to(0,0)}(\Delta x,\Delta y)=(dx,dy)$.
Using the formula for $\Delta A$, after taking the limit explained, yields:
$$dA=\left\Vert\vec{r_x}\times\vec{r_y}\right\Vert=\left\Vert\left(1,0,{\partial f\over\partial x}\right)dx\times\left(0,1,{\partial f\over\partial y}\right)dy\right\Vert $$
$$=\left\Vert\left(1,0,{\partial f\over\partial x}\right)\times\left(0,1,{\partial f\over\partial y}\right)\right\Vert dx\ dy$$
Note that $\Delta A$ converts to $dA$ because "deltas" converge to "differentials" when "deltas" are function of some parameter $\alpha$ and $\Delta\alpha\to 0$.
Computing that cross product and taking the norm gives:
$$=\left\Vert\left(1,0,{\partial f\over\partial x}\right)\times\left(0,1,{\partial f\over\partial y}\right)\right\Vert dx\ dy = \sqrt{\left(\partial f\over\partial x\right)^2+\left(\partial f\over\partial y\right)^2+1}\ dx\ dy $$
$$dA= \sqrt{\left(\partial f\over\partial x\right)^2+\left(\partial f\over\partial y\right)^2+1}\ dx\ dy$$
For last, How do we know that that formula is somewhat correct? Because the norm is a postive value, the same for the area of a surface. Also, if you parametrize the surface by a vector $\vec r$ and two parameters $u$ and $v$:
$$\vec r(u,v) = (x(u,v),y(u,v),z(u,v)) $$
$$dA=\left\Vert\vec{r_u}\times\vec{r_v}\right\Vert=\left\Vert\left({\partial x\over\partial u},{\partial y\over\partial u},{\partial z\over\partial u}\right)\times\left({\partial x\over\partial v},{\partial y\over\partial v},{\partial z\over\partial v}\right)\right\Vert du\ dv$$
and the formula for what you asked is just the case of parametrization by $x$ and $y$.
In general you have to transform the "volume form" of the three dimensional Euclidean space into the volume form of the embedded surface: there is a nice formula involving the determinant of the metric tensor that allows you to do that. However, to understand the general procedure you need some differential geometry.
What follows is not "rigorous" but gives you the idea, and depending on your definition of what is "rigorous" it may constitute a "visual proof".
First, given two 3-vectors $a$ and $b$, the area of the parallelogram of sides $a$ and $b$ is $|a \times b|$.
Now consider a unit square of sides $(1,0)$ and $(0,1)$ in the plane $(x,y)$, namely in the domain of the function $f$. This square is centered somewhere (say in $(x_0,y_0)$), not necessarily at the origin.
Imagine to project this square on the plane tangent to $f$ at the point $(x_0,y_0,f(x_0,y_0))$: the sides will be $(1,0,f_x(x_0,y_0))$ and $(0,1,f_y(x_0,y_0))$. Why? Because when you move in a direction you must move a bit also along the $z$-direction to remain on the tangent plane (its inclination is described by the gradient of $f$).
Therefore, the "unit square" in the domain draws an area $$|(1,0,f_x )\times (0,1,f_y )|=\sqrt{f_x^2 + f_y^2 +1}$$ on the tangent plane. This "unit area" on the tangent plane is used locally to rescale the natural measure $dx dy$ on the domain of $f$.