如何在Python的TASEP模拟中统计被占据格点的数量？

Hey there! I see you're running a TASEP simulation and thought you had to count occupied lattice sites from the generated plot—luckily, there's a much simpler and more accurate way using your existing simulation data. Let me walk you through how to do it.

Your state array already stores exactly whether each lattice site is occupied (1) or empty (0). Since the actual lattice sites are in state[1] to state[L] (note: state[0] is just an entry buffer that's never occupied), we can calculate the number of occupied sites with a single line of code.

Quick: Get the Occupied Count at Any Step #

To get the number of occupied sites at any point in the simulation, just sum the relevant part of the state array:

occupied_count = sum(state[1:L+1])

Adding This to Your Simulation #

Here are two common ways to integrate this into your existing code:

Option 1: Print the Final Occupied Count #

If you just want to know how many sites are occupied at the end of the simulation, add this line right after the main loop finishes:

final_occupied = sum(state[1:L+1])
print(f"Final number of occupied lattice sites: {final_occupied}")

Option 2: Track Occupancy Over Time (for Visualization) #

If you want to see how the number of occupied sites changes as the simulation runs, you can collect this data at the same intervals you generate your configuration plots. Here's the modified code:

import random, pylab, math
random.seed()
L=100 # Number of lattice sites
alpha=.2 # Rate of entry
beta=.4 # Rate of exit
Ntime=200000 # Simulation steps
state=[0 for k in range(L+1)]
# Initialize a list to store occupancy data over time
occupancy_history = []

for iter in range(Ntime):
    k=random.randint(0,L)
    if k==0:
        if random.random() < alpha:
            state[1]=1
    elif k==L:
        if random.random() < beta:
            state[L]=0
    elif state[k]==1 and state[k+1]==0:
        state[k+1]=1
        state[k]=0
    
    if iter%2000 == 0:
        # Calculate current occupancy and store it
        current_occupancy = sum(state[1:L+1])
        occupancy_history.append((iter, current_occupancy))
        
        # Your existing configuration plot code
        yaxis=[]
        for i in range(L):
            if state[i+1]==1:
                yaxis.append(i)
        xaxis=[iter for k in range(len(yaxis))]
        pylab.plot(xaxis,yaxis,'r.')

# Print the final occupied count
final_occupied = sum(state[1:L+1])
print(f"Final number of occupied lattice sites: {final_occupied}")

# Optional: Create a separate plot for occupancy over time
pylab.figure()
steps, counts = zip(*occupancy_history)
pylab.plot(steps, counts, 'b-')
pylab.xlabel('Number of steps')
pylab.ylabel('Number of occupied sites')
pylab.title('TASEP Occupancy Over Time')

pylab.show()

Why This Works #

The sum() function adds up all the 1s in state[1:L+1]—each 1 corresponds to an occupied site, so the total sum is exactly the count you need. This method is far more efficient and precise than manually counting from a plot!

内容的提问来源于stack exchange，提问作者Mr. Moe