Origin of Jacobian determinant What is the origin of the Jacobian determinant for changing variables in multiple integrals?
I mean, how to derive the formula for the Jacobian determinant?
I have seen the use of Jacobian in some books but could not find how it arises. 
 A: This is an intuitive explanation, not a rigorous explanation that you would see in an analysis class.  
Let's start in $\mathbb{R}.$ Take the interval $[a,b] \subset \mathbb{R},$ where $a \leq b.$ We measure volume of the interval by taking length $b-a$. Now let's say we have an invertible linear transformation $T: [a,b] \rightarrow \mathbb{R}$ defined $$T(x)=Ax,$$ in which $A$ is a $1 \times 1$ matrix with a nonzero entry $m.$ How do the volumes of the preimage $[a,b]$ and the image $T([a,b])$ compare to one another ? If $m>0,$ then we see that $$\text{Vol}(T([a,b]))=\text{Vol}([ma,mb])=m \text{Vol}([a,b]).$$ On the other hand, if $m<0,$ then
$$\text{Vol}(T([a,b]))=\text{Vol}([mb,ma])= -m \text{Vol}([a,b]).$$ In both cases, we see that $$\text{Vol}(T([a,b]))=|m| \text{Vol}([a,b])=|\det(A)|\text{Vol}([a,b]).$$ Now let's talking about change of variables in one dimension, which is commonly known as u-substitution. Consider 
$$\int_{I} f(x) \ dx,$$ in which $I$ is some interval. Let's make a substitution $$x=g(u).$$ Then our integral becomes $$\int_{J} f(g(u)) \left|\frac{dg}{du}\right| \ du,$$  where $J$ is the image of $I$ under the transformation $g.$ Some notes about using this formula. $J$ must be written in the form $[c,d]$ in which $c \leq d.$ Also intuitively, $\ dx$ is the length of an infinitesimal interval in the $x$ axis and that $\ du$ is the length of an infinitesimal interval in the $u$ axis, and $ \left|\frac{dg}{du}\right|$ plays the role of the scaling factor like $|\det(A)|$ does.  
Now let's talk about $\mathbb{R}^2.$ We take the rectangle $[a_1,b_1] \times [a_2,b_2] \subset \mathbb{R}^2.$ Here each $a_i \leq b_i.$ We measure volume through area, which is $(b_1-a_1)(b_2-a_2).$ Let's say we have an invertible linear transformation $T:[a_1,b_1] \times [a_2,b_2] \rightarrow \mathbb{R}^2,$ and $$T(x)=Ax$$ in which $A$ is a $2 \times 2$ matrix with a nonzero determinant.  If we compare the areas of the preimage $[a_1,b_1] \times [a_2,b_2]$ and the image $T([a_1,b_1] \times [a_2,b_2]),$ we will find that $$\text{Vol}(T([a_1,b_1] \times [a_2,b_2]))= |\det(A)|\text{Vol}([a_1,b_1] \times [a_2,b_2]).$$ Considering the 2 dimensional change of variables formula, if we have 
$$\int_{R} f(x,y) \ dy \ dx,$$ with $R$ being a rectangle, and make change of variables $x=g(u,v),y=h(u,v),$ then our integral becomes 
$$\int_{S} f(x(u,v),y(u,v)) \left| \det \left(  \frac{\partial(x,y)}{\partial(u,v)} \right ) \right| \ du \ dv,$$ where $S$ is the image of $R$ under the given transformation. As before, $S$ must be written in the form $[c_1,d_1] \times [c_2,d_2]$ in which each $c_i \leq d_i.$ Secondly $dy \ dx$ is the area of an infinitesimal rectangle in the $x-y$ plane, $du \ dv$ is the area of an infinitesimal rectangle in the $u-v$ plane. The determinant of the Jacobian $\left| \det \left(  \frac{\partial(x,y)}{\partial(u,v)} \right ) \right|$ acts as a scaling factor like $|\det(A)|$ does.  
The higher dimensional version of the rectangle and interval is the $n$ dimensional rectangle defined by $[a_1,b_1] \ times ... \times [a_n,b_n]$ which has hypervolume $\prod_{i=1}^{n} (b_i-a_i)$ If we consider an invertible transformation $T: [a_1,b_1] \times ... \times [a_n,b_n] \rightarrow \mathbb{R}^n,$ by which $T(x)=Ax$ and $A$ is some $n \times n$ matrix that has nonzero determinant, we have the general formula 
$$\text{Vol}(T([a_1,b_1] \times ... \times [a_n,b_n])) = |\det(A)| \text{Vol}([a_1,b_1] \times ... \times [a_n,b_n]).$$ Now consider the general change of variables formula. If we have 
$$\int_{R} f(x_1,...,x_n) \ dx_n ... \ dx_1,$$ with $R$ being an $n$ dimensional rectangle, and make change of variables $x_1=g_1(u_1,...,u_n),...,x_n=g_1(u_1,...,u_n)$ then our integral becomes 
$$\int_{S} f(x_1(u_1,...,u_n),...,x_n(u_1,...,u_n)) \left| \det \left(  \frac{\partial(x_1,...,x_n)}{\partial(u_1,...,u_n)} \right ) \right| \ du_n ... \ du_1.$$ As before, $S$ must be written in the form $[c_1,d_1] \times ... \times [c_n,d_n]$ in which each $c_i \leq d_i.$ Secondly $dx_n \ ... \ dx_1$ is the area of an infinitesimal rectangle in the $x_1...x_n$ hyperplane, $du_n \ ... \ du_1$ is the area of an infinitesimal rectangle in the $u_1...u_n$ plane. The determinant of the Jacobian $\left| \det \left(  \frac{\partial(x_1,...,x_n)}{\partial(u_1,...,u_n)} \right ) \right|$ acts as a scaling factor like $|\det(A)|$ does. 
A: I will try to give an answer that you can apply with little knowledge beyond what is taught in a standard multi-variable course.

What is the origin of the Jacobian determinant for changing variables in multiple integrals?

When dealing with multi-variable calculus, one often finds expressions such as
$$
\iiint_V f(x,y,z)\, dx\,dy\,dz \;.
$$
The first thing to understand is that $dx\,dy\,dz$ is not an "ordinary" product of infinitesimals. If you try to compute the quantity
$$
\iint_D \sqrt{x^2+y^2} \, dx\, dy \;,
$$
where $D=\{x,y \in \mathbb{R} : 1<\sqrt{x^2+y^2}<2 \}$, you can easily note that this integral is much simpler in polar coordinates. Going step by step, we have that $x=r\cos(\theta)$ and $y=r\sin(\theta)$ and therefore you will have the infinitesimals, $dx=\cos(\theta)dr - r\sin(\theta)d\theta$ and $dy=\sin(\theta)dr + r\cos(\theta) d\theta$.
If you change the variables naively you will get that
$$
dx\,dy= (\cos(\theta)dr - r\sin(\theta)d\theta)(\sin(\theta)dr + r\cos(\theta) d\theta) \;
$$
which is certainly not right because we are not preserving the area of the region $D$, as already explained in the other answers.
It turns out that the rigorous definition of the quantity $dx$ explains why the product $dx\,dy$ is actually omitting that this is an anticommutative product (meaning that $u*v=-v*u$) called Wedge product and in fact we say that we are integrating $f(x,y) \, dx \wedge dy = - f(x,y) \, dy \wedge dx$.
Going back to our example, we have that
$$
dx\,dy = (\cos(\theta)dr - r\sin(\theta)d\theta)\wedge(\sin(\theta)dr + r\cos(\theta) d\theta) \; ,
$$
and therefore, because $du\wedge du= - du\wedge du \Rightarrow du\wedge du=0 $, we have that
$$
dx\,dy = r \left(\cos^2(\theta) + \sin^2(\theta)\right) = r \, dr\, d\theta \; ,
$$
which is the right transformation for the variable. If one wishes to generalize this result, it is simple that from 
$$
dx=\frac{\partial x(r,\theta)}{\partial r} dr + \frac{\partial x(r,\theta)}{\partial  \theta} d\theta \; ,
$$
and the anticommutativity of the Wedge product, you will always get a determinant of the Jacobian matrix. The great advantage of realizing that the Jacobian arises from defining volumes from this products is that you can actually apply this not only to $\mathbb{R}^n$ but also to many other kinds of "well behaved" manifolds.

I mean, how to derive the formula for the Jacobian determinant?

Using the generalization in $\mathbb{R}^n$, although not necessarily needed as stated above, we have that
$$
dx_1\wedge dx_2 \wedge\dotsb\wedge dx_n \mathrm \quad \text{where} \quad x_k=x_k(y_1,...,y_n) \; ,
$$
and because we know that
$$ dx_k = \sum_{i=1}^n \frac{\partial x_k}{\partial y_i } dy_i \; , $$
from the definition of the Wedge product we have that
$$
dy_1\wedge dy_2 \wedge\dotsb\wedge dy_n = 
\det(J) dx_1\wedge dx_2 \wedge\dotsb\wedge dx_n \; ,
$$
where $J$ is the Jacobian matrix. Showing that the Wedge product is closely related to the definition of the Jacobian determinant goes as
\begin{align*}
dx_1\wedge dx_2 \wedge\dotsb\wedge dx_n &= \sum_{i_1=1}^n \frac{\partial x_1}{\partial y_{i_1} } dy_{i_1} \wedge \dotsb \wedge \sum_{i_n=1}^n \frac{\partial x_n}{\partial y_{i_n} } dy_{i_n}\\
&=\sum_{i_1,...,i_n} \epsilon^{i_1 ... i_n} J_{1,i_1}...J_{n,i_n} dy_1\wedge dy_2 \wedge\dotsb\wedge dy_n\\
&= \det(J) dy_1\wedge dy_2 \wedge\dotsb\wedge dy_n \; ,
\end{align*}
where $\epsilon^{i_1 ... i_n}$ is the Levi-Civita symbol.
One might argue that because of the anticommutativity of the Wedge product, we may have a minus sign on $\det(J)$ depending on how we order the product. This issue is solved by considering the absolute value of the Jacobian determinant when computing integrals in multi-variable calculus. Without aking the absolute value, we will have to understand the meaning of a signed area which is an interesting but not important topic for solving these integrals.
A: First of all, if $\Omega$ is a subset of $\mathbb R^n$ and $A \in \mathbb R^{n \times n}$ is a matrix, the measure of $A \Omega $ is equal to $|\det A|$ times the measure of $\Omega$. (One way to understand this fact clearly is to use the SVD of $A$. "Measure" just means "volume" if $n=3$ and "area" if $n=2$.)
Next, a differentiable function $f:\mathbb R^n \to \mathbb R^n$ can be approximated near a point $x_0$ by the function $x \mapsto f(x_0) + f'(x_0) (x-x_0)$. So, if $\Omega$ is a very small set containing $x_0$, the measure of the image of $\Omega$ is approximately $|\det f'(x_0) |$ times the measure of $\Omega$. We have discovered the Jacobian determinant and we see its significance.
(A comment about how to make this intuition rigorous: If you're worried about the precise definition of "measure", you could look at how "Lebesgue measure" is defined in real analysis textbooks. Technically I should assume the set $\Omega$ is "measurable" before referring to its "measure".)
