The stability analysis is made considering two close trajectories, if the difference between the state is guaranteed to reduce after one iteration, it means that the method is stable, the numerical errors introduced in one iteration (rounding errros) will not accumulate.
The numerical solutions satisfy the equation
$$ u_i^{k+1} = r\,(u_{i+1}^k + u_{i+1}^k) + (1 - 2r)u_{i}^k$$
Now consider the evolution of the same initial state with a perturbation $\Delta u_i^{k}$, at iteration $k$ in the update equation and you get
$$ u_i^{k+1} + \Delta u_i^{k+1} = r(u_{i+1}^k + \Delta u_{i+1}^k + u_{i-1}^k + \Delta u_{i-1}^k) + (1 - 2r)(u_i^k + \Delta u_i^k) $$
After one iteration the difference between the perturbed trajectory and the original trajectory satisfies
$$ \Delta u_i^{k+1} \le r(|\Delta u_{i+1}^k| + |\Delta u_{i-1}^k|) + |1 - 2r|(|\Delta u_i^k|) $$
Let $\Delta u^k = \max_{i} |\Delta u_i^k|$, then we can say
$$ \begin{eqnarray} |\Delta u^{k+1}| &\le& |r|(|\Delta u^k| + |\Delta u^k|) + |1 - 2r||\Delta u^k| \\ &\le& |2r| |\Delta u^k| + |1 - 2r||\Delta u^k|\end{eqnarray} $$
If $0 \le r \le 1/2$ we have
$$ |\Delta u^{k+1}| \le 2r |\Delta u^k| + (1 - 2r)|\Delta u^k| = |\Delta u^k| $$
And that means that the error will reduce each iteration.
When you compute using finite precision, at each iteration a new perturbation is added to $\Delta u_i^{k+1}$ and that is guaranteed to be small, if the method is stable it will keep small, if the method is not stable i.e. $|\Delta u^{k+1}|$ is not guaranteed to be smaller than $|\Delta u^k|$, then the rounding errors may be amplified in the subsequent iterations, being impossible to get an accurate solution.