At the end of this post, I have listed the code I wrote. But I'm not going to go through it, for two reasons. First, I know you would be bored (me too). Second, it isn't very polished, and has a few issues. But the results give a picture of the distribution that I think must be correct. Maybe one of you can tell me if that's true or not.
We can run the simulation with the option to report lots of information as it runs. I always use the same names for common variables. I use v for "verbose", L for "list", D for "dictionary", and so on.
Here is some output from the code listing below with v=True. We show the particle positions in 40 bins, updating only when positions change, and print detailed info when something more interesting happens (reversal, collision). The particles are numbered and identified with the mass type 'm' or 'M'.
How do we test the code? This is a hard problem. To make a start, we can run a simulation with X (box size) and the number of particles very small, and look carefully whenever things change. Two other tests: in a 1D box, the order of particles should never change, and also the total kinetic energy should not change with time. The latter requirement is met by this code. However, the former invariant is occasionally broken. The code reports this when it occurs. In the line that starts "96" above, we see "41235 14235", which indicates that the order switched as shown. On the following cycle it switched back, but this doesn't always happen. It gets harder when there are 3 or more particles within the same small interval. I have not solved this problem, but these "quantum"-like events are rare, and I believe they don't invalidate the statistic we are going to measure. (I can't prove that, of course, but the numbers are quite small).
I've run the simulation with more particles and a slightly larger space to make the velocities change more frequently. For these tests, I assigned the same velocity to every particle. I want to show you statistics for these runs, with various mass ratios as mentioned above. Two things are of interest: the standard deviation of the distribution of velocities, and the histograms. Here are the results (in R) for the stdevs:
I will have to look up the values that statistical mechanics theory predicts. It seems like a good guess that the sd of the small particles goes like the sd of the large particles * √M/n. But, I would have to do some other tests to be sure. For example, we have a number of m and M particules that is random, but we did random.seed(157), so it is actually set to one specific value. Finally, notice that the distributions are all normal. I hope this does not come as a surprise.
And here are the plots.
M/m = 2
M/m = 3
M/m = 4