Python,scipy-作业-如何保存网格中粒子位置的历史记录

Question

我有以下问题:

创建一个程序，其中粒子将在以下两种情况下执行N=1000步的随机漫步:i)在1D系统中ii)在2D系统中。程序必须计算平均值(S)，其中S是粒子至少访问一次的网格位置的数量。您将运行10000次并找到10个点(每100步一个点，从0到1000)，这将是10000次运行的平均值。绘制均值(S)与时间t的关系图。

我写了下面的代码:

import scipy as sc
import matplotlib.pyplot as plt
import random
plegma=1000
grid=sc.ones(plegma)   # grid full of available positions(ones)    
for p in range(10000):
    #-------------------Initialize problem-------------------------------------------------
    his_pos=[]                  # list which holds the position of the particle in the grid
    in_pos = int(sc.random.randint(0,len(grid),1))            #initial position of particle
    means=[]                                                    #list which holds the means 
    #--------------------------------------------------------------------------------------
    for i in range(0,1000,100):
        step=2*sc.random.random_integers(0,1)-1        #the step of the particle can be -1 or 1
        # Check position for edges and fix if required
        # Move by step
        in_pos += step
        #Correct according to periodic boundaries
        in_pos = in_pos % len(grid)  
        #Keep track of random walk
        his_pos.append(in_pos)
        history=sc.array(his_pos)
        mean_his=sc.mean(history) 
        means.append(mean_his)

plt.plot(means,'bo')
plt.show()

更新 -------------------------------------

import scipy as sc
import matplotlib.pyplot as plt
import random
plegma=1000
his_pos=[] # list which holds the number of visited cells in the grid
means=[] #list which holds the means
for p in range(10000):
    #-------------------Initialize problem-------------------------------------------------
    grid=sc.ones(plegma)   # grid full of available positions(ones)      
    in_pos = int(sc.random.randint(0,len(grid),1))            #initial position of particle
    num_cells=[]       # list which holds number of visited cells during run                         
    #--------------------------------------------------------------------------------------
    for i in range(1000):
        step=2*sc.random.random_integers(0,1)-1 #the step of the particle can be -1 or 1
        # Check position for edges and fix if required
        # Move by step
        in_pos += step
        #Correct according to periodic boundaries
        in_pos = in_pos % len(grid)  
        grid[in_pos]=0  # mark particle position on grid as "visited"
        if (i+1) % 100 == 0:
            number=1000-sc.sum(grid)  # count the number of "visited positions" in grid
            num_cells.append(number)  # append it to num_cells
      his_pos.append(num_cells) # append num_cells to his_pos
      history=sc.array(his_pos)

mean_his=history.mean(1)
means.append(mean_his)

UPDATE 2 -----------------------------…

    if (i+1) % 10 == 0:
            number=1000-sc.sum(grid)  # count the number of "visited positions" in grid
            num_cells.append(number)  # append it to num_cells
    his_pos.append(num_cells) # append num_cells to his_pos
    history=sc.array(his_pos)
mean_his=history.mean(0)
plt.plot(mean_his,'bo')
plt.show()

谢谢!

Answer 1

1)是的，10000步是指所需的精度-您必须在10000次随机行走中获得t = 100, 200 ... 1000时刻访问单元的平均数量。为了获得数据，您必须为您所做的每次随机漫步累积访问的单元格数量(10000)。要做到这一点，您必须在p循环的之外初始化您的问题(即初始化his_pos和means) 。在p循环中，你应该初始化你的随机游走——得到你的网格，初始位置和你要写结果的列表。因此，问题初始化看起来像

import scipy as sc
import matplotlib.pyplot as plt
plegma=1000
his_pos=[]  # list which holds the position of the particle in the grid
means=[]    #list which holds the means 
for p in range(10000):
    #-------Initialize problem--------
    in_pos = int(sc.random.randint(0,len(grid),1)) #initial position of particle
    grid=sc.ones(plegma)   # grid full of available positions(ones)    
    his_pos.append([])

init之后，我们需要执行随机漫步——i循环。目前你只做了10步随机漫步(len(range(0,1000,100)) == 10)，但漫步应该包含1000步。因此，i循环应该是

    for i in range(1000):

在行走过程中，你必须以某种方式标记访问过的单元格。最简单的方法是将grid[in_pos]更改为0。然后，每100步需要计算访问的单元格的数量。实现方法如下

        if i % 100 == 0:
            # count the number of 0s in grid and append it to his_pos[p]

最后，在你的1000次随机游走结束时，你的his_pos将是列表的(10000*10)列表，从中应该采取按列的方法。

更新:

为了在整个运行过程中不丢失信息，我们应该将存储p -th运行结果的列表添加到包含所有结果的列表中。后者是his_pos。为了实现这一点，我们可以将空列表附加到his_pos，并在p运行期间填充结果，或者在 p运行之前初始化空列表，并在 p运行之后将其附加到his_pos 。该算法将如下所示:

#-------Initialize problem--------
plegma=1000
his_pos=[]  # list which holds the number of visited cells in the grid
means=[]    #list which holds the means 
for p in range(10000):
    #-------Initialize run--------
    in_pos = int(sc.random.randint(0,len(grid),1)) #initial position of particle
    grid=sc.ones(plegma)   # grid full of available positions(ones)    
    num_cells = []  # list which holds number of visited cells during run  
    for i in range(1000):
        # make i-th step, get particle position
        # mark particle position on grid as "visited"
        if (i+1) % 100 == 0: # on each 100th step (steps count from 0, thus i+1)
            # count the number of "visited positions" in grid
            # append it to num_cells
    # append num_cells to his_pos
# find column-wise means for his_pos

Answer 2

我不确定理解问题是什么(在你的方程中考虑t的时间?， point指的是什么?)，但相对于您试图执行的操作，我理解您必须在最终的10000x10 means数组中包含每次迭代产生的每个或10-means mean_his数组。

每个mean_his数组都是通过100个步骤从一个1D数组中制备的。我将用一个包含十个步骤的数组来说明这一点，这些步骤必须每两步平均一次(而不是每100步平均1000步):

>>> his_pos = [1,2,3,4,5,6,7,8,9,10]    #the ten positions
>>> history = np.array(his_pos)
>>> history
array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10])
>>> ar2 = np.reshape(history, (-1,2))   # group in two in order to calculate mean
>>> ar2
array([[ 1,  2],
       [ 3,  4],
       [ 5,  6],
       [ 7,  8],
       [ 9, 10]])
>>> mean_his = ar2.mean(1)
>>> mean_his
array([ 1.5,  3.5,  5.5,  7.5,  9.5])
>>>  

然后你将mean_his附加到means 10000次，并计算类似的平均值(注意，意味着必须在外部循环之外初始化，以避免每次重复初始化)。