How to animate a randomly updated plot

Question

I have created a simple animation to explain to students how temperature affects the occupation of energy levels by using Animate to try to continually apply a Metropolis update to a system of particles in a three-level energy system. The code looks like this:

L = 5;
s = Table[{i, 0}, {i, 1, L}];
e = {0, 1, 4}; (*allowed energy levels*)
excite[n_, T_] := (
   x = RandomInteger[{1, 3}]; 
   de = e[[x]] - s[[n, 2]];
   If[de < 0, 
      s[[n, 2]] = e[[x]], 
      p = RandomReal[];  
      If[p < Exp[-de/T],  
         s[[n, 2]] = e[[x]]
        ]
     ]
  )

plt = Plot[e, {i, 0, L + 1}, AspectRatio -> 1.3, ImageSize -> Small];

Animate[
 Manipulate[
  excite[t, T];
  Show[plt, ListPlot[s, PlotStyle -> PointSize[0.05], Filling -> Axis]]
  , {T, 0.1, 11}]
 , {t, 1, L, 1}]

where the temperature can be controlled with Manipulate. However, this does not act the way I would expect whatsoever. Sometimes it works rather nicely, but it sometimes will just freeze and stop updating until I try to change the temperature, and often continues to update even after I have paused the animation. Also, the rate at which it updates does not seem to be influenced at all by the step size or refresh rate of the animation. I have also tried putting the Manipulate outside the Animate and it has similar behavior. Is there a better way to do this?

To be clear, I want the animation to run the excite routine for each particle at a reasonable rate. Even as I have it now, at each time step of the animation one particle should be excite'd once per update, but when I run it sometimes a single particle will jump up and down multiple times at a single time step, I have no idea what Mathematica is doing (also I would prefer if it applied excite to each one at once before updating the plot). When T is low, all the particles should stay at the bottom, but when T is high (above the maximum e), all the levels should on average have the same number of particles, but even when I crank T up they still seem to stay in the bottom level.

Here's an example of what it looks like:

Answer 1

You may use Manipulate and Scan without Animate.

Manipulate[
 Scan[excite[#, T] &, Range@Length@s];
 Show[plt, ListPlot[s, PlotStyle -> PointSize[0.05], Filling -> Axis]],
 {T, 0.1, 11, Appearance -> "Open"},
 TrackedSymbols :> {T}
 ]

Scan is used to recalculate the states of s at each T since the result of the calculation does not need to be returned as it is stored in s. Map could be used instead (e.g. excite[#, T]& /@ Range@Length@s;)

TrackedSymbols is used to explicitly identify which symbols should trigger an update. It is needed since T is in Scan which makes Mma uncertain on when it should update.

Hope this helps.

Answer 2

Here I added the variable maxEns this controls the number of excitations into the system, and thus the length of time that the simulation will run.

ClearAll["Global`*"]
L = 5;
s = Table[{i, 0}, {i, 1, L}];
e = {0, 1, 4}; (*allowed energy levels*)
maxEns = 150; (*This is the max number of excitations*)

Next I defined a function, which takes Temperature as its input, which will apply a random excitation to each of the elements in the list s. (Here I set each of the particles to recieve a different randomly chosen excitation energy from e)

excite[T_] := (
  randEnVals = Table[RandomChoice[e], Length[s]];
  For[n = 1, n <= Length[s], n++,
   de = randEnVals[[n]] - s[[n, 2]];
   If[de < 0,
    s[[n, 2]] = randEnVals[[n]],
    If[RandomReal[] < Exp[-de/T],
     s[[n, 2]] = randEnVals[[n]]]]];
  s)

The next function creates a list of all of the excitations as each additional excitation is applied

excitations[T_, numbEns_] := excitations[T, numbEns] = (
   s = Table[{i, 0}, {i, 1, L}];
   templist = {s};
   For[iEns = 1, iEns < numbEns, iEns++,
    AppendTo[templist, excite[T]]];
   templist)

We now define the plots

plt = Plot[e, {i, 0, L + 1},
   AspectRatio -> 1.3,
   ImageSize -> Small];

plt2[eVals_] := ListPlot[eVals,
   PlotStyle -> PointSize[0.05],
   Filling -> Axis];

Finally using manipulate we can plot the excitations over time (in units of number of excitations).

Manipulate[
 Show[{plt, Evaluate[plt2[excitations[T, maxEns][[t]]]]}], {T, 0.1, 
  11}, {t, 1, maxEns, 1}]

The desired result of having the t value run continuously as the T value is manually adjusted can be accomplished by clicking on the + sign next to the t slider. Here the speed can also be adjusted.