why divergence only cares about partial derivative of x, y, z, respectively? suppose that vector field F = [3xz, 2xy, -yz^2]
then divergence of F is 3z + 2x - 2yz
but I wonder why doesn't consider partial derivative of 3xz w.r.t z, 
partial derivative of 2xy w.r.t x, and partial derivative of -yz^2 w.r.t y
For example, when I change x, 3xz and 2xy changes together, I think, in divergence formula, it should contain the variation of 2xy when change x.
When vector field F(v1, v2, v3) is given, I think divergence formula should be like this.
$$
divF = \frac{\partial v_1}{\partial x} + \frac{\partial v_1}{\partial y} + \frac{\partial v_1}{\partial z} + \frac{\partial v_2}{\partial x} + \frac{\partial v_2}{\partial y} + \frac{\partial v_2}{\partial z} + \frac{\partial v_3}{\partial x} + \frac{\partial v_3}{\partial y} + \frac{\partial v_3}{\partial z}
$$
help me please
 A: To me, divergence at a point signifies the flux of a field out of a small volume (for simplicity, a cube) centered at that point. More specifically, if have a vector field $F$, and we define $C(r,x)$ to be the cube with sides of length $r$, centered at $x$ and with sides aligned with the coordinate axes, then the divergence of $F$ at $x$ is
$$
\nabla F(x)=\lim_{r\to 0}\frac{1}{r^3}\iint_{\partial C(r,x)}F\cdot \vec n\, dA
$$
I'm not going to do all the details of the derivation here, but we can see what happens for a small cube: Take the two faces that are normal to the $x$-axis, and see what happens to the contribution for the integral there. The normal vectors are $[-1,0,0]$ and $[1,0,0]$ respectively, and that makes the contribution equal to
$$
\frac1{r^3}\left(\iint_{\text{Face }1} F\cdot [-1,0,0]\,dA+\iint_{\text{Face }2} F\cdot [1,0,0]dA\right)\\
=\frac1{r^3}\left(\iint_{\text{Face }1} -F_x\,dA+\iint_{\text{Face }2} F_x\,dA\right)
$$
where $F_x$ is the $x$-component of $F$. What are these integrals? They are each done over a face with area $r^2$, and if the cube is small, that means each of them is equal to the area times the value of the integrand at the center point of the face. We get, approximately
$$
=\frac1{r^3}\left(-r^2F_x(x-[r/2,0,0])+r^2F_x(x+[r/2,0,0])\right)\\
=\frac{F_x(x+[r/2,0,0])-F_x(x-[r/2,0,0])}{r}
$$
If we let $r\to 0$ here, we recognize this as $\frac{\partial F_x}{\partial x}$. Similarly, the contribution from the other two pairs of faces become $\frac{\partial F_y}{\partial y}$ and $\frac{\partial F_z}{\partial z}$. Sum them up, and you have the entire divergence.
A: That is a very reasonable guess for what $\textrm{div}F$ should be, if you're thinking of it as a measure of ''how much $F$ is changing'' in some vague sense. 
However, there is a more precise quantity of ''change" which we actually want the divergence to measure: at a given point, the divergence measures the flux of $F$ through a small sphere centered at that point.
In other words, suppose at each point in space there is a (flying) ant, and the vector field $F$ describes the velocity the ant's motion at each point. Then we can ask, after 1 second passes will there be more or less (fewer?) ants in a given region?
The answer is that there will be less ants if $\textrm{div}F > 0$, and more ants if $\textrm{div}F<0$. [see https://ka-perseus-images.s3.amazonaws.com/733847896ca352cc429608496cfd1a46a025ad30.svg]
Now, how does this idea translate into the standard formula for divergence? Well, we can test it on two simpler examples: consider
$$ F_1 = (x,0,0) \quad\text{and}\quad F_2 = (y,0,0). $$
If we think of these two vector fields as describing configurations of ant movements, in each case do the ants become more or less concentrated at a given point? [It may help to sketch out the vector fields, which I haven't done out of laziness.] 
The standard divergence formula says $\textrm{div}F_1 = \frac{\partial x}{\partial x} = 1 > 0$ while $\textrm{div}F_2 = \frac{\partial y}{\partial x} =  0$. This should suggest that the ants moving according to $F_1$ will become more sparse, while the ants moving according to $F_2$ will not change their spatial density over time.
Does this match what you expect from looking at the two vector fields?
(Beyond looking at these two simple examples, I think the divergence formula comes out the way it does due to the magic of everything behaving linearly. But hopefully these two examples at least will illuminate the issue somwhat.)
