How To Create 'billiard Ball' Reflection Boundary Condition In Python?

According to Erwin Schrodinger (in What is Life?), diffusion can be explained entirely by the random motion of particles. I want to test this myself by create a program the creates

Solution 1:

You don't need to calculate the reflection angle, just decompose the problem in two: one for x and one for y. In both cases, you need the particle to "go back" when it exceeds the boundary.

I've done this time ago for an excercise studying particle density in fluids. The easiest thing is to consider a (0, 1) boundary in both directions. The following code should do it (tip: a proper use of abs will create the equivalent of a reflection):

x0 = [.1, .9]
delta = [-0.2, 0.3]
x1 = [(1-abs(abs(xi + di)-1)) for xi, di in zip(x0, delta)]
print(x1)
# 0.1, 0.8
#or using numpy:
x1 = 1-np.abs(np.abs(np.asarray(x0) + np.asarray(delta))-1)
print(x1)
>> [0.09999999999999998, 0.8]
   array([0.1, 0.8])

I'm assuming from your question that you are neglecting particle-particle collision and particle-particle "non-superposition"

Solution 2:

Here is a simple implementation. I change the movement vectors only every tenth step, that way one can visually check boundary reflections. Particles flash red when the movement vectors are updated.

The trick as described ħere is to "unroll" the bounding box. Instead we let the particles move unconstrained and then fold space into the bounding box.

import numpy as np
import pylab
from matplotlib.animation import FuncAnimation

xy = np.random.uniform(-1, 1, (2, 200))
xy[0, :160] = np.abs(xy[0, :160])
xy[0, 160:] = -np.abs(xy[0, 160:])
xy += 1

f, a = pylab.subplots()
pxy, = pylab.plot(*xy, 'o')

def init():
    a.set_xlim(0, 2)
    a.set_ylim(0, 2)
    return pxy,

def update(frame):
    global inc, xy
    if frame % 1 < 0.01:
        inc = np.random.normal(0, 0.01, xy.shape)
        pxy.set_markerfacecolor('red')
    elif frame % 1 < 0.11:
        pxy.set_markerfacecolor('blue')        
    xy += inc
    fxy = np.abs((xy+2)%4-2)
    pxy.set_data(*fxy)
    return pxy,

anim = FuncAnimation(f, update, frames=np.arange(1200) / 10,
                     init_func=init, blit=True)

pylab.show()

Solution 3:

So a couple things to keep in mind:

1. You need a friction component or else the particle will keep moving forever (conservation of energy). Friction happens as a function of velocity in this case, and friction will also happen on bounces.

2. If it's just a single particle, you can calculate it by defining bounding boxes like x is between 0 and 5, y is between 0 and 3 for example. Then you can calculate the intercept with a wall by plugging in the x = 5 value and then solving for y in the equation for the line.

For one particle, you don't have to do it parametrically with increments of t_0, you can calculate the intercepts and basically zoom it on over there. For multiple, you would have to calculate inter-molecular diffusion and collision forces...which is a much harder problem that should be done parametrically.

You would have to calculate collisions, which is when the centres of two molecules are 2*radius away from each other, then do a collision that conserves momentum.