Is it equivalent to derive Gauss's law from discrete and continuous source distributions? I've seen two derivations for Gauss's law in electrostatics.
The first assumes a discrete charge distribution, the second a continuous one:


*

*Use superposition
$$\vec{E}=\sum_{i=1}^n\vec{E}_i,$$
so that
$$\oint_{\partial\Omega}\vec{E}\cdot\vec{dA}=\sum_{i=1}^n\oint_{\partial\Omega}\vec{E}_i\cdot\vec{dA}=\sum_{i=1}^n(\frac{q_i}{\epsilon_0})=\frac{Q_{tot}}{\epsilon_0}.$$
Then use Divergence theorem.

*Start with
$$\vec{E} =
\frac{1}{4\pi\epsilon_0}\int_{\mathbb{R}^3} \frac{\vec r- \vec r'}{(r-r')^3}\rho(r')dV,$$
and use the fact that
$$\nabla\cdot\frac{\vec r- \vec r'}{(r-r')^3}=4\pi\delta^3(\vec{r}-\vec{r}')$$ 
to conclude that
$$\nabla\cdot\vec{E}=\frac{\rho}{\epsilon_0}.$$
Then use Divergence theorem.
My question is whether or not the two are equivalent or if there is a difference in assuming the source as discrete or continuous.
Maybe I'm wrong and you can use (1) for continuous distributions and (2) for discrete distributions. Also, wouldn't (1) be "more" correct in that there aren't really continuous charge distributions in nature, and only approximate them as such?
 A: For what concerns most calculations, the two forms are equivalent.
Infact, I'd say that you could safely use the identification
$$ \tag{1} \int d^3 x \rho(\textbf x) \sim \sum_i q_i,$$
in all those circumstances in which it makes sense to talk about localized charges (with some exceptions though, see the last paragraph).
You can of course mathematically justify this identification using Dirac delta functions, which gives you (1) definining the charge density $\rho(\textbf x)$ as
$$ \tag{2} \rho(\textbf x) = \sum_i q_i \delta^3(\textbf x - \textbf x_i),$$
and considering that
$$ \tag{3} \int d^3x \rho(\textbf x) f(\textbf x)
= \sum_i q_i \int d^3 x \delta^3(\textbf x - \textbf x_i) f(\textbf x)
= \sum_i q_i f(\textbf x_i).$$
The use of one approximation or the other does not rely much on the nature of the charge that is being analyzed as on the experimental setup used to probe/study it (or equivalently, on which properties of the charge we are interested in).

Now for a situation in which the two approaches give very different results: consider an array of $N$ charges $q_1,...,q_N$ located at the points $\textbf x_1,...,\textbf x_N$.
The total potential energy of this system is:
$$ \tag{4} W = \frac{1}{8\pi\epsilon_0} \sum_{i \neq j} \frac{q_i q_j}{|\textbf x_i - \textbf x_j|}, \qquad i,j=1,...,N$$
where it is important to notice the $i\neq j$ in the sum, which is due to the fact that we don't want to consider the energy coming from the interaction of a point charge with itself. Considering such a thing would be very problematic: you would have a distance $|\textbf x_i - \textbf x_i|=0$ and an obviously infinite energy.
But what about the continuous version of (4)? In this case the potential energy takes the form
$$ \tag{5}W = \frac{1}{8 \pi \epsilon_0} \int \int d^3 x  d^3 x'
\frac{\rho(\textbf x) \rho(\textbf x')}{| \textbf x - \textbf x'|}.$$
But now we have two integrals, how can we implement a condition like the $i \neq j$ above? We can't, and infact this last expression gives different results from (4): it includes self-interaction terms, i.e. it includes the potential energy coming from the interaction on the charges with themselves.
To understand this, consider for example the simple situation with two charges $q_1$ and $q_2$ at points $\textbf x_1$ and $\textbf x_2$. Using (4) you get
$$ \tag{6} W = \frac{1}{8\pi \epsilon_0} \frac{q_1 q_2}{| \textbf x_1 - \textbf x_2| },$$
which is what you naively would expect.
Using instead (5) with the charge density
$$ \tag{7} \rho(\textbf x) = q_1 \delta^3(\textbf x - \textbf x_1) + q_2 \delta^3(\textbf x - \textbf x_2)$$
gives:
$$ \tag{8} W = W_{11} + W_{22} + W_{12},$$
where $W_{12}$ is the interaction term (6), while $W_{11}$ and $W_{22}$ are the self-interactions of respectively the first and second charge, which have the expressions:
$$ \tag{9} W_{11} = \frac{1}{8\pi \epsilon_0} \int d^3 x \frac{q_1^2}{|\textbf x - \textbf x_1|} \delta^3(\textbf x - \textbf x_1),$$
$$ \tag{10} W_{22} = \frac{1}{8\pi \epsilon_0} \int d^3 x \frac{q_2^2}{|\textbf x - \textbf x_2|} \delta^3(\textbf x - \textbf x_2).$$
Not only these two additional terms are not vanishing, they are infinite.
Infact, it is easily seen that even the potential energy of a single point charge is infinite when calculated through (5). So what are we to do with all these clearly wrong (are they?) infinities? Do we have to throw away the whole theory as flawed?
Obviusly not: to see that these infinities are not really a problem we just have to remember what the potential energy is for: it gives us the amount of work required/released going from one state to another.
Given that we cannot break apart a point charge (by our very definition of it), those infinite self-energies will never be exchanged with other systems. They are just (very big) constant terms added to the interaction energy, and as we know, adding a constant to the potential energy doesn't change anything in the physics.
