Dhaval Adjodah final project
What did you do?
I wanted to simulate an ideal gas by creating a 2D box with bouncing balls undergoing elastic collisions to try to observe a Boltzmann distribution. The data obtained from the simulation would then be compared to an analytically derived expression of the Boltzmann distribution, and verified using non-linear fitting.
How did you do it?
Using Java, the balls were created as objects each with random x and y positions in a square box. the speed of each ball is initially chosen to be a set value, with the x and y velocity components of this speed being assigned from a using a randomly (uniform distribution) generated angle.
The balls are then allowed to bounce elastically with each other and with the walls. Since an elastic collision is defined by the conservation of momentum and energy, it was important to get the collision geometry right. To do so, at the time of collision, the velocity vector of each ball was resolved along the ball center-to-center (normal) axis and the tangent direction and a 1D elastic collision was performed along the normal direction. After normal velocity vectors were updated, the velocity vector of each ball was projected back to the x-y directions. Insights about how to achieve this was drawn from the mathematics and code of http://www.vobarian.com/bouncescope/index.html
Collisions were detected by checking each ball for every other ball to see if the distance between the centers of the balls were less than their diameters. This algorithm is O[n^2] but performed satisfactorily (less than 10 seconds to reach equilibrium) with 2500 balls of diameter 10 pixels in a box of size 4000x4000 pixels. Note that this leads to an ball to total volume occupation ratio of about 1% which is the within the ideal gas limit.
When equilibrium is reached (determined by looking at the GUI display for balls for maximum randomness, and checking the velocity histograms for convergence of mean and variance), the x and y velocity of each ball, and its speed are saved to file using Java's IO libraries.
Using python, the data is then loaded using numpy into an array and converted to a histogram using matplotlib. The x and y velocity correctly form a Gaussian distribution at equilibrium (even though they start out as a uniform distribution) and the speed seems to have the positive skew of a Boltzmann distribution.
Since the variance of either the x or y direction (they are very close), sqrt(kT/m) for a given system of gas particles is the only parameter to a Holzman distribution of speeds, the generic Boltzmann distribution function using the average of the x and y variance was plotted onto the simulation data. It was clear that the theoretical function did not fit the data.
After a lot of brainstorming, and trials, it seems that the problem came from the fact that the theoretical version of the Boltzmann distribution function assumes that each particle has three degrees of freedom, whereas my simulation assumes that a particle has only x and y kinetic energy.
Using the law of equipartition, and using the fact that the x and y velocities are correctly Gaussian-distributed, a 2D Boltzmann distribution function was analytically derived (pdf in repository). The analytical curve was plotted next to the simulation data and seemed to fit perfectly. Insights for the derivation were found at http://galileo.phys.virginia.edu/classes/252/kinetic_theory.html
To verify the fit and the validity of the correspondence between the derived and simulation data, a non-linear fit was sought for the constants and exponents of the general form of a Boltzmann distribution function. They matched the analytical solution correctly within 7%, 1% and 2% error for the three parameters - seed attached document)
How did you select the approach?
Given my previous experience with object-oriented programming, my natural first instinct was to model the system using bouncing ball objects. The collision detection mechanism was supposed to start out being the simplest (and least efficient one O[n^2]) and then upgraded to one that exploits the power of binary tree partitioning, or quad-trees, but since the simulation reached equilibrium with more than enough balls under 10 seconds, this was not updated.
I wanted to do the graphing and possibly data analysis of the data using Java libraries with the aim that one could observe a Gaussian distribution form in real-time as the balls bounce to equilibrium. After many days of frustration trying to get the JFreeChart libraries to update in real-time, I gave up and switched up writing the ball data to file and reading them using python numpy, matplotlib and scipy.
matplolib, numpy and scipy were like a breeze to load and plot the data.
I tried to learn how to write my own Levenberg-Marquardt routine, but eventually decided (for the same of time, and learning one step at a time - not being used to python's scientific environment) to use scipy.optimize's leastsq routine which performs least square minimization given and error function over a vector of parameters.
How does it relate to what's been done before?
Entire ideal gas simulations exists written as Cellular Automata lattice gas models. The conceptual difficulty of CA made it not the preferable choice for me to use this as my model.
Collision detection schemes have been used that employ space-partitioning and binary/quad-tree algorithms leading, for example, to O[n log n] time, but again, performance was not an issue.
Various derivations of Boltzmann gasses exists using different routes in 3D. The extension to 2D - once the right transformations were identified - was relatively straightforward.
How did you validate the model?
The model was validated through three means:
The verified zero mean and Gaussian of the equilibrium x and y velocities.
The fact that, if started closed enough, the parameters of the non-linear fit to simulated data match those of the analytical derivation.
How did you evaluate the results?
Using the metrics above, as long as errors were within a margin of 5-10% for Gaussian and Boltzmann distribution, the correspondence between the analytically derived distribution and the simulation data was verified.
What did you learn?
I learned a great deal about how to combine the Model-View-Controller design with file-writing, vector algebra and data visualization extensions.
I learned how to used the python scientific analysis environment.
I learned how to do non-linear fitting and the beauty and grace of python scripting.
I learned how to work with large data sets and match data to analytically-derived expressions.
What were the problems?
Getting the collision geometry and momentum-energy conservation to work took a lot of nit-picky debugging. Eventually, very specific test cases (for example along certain axes, etc) had to be chosen and run to solve problems.
Getting the Java real-time graphing proved to be quite a challenge. It is fairly easy to write my own histogram data and visual rendering routine, but it would take a long time and take away from the actually simulation data analysis. Eventually, the real-time rendering extension was dropped and writing to file from Java and reading from python was chosen. However, it is possible to have a python script watching a file to see any change to the data and having it plot every time Java writes to the file - this would have been to time consuming and probably would take too much processor time anyway.
Writing my own Levenberg-Marquardt routine proved to be a challenge given my poor understanding of it - hence the scipy least square method was used.
Semester Assignments:I wasn't able to complete all the coding assigments, but here are all the pdfs of the math parts I could find!