Ich versuche, eine Monte-Carlo-Simulation auf einem eindimensionalen Gitter auszuführen. Auf dem Gitter gibt es Partikel, die mit bestimmten Raten nach rechts oder links hüpfen können, die von den Populationen an der entsprechenden Gitterstelle und seinem Nachbarn abhängen. Auch links und sehr rechte Seiten gibt es ein Reservoir mit einer festen Partikelnummer. System der ODEs und erhalten Sie die erwarteten stationären Populationen an jedem Gitterstandort. Ich versuche jedoch auch, mit der dynamischen/kinetischen Monte-Carlo-Methode (Gillespie-Algorithmus) ähnliche stationäre Werte zu erzielen. Ich stecke gerade fest und verstehe nicht, wo ich mit dem Algorithmus falsch lieg. Wenn jemand eine Erfahrung hat, könnten Sie möglicherweise gemeinsame Fallstricke bestimmen? Ich teile meinen Code unten: Der erste Teil ist der ODES -Solver zum Benchmarking der Ergebnisse, und der zweite Teil ist die Monte -Carlo -Simulation, die eine falsche Lösung gibt. Schematische Darstellung des Problems, das aus Scipost Phys entnommen wurde. 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[i]) #particle jump rate from i+1 to i site
rate_R_p=GammaR*n[i]*(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[i] #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[i]+=n[i]
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()
Ich versuche, eine Monte-Carlo-Simulation auf einem eindimensionalen Gitter auszuführen. Auf dem Gitter gibt es Partikel, die mit bestimmten Raten nach rechts oder links hüpfen können, die von den Populationen an der entsprechenden Gitterstelle und seinem Nachbarn abhängen. Auch links und sehr rechte Seiten gibt es ein Reservoir mit einer festen Partikelnummer. System der ODEs und erhalten Sie die erwarteten stationären Populationen an jedem Gitterstandort. Ich versuche jedoch auch, mit der dynamischen/kinetischen Monte-Carlo-Methode (Gillespie-Algorithmus) ähnliche stationäre Werte zu erzielen. Ich stecke gerade fest und verstehe nicht, wo ich mit dem Algorithmus falsch lieg. Wenn jemand eine Erfahrung hat, könnten Sie möglicherweise gemeinsame Fallstricke bestimmen? Ich teile meinen Code unten: Der erste Teil ist der ODES -Solver zum Benchmarking der Ergebnisse, und der zweite Teil ist die Monte -Carlo -Simulation, die eine falsche Lösung gibt. Schematische Darstellung des Problems, das aus Scipost Phys entnommen wurde. 16, 029 (2024) [code]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)
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)
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
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
#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[i] #save lattice populations at each timestep return n, nInTime
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()
Ich verwende Node JS für das Backend und EJS-Vorlage für mein Frontend. Im Grunde erstelle ich ein Salon-Terminsystem, bei dem Benutzer ein Datum auswählen und auf der Grundlage des ausgewählten...