Microscopic and macroscopic description of spin waves Hamiltonian
Consider the one-dimensional Heisenberg ferromagnet specified by the Hamiltonian
$$H = -\frac{|J|}{2}\sum_{i,\delta} \mathbf{S}_i\cdot \mathbf{S}_{i+\delta}.$$
Here $i$ labels the spin sites and $\delta$ connects a site to its nearest neighbour. $\mathbf{S}_i$ is the spin operator living at site $i$.
Microscopic approach
By utilizing the Holstein-primakoff transformation and performing a Fourier transformation, we find that the eigenstates are magnons. We usually say that when magnons are present in the system the spins are precessing around some fixed axis (let us choose the z-axis) giving rise to a spin wave.
Numerical approach
Now, say we wish to study the spin wave numerically. I guess we can do this by solving the (undamped) LLG equation
$$\frac{d\mathbf{S}_i}{dt} = \mathbf{S}_i \times \mathbf{\mathcal{H}},$$
where the effective field is defined as $\mathbf{\mathcal{H}} = -\delta H/\delta \mathbf{S}_i$. I would then expect to get an effective field that points in the $z-$direction such that the spins $\mathbf{S}_i$ are precessing around the $z-$axis. However, when I compute the effective field I obtain
\begin{equation}
\begin{split}
&\mathcal{H}_i^x = |J|(S_{i-1}^x + S_{i+1}^x),\\
&\mathcal{H}_i^y = |J|(S_{i-1}^y + S_{i+1}^y),\\
&\mathcal{H}_i^z = |J|(S_{i-1}^z + S_{i+1}^z).
\end{split}
\end{equation}
Clearly, in an excited state, the effective field has components in all three cartesian directions not just the $z$-direction. In addition the effective field seems to have some time dependence because the components of $\mathbf{S}_i$ are time dependent.
The question
So to sum up: when solving the equation of motion numerically it seems that the spins are precessing around some time-dependent effective field. However, the usual microscopic picture is that the spins are precessing around the $z$-axis. How can I reconcile these two pictures?
Attempt 1
I can get the effective field to point in the $z$-direction if I choose that
\begin{equation}
\begin{split}
&S_{i-1}^x + S_{i+1}^x = 0,\\
&S_{i-1}^y + S_{i+1}^y = 0.
\end{split}
\end{equation}
But then it seems that neighbouring spins are always out of phase by $\pi/2$ radians, which I do not think is in accordance with the microscopic description.
Attempt 2
One could argue that the $x$ and $y$ components of $\mathbf{S}_i$ are very small compared to the $z$-component. In that case the effective field will approximately point in the $z$-direction. Is this a correct way of looking at the situation?
 A: Firstly, because the hamiltonian is invariant to rotations (in this isotropic case), there is no predefined $z$-axis. However, the LLG equations exactly conserve the total magnetization
$$
\mathbf{S} \equiv \sum_{n=1}^N \mathbf{S}_n
$$
so in any finite-sized system there is de facto a unique axis defined by $\mathbf{S}$.
The conservation law can be easily established from
$$
\dot{\mathbf{S}}=\sum_{n=1}^N \mathbf{S}_n \times 
( \mathbf{S}_{n-1}+\mathbf{S}_{n+1} ) = \mathbf{0}
$$
noting that every cross product turns up twice in the sum,
with opposite signs,
assuming periodic (wraparound) boundary conditions.
This leads us to a couple of easy ways to break the symmetry
(and the conservation law).
Firstly, if we apply a small external field in the $z$-direction,
it is easy to see that a uniform precessional motion about $z$ is added 
to the above equations of motion.
Alternatively,
if we abandon periodic boundaries and have two fixed spins 
$\mathbf{S}_0$ and $\mathbf{S}_{N+1}$ at the ends of our chain
and (for the sake of argument)
have them both pointing in the $z$ direction,
then the total magnetization $\mathbf{S}$ will no longer be conserved
and will rotate about the $z$ axis:
$$
\dot{\mathbf{S}}= \mathbf{S}_1 \times \mathbf{z} + \mathbf{S}_{N}\times \mathbf{z}
$$
The effect will be small, because it comes from the boundaries
(this is a general consequence of the conservation law)
and it depends on the fluctuating $x$ and $y$ components of spins $1$ and $N$.
However, you were specifically interested in the spin waves.
Spin-wave-like solutions of the 
classical Heisenberg spin-chain equations have been
discussed in the literature.
For example in JAG Roberts and CJ Tompson J Phys A, 21, 1769 (1988) the following exact solutions are presented
$$
\mathbf{S}_n(t) = [\mathbf{a}\cos\phi_n(t)+\mathbf{b}\sin\phi_n(t)]\sin\theta
+ \mathbf{c}\cos\theta
$$
where 
$\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$ 
are a set of basis vectors for the spins,
$\phi_n(t)= np-\omega t$ and $\omega=2(1-\cos p)\cos\theta$,
and I have taken the liberty of redefining their angles
so that they are the usual polar ones.
These go over to the continuum spin-wave solutions by replacing
$p\rightarrow pa$ and $\omega\rightarrow\omega/a^2$ where $a$ is the 
lattice spacing, and allowing $a\rightarrow 0$.
Importantly, $\theta$ is a constant and they point out

These solutions describe the precessional motion of each spin about
  the $\mathbf{c}$ axis.

It's clear that $\mathbf{c}$ is the direction of the total (conserved)
magnetization. 
If we do the sum $\sum_n\mathbf{S}_n(t)$,
and impose the restriction that $p$ is an integer multiple of $2\pi/N$
for consistency with the periodic boundaries,
the time dependent terms in the $\mathbf{a}$ and $\mathbf{b}$
directions sum up to zero.
Hence
$$
\sum_{n=1}^N\mathbf{S}_n(t) \equiv \mathbf{S} = N\cos\theta \, \mathbf{c}
$$
So there is a precession of each of the spins around a fixed axis,
although the total magnetization vector is conserved (and lies along that axis).
CAVEAT: in that paper they also state 

Unlike the continuum isotropic case, however, it is almost certain
  that the above spin-wave solutions of the discrete equations  are not
  the most general solutions of the form 
  $\mathbf{S}_n=\mathbf{s}(pn-\omega t)$.

But this example does illustrate how precessing spin-waves can arise
in the numerical solutions,
without imposing restrictions on the values of individual spin components. There may be other kinds of solutions too, of course.
Incidentally, if you are interested in the numerical solution of this
kind of model on a lattice,
one kind of algorithm is described by 
J Frank, W Huang and B Leimkuhler, J Comp Phys, 133, 160 (1997)
and by
M Krech, A Bunker and DP Landau, Comp Phys Commun, 111, 1 (1998)
which is also available here. 
Basically,
the spins are divided into even and odd subsets,
and each subset is advanced in turn 
(in the fixed field of the neighbouring spins).
The method generalizes straightforwardly to 2D and 3D lattices.
It allows quite long time steps,
and exactly conserves spin norm and total energy, 
but not total magnetization, which fluctates.
If magnetization conservation is important, a more traditional predictor-corrector method may be preferable.
