What is the evolution equations of spins with an Ising model hamiltonian? Given an Hamiltonian $H$ defined with respect to coordinates $q_i$ and momenta $p_i$ we can write the Hamilton equations:
$ \begin{cases}
   \dot{q}_i = -\frac{\partial H}{\partial p_i}\\
   \dot{p}_i = \frac{\partial H}{\partial q_i}
\end{cases} $
Can we construct similar equations for the spins in an Ising model? To be more precise, the hamiltonian is the following:
$ H(\{S_i\}) = -J\sum_{<ij>} S_iS_j \qquad S_i = \pm 1 $
This should be something of the kind:
$ S_i(t) = f(H(\{S_i(t)\}),\{S_i(t)\}) $
I would also insert the derivative but it will always be zero or undefined where $S_i(t)$ changes sign.
Here is how I tried to find the equation, or rather, managed not to find it. Let's make things easier by considering the 1D Ising model on a chain, so that we can write explicitly:
$ H = -J\sum_{i=1}^N S_iS_{i+1} $
If we change a spin then the energy will necessarily change as well. If the system is isolated then we will have to change at least two spins. This is already problematic for the ground states, where they are all spin up or spin down. If one spin changes, in order to have the same energy, all other spins will have to change as well. This means that the time evolution has to be such that there is no evolution for the ground state or the evolution is periodic between the two. However, this motion would have to be periodic even without a given time-scale, there is in fact no $\hbar$ to get $E/\hbar$ as a time.
The first excited state (suppose we have free boundary conditions) will be given by one broken link, and we can place this anywhere along the chain. If the time evolution is local, meaning that a spin can change only its neighbours, then there is no preference between its left and right neighbours. This clearly means that the evolution can't be deterministic and we can't find $S_i(t)$, but it might makes sense to have an evolution equation for a probability $p(S_i,t)$. This is still problematic since $S_i$ can have only discrete values, but something akin to the Liouville equation for $\rho(q,p)$ might make sense.
Am I wrong? Is there a possible equation we can have which determines the time evolution of spins?
What happens if we consider O(N) models, where have continuous spins?
To study numerically thermodynamic properties of systems, one often does simulations where the system evolves isolated or in contact with a bath. Has this been done for Ising models (or other similar magnetic models)?
 A: Classical two-level system is problematic
The problem is not specifically related to the Ising model, but to considering a two-level system using classical Hamiltonian equations of motion. Since the states are discrete, using continuous formulation in terms of coordinates and momenta doesn't work. There are however various situations where the model would make perfect sense dynamically, but this requires some adding additional information:

*

*From the quantum point of view two-level systems do not represent anything particularly difficult. Associating Ising spins with Pauli matrices makes writing the dynamic equations straightforward.

*Alternative approach is considering spins as rotating arrows, which interact only via their z-components.

Dynamics is not important
The point of the Ising model and various other models inspired by it is describing phase transitions/critical phenomena, which are thermodynamic rather than dynamic phenomena. The whole point of thermodynamic/statistical physics approach to modeling is considering thermodynamic equilibrium, while neglecting the interactions/processes responsible for establishing this equilibrium. This is why most models used in statistical physics simply do not contain these interactions: the ideal gas or a non-interacting collection of two-level systems are obvious examples.
Numerical modeling
In some numerical modeling approaches the interactions are added to Ising or similar models to permit the system to achieve thermodynamic equilibrium. Note that the interactions added in such cases by no means have to be realistic, but simply allowing the transitions between the states. A Monte Carlo approach randomly flipping spins works quite well for the Ising model.
Note that such artificial interactions are not reserved solely to numerical approaches. E.g., the concept of a black body which absorbs all the radiation incident on it is a trick to be able to obtain the thermal radiation spectrum (although one can obtain it completely neglecting the interactions, by considering free photon gas in thermodynamic equilibrium) - see discussion in this thread.
