Правильно выполнить кинетическую Монте -Карло на 1D решетке ⇐ Python

[Anonymous]

Я пытаюсь выполнить симуляцию Монте-Карло на одномерной решетке. На решетке существуют частицы, которые могут прыгать вправо или слева с определенными скоростями, которые зависят от популяций на соответствующем участке решетки и его соседа. Также с самой левой и очень правой стороны есть резервуар с фиксированным номером частиц. Система ODES и получить ожидаемую стационарную популяцию на каждом участке решетки. Тем не менее, я также пытаюсь достичь аналогичных стационарных значений с помощью динамического/кинетического метода Монте-Карло (алгоритм Gillespie). В настоящее время я застрял и не понимаю, где я ошибаюсь с алгоритмом. Если у них есть опыт, не могли бы вы определить общие ловушки? Я делюсь своим кодом ниже: первая часть - решатель ODES для сравнения результатов, а вторая часть - моделирование Монте -Карло, которое дает неверное решение. < /P>
Я также прикрепляю Схематический сюжет проблемы, взятой из Scipost Phys. 16, 029 (2024).import numpy as np
import matplotlib
from matplotlib import pyplot as plt
from matplotlib.pyplot import cm
from scipy.integrate import solve_ivp
matplotlib.rcParams['figure.figsize'] = (10, 7)

###--------DEFINING PARAMETERS-----------###########
lattice_length=10
cmap = cm.terrain(np.linspace(0, 0.8, lattice_length))

lattice=np.linspace(1, lattice_length, lattice_length)
initial_condition=np.zeros(lattice_length)

GammaL=1 #hopping rate to the left
GammaR=0 #hopping rate to the right
kappaL=kappaR=1 #hopping rates between the very left and right lattice sites and their corresponding left and right reservoirs

nRight=1 #right reservoir population
nLeft=0 #left reservoir population

######-----------SOLVING RATE EQUATIONS WITH ODE---------############
time_ODE=np.linspace(0,100,10000)

def sol_fun_eq():

def dndt(t,V):

dndt=np.zeros(lattice_length)

J_01=kappaL*V[0]*(1+nLeft)-kappaL*nLeft*(1+V[0])
J_12=GammaL*V[1]*(1+V[0])-GammaR*V[0]*(1+V[1])

J_Last=kappaR*nRight*(1+V[-1])-kappaR*V[-1]*(1+nRight)
J_LastButOne=GammaL*V[-1]*(1+V[-2])-GammaR*V[-2]*(1+V[-1])

dndt[0]=J_12-J_01
dndt[-1]=J_Last-J_LastButOne

for p in range(1,lattice_length-1):
J_in = GammaL*V[p+1]*(1+V[p])-GammaR*V[p]*(1+V[p+1])
J_out = GammaL*V[p]*(1+V[p-1])-GammaR*V[p-1]*(1+V[p])
dndt[p]= J_in - J_out
return [dndt]

sol = solve_ivp(dndt, [time_ODE[0], time_ODE[-1]], initial_condition, method='LSODA', t_eval=time_ODE)
return sol

solution_eq=sol_fun_eq()

steady_state=np.zeros(lattice_length)
for p in range(lattice_length):
steady_state[p]=solution_eq.y[p][-1]

# Plot the results in time
plt.figure(figsize=(12, 8))
for p in range(lattice_length):
plt.plot(solution_eq.t, solution_eq.y[p], label=f'$n_{{{p}}}$', color=cmap[p])
plt.title("Time traces")
plt.xlabel("Time")
plt.ylabel("Populations")
plt.show()

# Plot the steady state
plt.figure()
plt.plot(lattice, steady_state)
plt.scatter(lattice, steady_state)
plt.xlabel("Lattice sites")
plt.ylabel("Populations")
plt.title("Steady state populations")
plt.show()

###--------------------MONTE CARLO SIMULATION----------##########
timeSize=500
time_MC=np.linspace(0, timeSize, timeSize)
SampleRange=2000 #number of sampling

def MC_sim(dummy):

n=np.zeros(lattice_length)
nInTime=np.zeros((timeSize, lattice_length))

for tIdx in range(timeSize):

rates=[]
event_list=[]
#calculate boundary rates on the very left and right sites
rate_01=kappaL*n[0]*(1+nLeft) #particle jump rate from site #1 to left reservoir
rate_10=kappaL*nLeft*(1+n[0]) #particle jump rate from left reservoir to site #1

rate_LastR=kappaR*nRight*(1+n[-1]) #particle jump rate from right reservoir to last lattice site
rate_RLast=kappaR*n[-1]*(1+nRight) #particle jump rate from last lattice site to right reservoir

rates.extend([rate_01, rate_10, rate_LastR, rate_RLast])
event_list.extend(["rate_01", "rate_10", "rate_LastR", "rate_RLast"])

for i in range(lattice_length-1):

rate_L_p=GammaL*n[i+1]*(1+n) #particle jump rate from i+1 to i site
rate_R_p=GammaR*n*(1+n[i+1]) #particle jump rate from i to i+1 site

rates.extend([rate_L_p, rate_R_p])
event_list.extend([f"L{i}", f"R{i}"])

sumRates = np.sum(rates)
probabilities=rates/sumRates
events=np.random.choice(len(rates), p=probabilities)

#update the lattice site
if events==0:
n[0]-=1
elif events==1:
n[0]+=1
elif events==2:
n[-1]+=1
elif events==3:
n[-1]-=1
elif events>3:
event_idx=event_list[events]
transition=event_idx[0]
lattice_idx=int(event_idx[1])
if transition=="L":
n[lattice_idx]+=1
n[lattice_idx+1]-=1
elif transition=="R":
n[lattice_idx]-=1
n[lattice_idx+1]+=1
else:
raise Exception("Coding error")
for i in range(lattice_length):
nInTime[tIdx, i]=n #save lattice populations at each timestep
return n, nInTime

nSummed=np.zeros(lattice_length)
nInTimeSummed=np.zeros((timeSize, lattice_length))

for s in range(SampleRange):

n, nInTime = MC_sim(0)

for i in range(lattice_length):

nSummed+=n
nInTimeSummed[:,i]+=nInTime[:,i]

n_avg=nSummed/SampleRange
nInTimeAvg=nInTimeSummed/SampleRange

plt.figure(figsize=(12,7))
for i in range(lattice_length):
plt.plot(time_MC, nInTimeAvg[:, i], label=f"# {i}")
plt.legend()
plt.title("Averaged time traces")
plt.xlabel("Time (arbitrary)")
plt.ylabel("Populations")
plt.show()

plt.figure()
plt.plot(lattice, n_avg)
plt.xlabel("Lattice sites")
plt.title("Averaged steady state populations")
plt.ylabel("Populations")
plt.show()


Подробнее здесь: https://stackoverflow.com/questions/794 ... 1d-lattice