First some thoughts and not an answer. It is helpful to re-think the problem in terms of the gaps between the stones. You have $n$ slots with $k$ stones. So in the case with boundaries you have $(k+1)$ stretches of "in between space", each of which we will denote $a_i = x_i - x_{i-1} - 1$. The sum of all these $a_i$ add up to $n-k$. The dynamical system in terms of $a_i$ can be written as follows: let $a_0 = a_1$ and $a_{k+2} = a_{k+1}$ (translating the Dirichlet conditions of your original system to a Neumann boundary condition), and demand
$ a_i(t+1) - a_i(t) = \lfloor \frac{a_{i+1}(t) - a_i(t)}{2}\rfloor - \lfloor \frac{a_i(t) - a_{i-1}(t)}{2}\rfloor $
which is in fact (more-or-less) a discrete heat equation for $a:\mathbb{Z}\to\mathbb{Z}$. (The discrete Laplacian may also be written as $\lfloor (a_{i+1} + a_{i-1})/2\rfloor - a_i$ which differs from the RHS about by at most 1. The version above has instead the nice property that it preserves the $\ell_1$ norm of $a$, that is, the sum of all $a_i$ is preserved under this evolution.)
So we can apply similar intuition behind the heat equation: namely that one can expect convergence of the solution to a harmonic function. (The discretized setting however suggests that the limiting solution will not necessarily have $a_i(t+1) = a_i(t)$, but small fluctuations like $|a_i(t+1) - a_i(t)| \leq 1$ should be expected.)
Now, the harmonic functions on an intervale are the constant functions and the linear functions. The latter is ruled-out by the Neumann boundary condition. However, because of your discretization and the preferential placement of the stones to the left spot when two are equal, the modified Laplacian given above admits a harmonic function as described by mjqxxxx in his comment: since differences of +1 between consequentive $a_i$s are treated as if they don't exist. So you can have solutions which locally looks like it is equidistributed (left side and right side differs by a small amount) but globally isn't.
The periodic model, however, doesn't necessarily help. Consider the case with 4 stones and 12 slots. Take the initial distribution (O for stone, u for empty slot)
OOuuuuOOuuuu
this becomes
uuuOOuuuuOOu
and the pattern oscillates, never settling down. This example can be generalized to large $n$ also. This oscillation represents a failure of the intuition for passing from the continuous heat equation to the discrete case. When the spatial wavelength of the oscillation is small enough (on the order of the lattice size), artefacts from the discretisation starts to interfere with the intuition. We will see below that this artefact can only manifest for when there are even number of stones in the periodic case, can cannot manifest in the interval case.
Now more concrete answers. Using the property of the heat equation, one can be inspired to check for dissipation laws (in other words, apply the maximum principle for parabolic differential equations). And indeed, by the definition of the system, in either the interval case or the periodic case, you can easily check that
- The value $a_{min} = \min a_i$ is non-decreasing, and $a_{max} = \max a_i$ is non-increasing.
- Suppose $a_i$ attains the minimum at time $t$. Then $a_i(t+1) = a_i(t)$ if and only if $a_{i-1}(t) = a_i(t)$ and either $a_{i+1}(t) = a_i(t)$ or $a_{i+1}(t) = a_i(t) + 1$. Similarly for the maximum.
Reasoning directly from point 2 you obtain that, in the interval case (your original question), any distribution will converge toward a distribution that looks like the following: there exists $k_m$ and $k_M$ such that for all $i \leq k_m$, you have $a_i = a_{min}$. For all $i \geq k_M$, you have $a_i = a_{max}$. In between $k_m$ and $k_M$, you have $a_{i+1} - a_i = 1$. The values $k_m$ and $k_M$ can be $0$ or $k+1$, reflecting a limiting distribution that is a constant.
Parity
Let us re-write the discrete heat equation. A bit of computation shows that we can write it in the following way:
$ a_i(t+1) = \lfloor \frac{a_{i+1}(t) + a_{i-1}(t) + \left( a_{i-1}(t) - a_i(t) \mod 2\right)}{2}\rfloor$
which says that the value of $a_i$ at $t+1$ is the average of its neighbors at $t$, rounded up if $a_{i-1}-a_i$ is odd, and rounded down if it is even. So in the periodic case with $k$ even, there are even number of $a_i$s, which can be then partitioned into two (effectively non-interacting) subsets (the evens and the odds) that results in the oscillating behaviour noted above. In the case of the interval by the Neumann boundary condition, and in the periodic case with the odd number of stones, there cannot by such partition of the system into two independent subsystems.
(to be cont.)