Many Particle System - Setting up the Initial State

Simulation and Modeling (CSCI 3010U)

In-class Exercise

Introduction

Consider the code below. It places particles on a regular grid.

import numpy as np
from numpy import random
from scipy.integrate import ode
from matplotlib import pyplot as plt
from matplotlib import animation

def rect_pos(nx, ny, Lx, Ly):
    '''
    Create 4 x n_particles state vector
    '''
    n_particles = nx * ny
    state = np.zeros([4, n_particles])
    dx = Lx / float(nx)
    dy = Ly / float(ny)
    for ix in range(nx):
        for iy in range(ny):
            i = ix + iy * nx
            state[0,i] = dx * (ix + 0.5) 
            state[1,i] = dy * (iy + 0.5) 
    return state

nx, ny, Lx, Ly = 15, 20, 10, 10

state = rect_pos(nx, ny, Lx, Ly)

plt.figure(1)
plt.axes(xlim=(-1,Lx+1), ylim=(-1,Ly+1))
plt.title('Initial positions')
plt.xlabel('x')
plt.ylabel('y')
plt.plot(state[0,:], state[1,:], '.')
plt.show()

Task 1

Your goal is to complete a function that places particles at random locations between $ [ 0, L_x ] \times [ 0, L_y ]$. Please ensure that no two particles are within radius $r$ of each other.

def rand_pos(nx, ny, Lx, Ly, r):
    # TO DO
    n_particles = nx * ny
    state = np.zeros([4, n_particles])
    return state

nx, ny, Lx, Ly, r = 15, 20, 10, 10, 2

state = rand_pos(nx, ny, Lx, Ly, r)

plt.figure(1)
plt.axes(xlim=(-1,Lx+1), ylim=(-1,Ly+1))
plt.title('Initial positions')
plt.xlabel('x')
plt.ylabel('y')
plt.plot(state[0,:], state[1,:], '.')
plt.show()

Task 2

Can you comment on sample efficiency of your approach, i.e., what fraction of samples are wasted? Can you plot sample efficiency as a function of radius $r$?