# Simulating Brownian motion for N particles

The following is the web exercise 3.2.27. from the book Computer Science An Interdisciplinary Approach by Sedgewick & Wayne:

Write a data type that represents a particle undergoing a Brownian motion in two dimensions. Create a client program that takes a command line integer N, creates N particles at the origin, and simulates a Brownian motion for the N particles.

One thing to note: For the initial position of particles I chose random points throughout the blank page instead of the origin.

Here is my data-type implementation for Brownian particles:

public class BrownianParticle {
private double xPosition;
private double yPosition;

public BrownianParticle(double xPosition, double yPosition, double particleRadius) {
this.xPosition = xPosition;
this.yPosition = yPosition;
}

public double getXPosition() {
return xPosition;
}

public double getYPosition() {
return yPosition;
}

}

public BrownianParticle updatePosition(double xIncrement, double yIncrement) {
double random = Math.random();
if      (random < 1.0/8.0) {xPosition -= xIncrement; yPosition += yIncrement;}
else if (random < 2.0/8.0) {yPosition += yIncrement;}
else if (random < 3.0/8.0) {xPosition += xIncrement; yPosition += yIncrement;}
else if (random < 4.0/8.0) {xPosition += xIncrement;}
else if (random < 5.0/8.0) {xPosition += xIncrement; yPosition -= yIncrement;}
else if (random < 6.0/8.0) {yPosition -= yIncrement;}
else if (random < 7.0/8.0) {xPosition -= xIncrement; yPosition -= yIncrement;}
else if (random < 8.0/8.0) {xPosition -= xIncrement;}
}
}


Here is my data-type implementation for drawing Brownian particles:

import java.awt.Color;

public class BrownianParticleDraw {
private final BrownianParticle particle;

public BrownianParticleDraw(BrownianParticle particle) {
this.particle = particle;
}

public void draw() {
StdDraw.setPenColor(StdDraw.GRAY);
StdDraw.filledCircle(particle.getXPosition(),
particle.getYPosition(),
}

public void drawWithSpotlight() {
Color lightYellow = new Color(255,255,51);
StdDraw.setPenColor(lightYellow);
StdDraw.filledCircle(particle.getXPosition(),
particle.getYPosition(),
StdDraw.setPenColor(StdDraw.GRAY);
}
}


Here is my test client for simulating Brownian motion:

public class BrownianMotion {

public static double pickRandomlyBetween(double a, double b) {
return a + Math.random()*(b-a);
}

public static BrownianParticle[] InitializeRandomly(int numberOfParticles) {
BrownianParticle[] particles = new BrownianParticle[numberOfParticles];
for (int i = 0; i < numberOfParticles; i++) {
// StdDraw creates a 1-by-1 blank page (from 0 to 1) by default
double xPoistion = pickRandomlyBetween(0.05, 0.95);
double yPosition = pickRandomlyBetween(0.05, 0.95);
particles[i] = new BrownianParticle(xPoistion, yPosition, particleRadius);
}
return particles;
}

public static void animate(int trials, BrownianParticle[] particles) {
int numberOfParticles = particles.length;
StdDraw.enableDoubleBuffering();

for (int t = 1; t <= trials; t++) {
BrownianParticleDraw specificParticleDrawn = new BrownianParticleDraw(particles);

StdDraw.clear();

specificParticleDrawn.drawWithSpotlight();
specificParticleDrawn.draw();

for (int i = 1; i < numberOfParticles; i++) {
BrownianParticleDraw particleDrawn = new BrownianParticleDraw(particles[i]);
particleDrawn.draw();
}
StdDraw.pause(100);
StdDraw.show();
}
}

public static void main(String[] args) {
int numberOfParticles = Integer.parseInt(args);
int trials = Integer.parseInt(args);
BrownianParticle[] particles = InitializeRandomly(numberOfParticles);

animate(trials, particles);
}
}


StdDraw is a simple API written by the authors of the book. I checked my program and it works. Here are 2 different instances of it.

Input: 2 50

Output: Input: 200 50

Output: Is there any way that I can improve my programs (especially their performance)?

• Do particles pass through each other? – Ross Presser Sep 29 '20 at 21:06
• Unfortunately yes. There are no elastic collisions. – Khashayar Baghizadeh Oct 2 '20 at 10:14

Summary: You need to decide whether you are doing a random walk on a lattice (with int) or a random walk on a continuum (with double).

I'm going to discuss the physics and scientific computing aspects.

## Brownian motion on a square lattice

Contrary to what was said in another answer, the physics is not necessarily wrong. First, you might be simulating diffusion of molecules on top of a crystal plane with a square lattice (for example, the surface of gold called Au(100)). For this, a square grid would be completely appropriate. Furthermore, even in the absence of a physical lattice in what you are studying, your position distributions are isotropic in the limit of many steps. Your eventual xPosition is the sum of many random deviates from the set {-xIncrement, 0, +xIncrement} and yPosition is the sum of random deviates from the set {-yIncrement, 0, +yIncrement}. UpdatePosition() applies the displacements in x and the displacements in y in a way that they are independent. By the central limit theorem, xPosition and yPosition approach normal distributions, so that the position distribution after a large number of steps is:

$$P(x,y) \sim \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left( -\frac{x^2+y^2}{2\sigma^2} \right).$$

This has a circularly symmetric form: that is, it can be expressed as a function of distance from the initial point $$\r=\sqrt{x^2+y^2}\$$ only.

Here's an example of the approach to circular symmetry in Wolfram's New Kind of Science. I'm not sure how trustworthy this book is on everything, but I think he's probably right here. You are using a fixed increment, which means your particles stay on a fixed grid (lattice). Studying molecules on a lattice can be a useful approximation to allow efficient simulation. However, if you are doing a random walk on a lattice, xPosition and yPosition should be int, since all of the advantages of doing simulation on a lattice are best implemented with integers, which are discrete like the lattice. If xPosition and yPosition are supposed to be on a grid, you don't want to accumulate small floating-point errors that make all the particles slightly off the lattice. In the future, you may want to make it so two particles can't occupy the same positions. Checking whether a lattice position is occupied is easy with int positions. With double, checking whether two particles occupy the same lattice point could be error-prone due to the limited precision of floating-point representations.

## Brownian motion on a continuum

If you want to do this off-lattice with floating-point numbers, then you should implement the physics differently. @vnp's suggestion of randomly selecting an angle from 0 to 2π and incrementing by xPosition += increment*cos(angle) or yPosition += increment*sin(angle) is isotropic by construction and you wouldn't have to run many steps to see circular symmetry. However, the magnitude of the steps is wrong and you have to run many steps to approach true Brownian motion. A better approach that gives you true Brownian motion even for a single step is the algorithm of Ermak and McCammon.

If we ignore their hydrodynamic interactions, position-dependent diffusion coefficients and forces, the algorithm simplifies to:

$$x_{i+1} = x_i + g_t \sqrt{2 D \Delta t},$$

where $$\D\$$ is the diffusion coefficient (in units of $$\\mathsf{distance}^2/\mathsf{time}\$$), $$\\Delta t\$$ is the time step, and $$\g_t\$$ is a random number from a normal distribution with a mean of 0 and a variance of 1. You can get such normal deviates from the Java Random class using the nextGaussian() function or by generating them from uniformly distributed random numbers using a method like the Box–Mueller transform. The Ermak–McCammon algorithm is exact for your system, although it becomes approximate if you have forces between the particles. In the continuum, you can easily add interparticle forces and/or external forces between the particles:

$$x_{i+1} = x_i + \frac{D}{kT}F^{(x)}_i \Delta t + g^{(x)}_t \sqrt{2 D \Delta t},$$ $$y_{i+1} = y_i + \frac{D}{kT}F^{(y)}_i \Delta t + g^{(y)}_t \sqrt{2 D \Delta t},$$

where $$\F^{(x)}_i\$$ is the net force on particle $$\i\$$ in the $$\x\$$-direction. $$\\Delta t\$$ must be small enough that the force on each particle changes little between steps.

• Thank you very much for the great answer. – Khashayar Baghizadeh Oct 2 '20 at 10:16
• Physics is wrong.

The Brownian motion is isotropic. Your model is not. The steps in major directions are xIncrement and yIncrement ($$\I_x\$$ and $$\I_y\$$ respectively), but the steps along the diagonals are larger, $$\\sqrt{I_x^2 + I_y^2}\$$.

The tests hide it because of the random initial positions. Shall you follow instructions and set the particles at the origin, you would immediately see that the distribution of particles is not radially symmetrical.

Consider drawing a direction uniformly from a $$\[0, 2\pi)\$$ range, and a step from a geometric distribution. That will give much more realistic model.

• BrownianParticle.updatePosition returns new BrownianParticle. This is strange, especially considering that the return value is always ignored.

• I don't see how the performance can be improved (beside not creating new particles continuously).

• The fact that they have separate variable for x and y increments suggest that they want to be able to modify the code to allow for further anisotropy. Weird. – Acccumulation Sep 30 '20 at 3:53

## Naming

The names are not wrong, but I find them very verbose, which makes the code longer and harder to read and write.

private double xPosition;
private double yPosition;


It would be fine with x y and radius , and I see no risk of confusion as to the meaning of those names.

This spreads further through your code when using functions like .getParticleRadius() and getXPosition()

## UpdatePosition

 public BrownianParticle updatePosition(double xIncrement, double yIncrement) {
double random = Math.random();
if      (random < 1.0/8.0) {xPosition -= xIncrement; yPosition += yIncrement;}
else if (random < 2.0/8.0) {yPosition += yIncrement;}
else if (random < 3.0/8.0) {xPosition += xIncrement; yPosition += yIncrement;}
else if (random < 4.0/8.0) {xPosition += xIncrement;}
else if (random < 5.0/8.0) {xPosition += xIncrement; yPosition -= yIncrement;}
else if (random < 6.0/8.0) {yPosition -= yIncrement;}
else if (random < 7.0/8.0) {xPosition -= xIncrement; yPosition -= yIncrement;}
else if (random < 8.0/8.0) {xPosition -= xIncrement;}
}


This is wrong. This method should be public void and it should not return anything at the end. It changes the position of an existing particle and there is no need to create a new particle, like you do at the end. That is a waste of memory. This is confirmed by looking at where you call this function, and you are not assigning the returned new BrownianParticle to anything, so it just takes up some CPU and memory when created and then later it will be garbage collected and removed again.

## Constants

public void drawWithSpotlight() {
Color lightYellow = new Color(255,255,51);
StdDraw.setPenColor(lightYellow);
StdDraw.filledCircle(particle.getXPosition(),
particle.getYPosition(),
StdDraw.setPenColor(StdDraw.GRAY);
}


This creation of a new Color is wrong. Just like your previous chaos drawing program, you should define and use constants ( final Color ... ) for anything that you are re-using. If you don't use constants, this code will create a new color object that takes up memory every time it iterates, which is many times per seconds.

You need to understand what new means and be careful about using it in wasteful ways. It creates new objects and uses memory every time it is called.

• This is a nice example of where getters and setters make readability worse, because at best they litter the source code with useless pairs of (). If you must use them (e.g. for ideological reasons, not because they are necessary) x() and y() are better names than getXPosition() and getYPosition(). – alephzero Oct 1 '20 at 15:25

First of all, a general rule is that if you have more than a few else-ifs, you should start looking for a different control structure. For instance, you could replace your if-then-else-if block with xPosition += (round(3*random)%3-1)*xIncrement; yPosition += (round(9*random)%3-1)*yIncrement; But it's probably better to have two separate random numbers for x and y. And having different x and y increments is a bit weird. You end up plugging in the same numbers when you use the method, and Brownian motion generally has the same magnitude in both axes. And since the increment appears to depend on the radius, you might want to have the method simply access that attribute of the object, rather than having the external function access the attribute, then feed it back into the method.

Furthermore, Brownian motion is Gaussian, so you should generate random numbers from a normal distribution. This addresses the isotropy issue raised in another answer, as two-dimensional Gaussian motion can be decomposed into independent Gaussian motion along the two axes.

I also agree that the naming is overly wordy; having particle.getParticleRadius() is redundant. When you have a variable inside a class, you don't need to include class information in the variable name; that information is conveyed by the fact that it's part of the class.