Kinetic Theory of Gases#
In your lab you made a simulation of a gas of hard discs bouncing around in a billiard. From an algorithm standpoint, you learned about event driven simulations, and some of the tricks of making these efficient.
Today we will discuss some of the physics around this system. It is one of the iconic many-body systems, which is used to think about many physical processes.
Boltzmann distribution#
The first thing that we saw in this simulation is that after a few collsions, the state of the system was essentially random. It is chaotic, and small changes in initial conditions lead to completely different configurations later. Thus we are naturally led to a statistical description of the system. For example, if you look at a random frame of the simulation, and look at the speed of one of the particles, what distribution do you expect the speed to be drawn from?
The generic argument here starts with the premise of ergodicity. We posit that things are so random that any state of the system, consistent with the conservation laws, is equally likely. Thus if we had 1 particle, it will have a random velocity whose magnitude is fixed, but whose direction is uniformly distributed. With 2 particles, \(v_{x1},v_{x2},v_{y1},v_{y2}\) are uniformly distributed with the constraint that \(v_{x1}^2+v_{x2}^2+v_{y1}^2+v_{y2}^2\) is fixed. In the limit of a large number of particles this leads to a universal result.
Working with continuous probability distributions is a bit of a pain, so lets imagine that we discretize velocity-space. We define \(\Omega(N,E)\) be the number of ways of choosing the velocities of \(N\) particles, such that the total energy is \(E\).
Lets now consider a system of \(N+1\) particles. We want to know what the probability that the extra particles has velocity \(\vec v\). To find that, we just count how many ways of arranging all of the other particles:
here \(\epsilon_v=mv^2/2\) is the energy of the particle we are caring about.
The next stage of the argument is to note that \(\Omega\) is exponentially big in \(N\). Thus it is natural to write
where \(S\) should be proportional to \(N\). We call \(S\) the entropy. The factor of \(k_B\) is just convention – it gives the units of \(S\). The only physical quantity is \(k_B S\). I often take \(k_B=1\).
Next we note that generically \(\epsilon_v\ll E\), so we can Taylor expand,
The derivative \(\partial S/\partial E\) is some function of \(N\) and \(E\). We want to give this function a name. By convention, we define
where the function \(T(N,E)\) is known as the temperature.
Putting this together, we have
which is known as the Boltzmann distribution. To find the distribution for the speed, we note that
As a final step we normalize the distribution by noting that the total probability is equal to 1,
and hence
Note, this result depends on dimension – so the 3D results look slightly different.
using Plots,Logging,ProgressMeter,LinearAlgebra,LaTeXStrings
abstract type Event end
struct NullEvent <: Event
t ::Float64
end
NE=NullEvent(Inf)
mutable struct Ball
x ::Float64
y ::Float64
vx ::Float64
vy ::Float64
radius ::Float64
nextevent ::Event
end
function set_xv!(ball ::Ball,xy,vxy)
ball.x,ball.y=xy
ball.vx,ball.vy=vxy
return ball
end
function set_r!(ball ::Ball,x,y)
ball.x=x
ball.y=y
return ball
end
function set_r!(ball ::Ball,r)
ball.x,ball.y=r
return ball
end
function set_v!(ball ::Ball,v)
ball.vx,ball.vy=v
return ball
end
function set_v!(ball ::Ball,vx,vy)
ball.vx=vx
ball.vy=vy
return ball
end
function get_r(ball ::Ball)
return [ball.x,ball.y]
end
function get_v(ball ::Ball)
return [ball.vx,ball.vy]
end
Base.show(io::IO,b::Ball)=print(io,
"Ball << r="*repr((b.x,b.y))*
" v="*repr((b.vx,b.vy))*">>)")
Ball(;x,y,vx,vy,radius,nextevent=NE)=Ball(x,y,vx,vy,radius,nextevent) # keyword constructor
function Ball(r,v,radius,nextevent=NE) # constructor using containers
x,y=r
vx,vy=v
Ball(x,y,vx,vy,radius,nextevent)
end
Ball(;r,v,radius,nextevent=NE)=Ball(r,v,radius,nextevent) # keyword constructor using containers
abstract type Wall end
struct Hwall <: Wall
y ::Float64
end
struct Vwall <: Wall
x ::Float64
end
function impact_time(w::Hwall,b::Ball,t)
dy=w.y-b.y
return (dy-sign(b.vy)*b.radius)/b.vy+t
end
function impact_time(w::Vwall,b::Ball,t)
dx=w.x-b.x
return (dx-sign(b.vx)*b.radius)/b.vx+t
end
impact_time(b::Ball,w::Wall,t)=impact_time(w,b,t)
function circ_time(r0,v,r)
# Calculate the various terms in our expression
r0dotv=r0⋅v
if r0dotv>0 #moving away from one another
return Inf
end
vsq=v⋅v
r0sq=r0⋅r0
if r0sq<=r^2 #balls are interpenetrating
return Inf
end
discriminant=r0dotv^2+(r^2-r0sq)*vsq
if discriminant<0 # balls will never hit one-another
return Inf
end
times=[(-r0dotv-sqrt(discriminant))/(vsq),(-r0dotv+sqrt(discriminant))/(vsq)]
t1,t2=minmax(times...)
if t2<=0
return Inf
end
if t1>0
return t1
end
return t2
end
impact_time(b1::Ball,b2::Ball,t)=
t+circ_time(get_r(b1)-get_r(b2),
get_v(b1)-get_v(b2),b1.radius+b2.radius)
struct Collision <: Event
object1 :: Union{Ball,Wall}
object2 :: Union{Ball,Wall}
t ::Float64
end
mutable struct Billiard
walls ::Array{Wall}
balls ::Array{Ball}
t ::Float64
nextevent ::Event
boundingbox ::Array{Float64,1}# x1,y1,x2,y2 = corners of a rectangle
end
Billiard(;walls,balls,t=0.,nextevent=NE,boundingbox)=
Billiard(walls,balls,t,nextevent,boundingbox)
function next_event(b::Ball,billiard::Billiard)
t=billiard.t
t_event=Inf
result=NE
for obj in union(billiard.balls,billiard.walls)
t_try=impact_time(b,obj,t)
if t_try<t_event && t_try>t
t_event=t_try
result=Collision(b,obj,t_try)
end
end
return result
end
function drawcircles!(centers,radii;opts...)
c=[cos(θ) for θ in 0:2*pi/100:2pi]
s=[sin(θ) for θ in 0:2*pi/100:2pi]
makecircle(center,radius)=Shape(radius*c.+center[1],radius*s.+center[2])
circles=[makecircle(c,r) for (c,r) in zip(centers,radii)]
plot!(circles;aspect_ratio=:equal,label="",opts...)
end
function drawcircles!(plt::Plots.Plot,centers,radii;opts...)
c=[cos(θ) for θ in 0:2*pi/100:2pi]
s=[sin(θ) for θ in 0:2*pi/100:2pi]
makecircle(center,radius)=Shape(radius*c.+center[1],radius*s.+center[2])
circles=[makecircle(c,r) for (c,r) in zip(centers,radii)]
plot!(plt,circles;aspect_ratio=:equal,label="",opts...)
end
function drawcircles(centers,radii;opts...)
plt=plot()
drawcircles!(plt,centers,radii;opts...)
end
# Master function for visualizing
# -- creates new plot object, and then calls visualize!
function visualize(opts...;namedopts...)
viz=plot()
visualize!(viz,opts...;namedopts...)
return viz
end
function visualize!(viz::Plots.Plot,w::Hwall,bbox;opts...)
x1,y1,x2,y2=bbox
y=w.y
plot!(viz,[x1,x2],[y,y];opts...)
end
function visualize!(viz::Plots.Plot,w::Vwall,bbox;opts...)
x1,y1,x2,y2=bbox
x=w.x
plot!(viz,[x,x],[y1,y2];opts...)
end
function visualize!(viz::Plots.Plot,b::Billiard;wallopts=[],ballopts=[])
walls=b.walls
balls=b.balls
bbox=b.boundingbox
for w in walls
visualize!(viz,w,bbox;label="",wallopts...)
end
centers=[(b.x,b.y) for b in balls]
radii=[b.radius for b in balls]
drawcircles!(viz,centers,radii;ballopts...)
return viz
end
function setevents!(billiard::Billiard)
t_global=Inf
for b in billiard.balls
next=next_event(b,billiard)
b.nextevent=next
if next.t<t_global
t_global=next.t
billiard.nextevent=next
end
end
return billiard
end
reflect(vec,normal)=vec-2*(vec⋅normal)*normal/(normal⋅normal)
normal(w::Hwall,r)=[0,1]
normal(w::Vwall,r)=[1,0]
function newvelocity(w::Wall,b::Ball)
r=get_r(b)
n=normal(w,r)
v=get_v(b)
return reflect(v,n)
end
function set_newvelocity!(w::Wall,b::Ball)
v=newvelocity(w,b)
set_v!(b,v)
return b
end
set_newvelocity!(b::Ball,w::Wall)=set_newvelocity!(w,b)
function set_newvelocity!(b1::Ball,b2::Ball)
r1=get_r(b1)
r2=get_r(b2)
v1=get_v(b1)
v2=get_v(b2)
norm=r1-r2
vcm=(v1+v2)/2
r12=r1-r2
newv1=reflect(v1-vcm,r12)+vcm
newv2=reflect(v2-vcm,r12)+vcm
set_v!(b1,newv1)
set_v!(b2,newv2)
return (newv1,newv2)
end
function time_evolve!(b::Ball,dt)
r=get_r(b)
v=get_v(b)
newr=r+v*dt
set_r!(b,newr)
end
function update!(billiard::Billiard,event::Event)
dt=event.t-billiard.t
for b in billiard.balls
time_evolve!(b,dt)
end
set_newvelocity!(event.object1,event.object2)
billiard.t=event.t
return billiard
end
struct State
xvalues ::Array{Float64,1}
yvalues ::Array{Float64,1}
vxvalues ::Array{Float64,1}
vyvalues ::Array{Float64,1}
radii ::Array{Float64,1}
t ::Float64
end
struct Simulationdata
walls :: Array{Wall}
statelist :: Array{State}
boundingbox ::Array{Float64,1}
end
Base.length(s::Simulationdata)=length(s.statelist)
function State(b::Billiard)
xvalues = [ball.x for ball in b.balls]
yvalues = [ball.y for ball in b.balls]
vxvalues = [ball.vx for ball in b.balls]
vyvalues = [ball.vy for ball in b.balls]
radii = [ball.radius for ball in b.balls]
t=b.t
State(xvalues,yvalues,vxvalues,vyvalues,radii,t)
end
function randomgrid(nx,ny,r,Lx,Ly)
xstep=Lx/nx
ystep=Ly/ny
function randomball(x,y)
θ=2*pi*rand()
Ball((x,y),(cos(θ),sin(θ)),r)
end
Billiard([Hwall(Ly/2),Hwall(-Ly/2),Vwall(Lx/2),Vwall(-Lx/2)],
[randomball(x,y)
for x in -Lx/2+xstep/2:xstep:Lx/2-xstep/2
for y in -Ly/2+ystep/2:ystep:Ly/2-ystep/2],
0., NE,[-Lx/2,-Ly/2,Lx/2,Ly/2])
end
randomgrid(;nx,ny,r,Lx,Ly)=randomgrid(nx,ny,r,Lx,Ly)
function setevents!(billiard::Billiard,lastevent::Collision)
# The first thing we do is make a list of all the `Ball`
# objects in `lastevent`. We store this in `cballs
cballs=Ball[]
if typeof(lastevent.object1)==Ball
push!(cballs,lastevent.object1)
end
if typeof(lastevent.object2)==Ball
push!(cballs,lastevent.object2)
end
# Next unpack the data from `billiard`
balls=billiard.balls
walls=billiard.walls
t=billiard.t
# set null events for unknown collisions
billiard.nextevent=NE
for cb in cballs
cb.nextevent=NE
end
# Loop over balls
for b in balls
for cb in cballs
ct=impact_time(b,cb,t)
if ct<=t ;continue; end
if ct<b.nextevent.t
b.nextevent=Collision(b,cb,ct)
end
if ct<cb.nextevent.t
cb.nextevent=Collision(b,cb,ct)
end
if ct<billiard.nextevent.t
billiard.nextevent=Collision(b,cb,ct)
end
if b.nextevent.t<billiard.nextevent.t
billiard.nextevent=b.nextevent
end
end
end
# loop over walls
for w in walls
for cb in cballs
ct=impact_time(w,cb,t)
if ct<=t
continue
end
if ct<cb.nextevent.t
cb.nextevent=Collision(w,cb,ct)
end
if ct<billiard.nextevent.t
billiard.nextevent=Collision(w,cb,ct)
end
end
end
return billiard
end
function simulate(billiard::Billiard,n)
result=Array{State}(undef,n+1)
result[1]=State(billiard)
setevents!(billiard)
for j in 2:n+1
setevents!(billiard,billiard.nextevent)
update!(billiard,billiard.nextevent)
result[j]=State(billiard)
end
return Simulationdata(billiard.walls,result,billiard.boundingbox)
end
"""
speedlist(sim::Simulationdata,frame::Integer)
Makes a list of the speeds of all of the particles in frame number
`frame` of the simulation data. You can also call it with a list of
frames and it will concatentate those lists.
For example, you can call with
speedlist(s,[1,3,5])
and get the speeds in frames 1, 3, 5. We can even use
speedlist(s,1::10::100)
which will do every 10'th frame between 1 and 100.
"""
speedlist(sim::Simulationdata,frame::Integer) =
sqrt.(sim.statelist[frame].vxvalues.^2
+sim.statelist[frame].vyvalues.^2)
speedlist(sim::Simulationdata,frames)=
vcat([speedlist(sim,j) for j in frames]...)
# Visualization
function drawcircles!(plt::Plots.Plot,xvalues,yvalues,radii;opts...)
c=[cos(θ) for θ in 0:2*pi/100:2pi]
s=[sin(θ) for θ in 0:2*pi/100:2pi]
makecircle(x,y,radius)=Shape(radius*c.+x,radius*s.+y)
circles=[makecircle(x,y,r) for (x,y,r) in zip(xvalues,yvalues,radii)]
plot!(plt,circles;aspect_ratio=:equal,label="",opts...)
end
function visualize!(plt::Plots.Plot,s::State;opts...)
drawcircles!(plt,s.xvalues,s.yvalues,s.radii;opts...)
end
function drawframe(s::Simulationdata,j)
walls=s.walls
bbox=s.boundingbox
plt=plot()
for w in walls
visualize!(plt,w,bbox;label="")
end
visualize!(plt,s.statelist[j])
end
#
#
# function eventanimate(s::Simulationdata,skip=1)
# # print out a message to let the user know things are starting
# @info "Animating Simulation"
# len=length(s)
# progress_meter=Progress(len,desc="Generating frames: ")
# anim = @animate for i in 1:skip:len
# ProgressMeter.update!(progress_meter,i)
# drawframe(s,i)
# end
# Plots.gif(anim, fps = 12)
# end
drawframe (generic function with 1 method)
While the Plots.jl visualizations work, they are kind of slow. Here is a quick GLMakie script to animate the Simulationdata
using GLMakie
function drawcircles!(ax::Axis,centers,radii;opts...)
circles=[Circle(Point2f(c),r) for (c,r) in zip(centers,radii)]
poly!(ax,circles;opts...)
end
function drawcircles!(plt::Axis,xvalues,yvalues,radii;opts...)
circles=[Circle(Point2f(x,y),r) for (x,y,r) in zip(xvalues,yvalues,radii)]
poly!(plt,circles;opts...)
end
function newplot()
f=Figure()
ax=Axis(f[1,1],aspect=1.0)
display(f)
return f,ax
end
function visualize!(viz::Axis,w::Hwall,bbox;opts...)
x1,y1,x2,y2=bbox
y=w.y
lines!(viz,[x1,x2],[y,y];opts...)
end
function visualize!(viz::Axis,w::Vwall,bbox;opts...)
x1,y1,x2,y2=bbox
x=w.x
lines!(viz,[x,x],[y1,y2];opts...)
end
function visualize!(viz::Axis,b::Billiard;wallopts=[],ballopts=[])
walls=b.walls
balls=b.balls
bbox=b.boundingbox
for w in walls
visualize!(viz,w,bbox;label="",wallopts...)
end
centers=[Point2f(b.x,b.y) for b in balls]
radii=[b.radius for b in balls]
drawcircles!(viz,centers,radii;ballopts...)
return viz
end
function setcircles!(circles,centers,radii)
for j in eachindex(circles)
circles[j]=Circle(Point2f(centers[j]),radii[j])
end
end
function setcircles!(circles,xs,ys,radii)
for j in eachindex(circles)
circles[j]=Circle(Point2f(xs[j],ys[j]),radii[j])
end
end
setcircles(centers,radii)=[Circle(Point2f(c),r) for (c,r) in zip(centers,radii)]
setcircles(xs,ys,radii)=[Circle(Point2f(x,y),r) for (x,y,r) in zip(xs,ys,radii)]
function eventanimate(s::Simulationdata,timestep=0.001)
# Create Plot
f,ax=newplot()
Makie.deactivate_interaction!(ax, :rectanglezoom)
# Draw Walls
for w in s.walls
visualize!(ax,w,s.boundingbox)
end
#Draw Circles
state=s.statelist[1]
circles=setcircles(state.xvalues,state.yvalues,state.radii)
obscircles=Observable(circles)
poly!(ax,obscircles)
# Create Slider
controls=Slider(f[2,1],range=1:1:length(s.statelist),startvalue=1)
on(controls.value) do j
state=s.statelist[j]
setcircles!(circles,state.xvalues,state.yvalues,state.radii)
notify(obscircles)
end
# Create Run button
isrunning = Observable(false)
label = map(cond -> cond ? "Stop" : "Run", isrunning)
run = Button(f[3,1]; label = label, tellwidth = false)
on(run.clicks) do clicks
isrunning[] = !isrunning[]
end
isrunning_notifier = Condition()
on(cond -> cond && notify(isrunning_notifier), isrunning)
errormonitor(@async while true
if isrunning[]
j=controls.value[]
maxj=controls.range[][end]
while controls.value[]==maxj
isrunning[]=false
wait(isrunning_notifier)
end
isopen(f.scene) || break # stops if window is closed
set_close_to!(controls,j+1)
sleep(timestep)
else
wait(isrunning_notifier)
end
end)
end
# Make Plots the default package for some of the
# typical commands
scatter=Plots.scatter
scatter! =Plots.scatter!
heatmap=Plots.heatmap
heatmap! =Plots.heatmap!
heatmap! (generic function with 1 method)
Running the simulations#
Lets start with a random grid of 100 particles.
It is useful to name our variables in such a way that we can figure out what parameters they correspond to.
g_N100_r01=randomgrid(nx=10,ny=10,r=0.01,Lx=1,Ly=1)
visualize(g_N100_r01)
WARNING: both GLMakie and Plots export "plot"; uses of it in module Main must be qualified
WARNING: both GLMakie and Plots export "plot!"; uses of it in module Main must be qualified
UndefVarError: `plot` not defined
Stacktrace:
[1] visualize(opts::Billiard; namedopts::@Kwargs{})
@ Main ./In[1]:177
[2] visualize(opts::Billiard)
@ Main ./In[1]:176
[3] top-level scope
@ In[3]:2
g_N100_r01.balls
100-element Vector{Ball}:
Ball << r=(0.09736651127562534, -0.38473254411470137) v=(0.3771157618524579, -0.06557828744267591)>>)
Ball << r=(0.15763542829785052, -0.24765240087742543) v=(0.6712329981557774, -0.5018738247904133)>>)
Ball << r=(-0.13170764607361318, 0.08575819109140606) v=(-0.24551939216299357, -0.5424908591919243)>>)
Ball << r=(0.030692047330353, -0.10961413115954946) v=(-0.4204863282573531, 0.7553085474861749)>>)
Ball << r=(-0.44411361627827484, -0.14760786041012697) v=(0.43805588275968055, -0.27320369336858596)>>)
Ball << r=(-0.44707909016646735, 0.4855816259466236) v=(-0.17391893757044552, -0.2807754919884099)>>)
Ball << r=(0.3326702785584588, -0.35915448670630457) v=(0.6647476563974672, 0.08685753834903664)>>)
Ball << r=(-0.17521396344517276, 0.2893111284110229) v=(1.3091055646501393, 0.2601140896915255)>>)
Ball << r=(-0.2244984731869328, 0.2858647103061798) v=(0.4465770725243178, -1.0530926089699857)>>)
Ball << r=(-0.11093775721628042, -0.26398585648000555) v=(-0.5531097652470542, -0.1248026808184827)>>)
Ball << r=(-0.3455295146581722, 0.2161855178320665) v=(-0.3196864933640767, 0.5469643364150573)>>)
Ball << r=(-0.17782558459748599, 0.23593407363440624) v=(-0.2339731744552337, -0.4947666808914324)>>)
Ball << r=(-0.27018501321133087, -0.40299780832534016) v=(0.3310614039588924, 0.4255492474160527)>>)
⋮
Ball << r=(0.01494445249540023, -0.44618028756190686) v=(-0.33108658001019664, -1.9881868387116122)>>)
Ball << r=(-0.07984232652391424, 0.010462336861435608) v=(-1.1672863022848354, 0.3771511027641144)>>)
Ball << r=(-0.4265275218885858, -0.27782105148920994) v=(0.8607228306379522, 1.0485527650846782)>>)
Ball << r=(-0.4817968245092172, -0.2583997756074282) v=(1.2230197525032247, 0.16051486005252308)>>)
Ball << r=(-0.22795162474559003, 0.010807485629092575) v=(-0.35176129391356226, -0.1544870692453283)>>)
Ball << r=(0.03500815940466245, 0.2906454254364328) v=(1.3260277049681104, -0.5206319972411695)>>)
Ball << r=(0.38301652076152576, -0.3931406965214097) v=(0.24831143397672029, 1.4720115684077721)>>)
Ball << r=(0.386128447155886, 0.06683307920977116) v=(-0.14296745753487428, 0.25385362527624755)>>)
Ball << r=(0.15079679051890998, -0.40684133251586724) v=(0.18813335278761312, -1.5553861539480542)>>)
Ball << r=(0.06276931628996894, 0.2715326821863609) v=(-0.016865197019723777, -0.6987227879994578)>>)
Ball << r=(0.28323186917054616, -0.29531733487316497) v=(-0.8864065141703446, 1.0382373299481542)>>)
Ball << r=(-0.1341491535039915, -0.03525747031662463) v=(-0.12481322784847076, -0.24113095209421737)>>)
Now lets do a “equilibration” run, where we evolve for enough time for the initial conditions to be forgotten. Later we will figure out how to determine when things have equilibrated.
@time s_N100_r01=simulate(g_N100_r01,5000);
visualize(g_N100_r01)
0.630967 seconds (15.02 M allocations: 716.513 MiB, 13.57% gc time, 12.37% compilation time)
eventanimate(s_N100_r01)
Task (runnable) @0x000000036ca407e0
Now we can do a “production” run, where we extract the speeds of the particles after every 50 collisions
@time sfull_N100_r01=simulate(g_N100_r01,50000);
histogram(speedlist(sfull_N100_r01,1:50:50000),normalize=true,
label="100 particles, r=0.01")
plot!(v->2*v*exp(-v^2),0,3,label="Maxwell-Boltzmann")
4.338162 seconds (148.81 M allocations: 6.891 GiB, 23.43% gc time)
eventanimate(sfull_N100_r01,0.000001)
Task (runnable) @0x000000035430f9e0
Correlation time, and correlation functions#
To further understand the structures in this gas, it is useful to think about correlations.
We say that random variables \(X\) and \(Y\) are uncorrelated if the value of \(X\) does not depend on the value of y. For example, here are some uncorrelated random variables
xvals=rand(1000)
yvals=rand(1000)
scatter(xvals,yvals,xlabel="x",ylabel="y",label="",title="Uncorrelated variables")
Here is another, less trivial, uncorrelated example
xvals2=xvals.*(1 .-xvals)
yvals2=exp.(-yvals)
scatter(xvals2,yvals2,
xlabel="x",ylabel="y",label="",title="Uncorrelated variables")
there the x-value and y-value are drawn from non-trivial distributions, but the value of x does not depend on the value of y.
Here is an example of correlated random variables
xvals3=xvals
yvals3=xvals-0.1*yvals
scatter(xvals3,yvals3,
xlabel="x",ylabel="y",label="",title="Correlated variables")
Here the values of \(y\) depend on the values of \(x\).
Here is another example:
xvals4=@. sqrt(xvals)*cos(2*pi*yvals)
yvals4=@. sqrt(xvals)*sin(2*pi*yvals)
scatter(xvals4,yvals4,
xlabel="x",ylabel="y",label="",title="Correlated variables (but nonlinear correlations)",aspectratio=1)
One common way to check to see if two random variables are correlated is to calculate the correlation function \(\langle XY\rangle-\langle X\rangle \langle Y \rangle\), also referred to as the covariance.
using Statistics
cfun(xvals,yvals)=mean(xvals.*yvals)-mean(xvals)*mean(yvals)
cfun (generic function with 1 method)
@show cfun(xvals,yvals)
@show cfun(xvals2,yvals2)
@show cfun(xvals3,yvals3);
cfun(xvals, yvals) = -0.000836323667062544
cfun(xvals2, yvals2) = 0.0003604819598582998
cfun(xvals3, yvals3) = 0.08439312831500959
The correlations in the first two are much smaller than those in the third. They are non-zero because we are working with a finite sample. Any particular sample will have some correlations, just by chance. Later we will work out how large we expect those to be. One of the primary uses of “statistics” is to determine if the random fluctuations are consistent with uncorrelated noise.
There is a small caution:
@show cfun(xvals4,yvals4);
cfun(xvals4, yvals4) = -0.0030791457090725547
Even though xvals4 and yvals4 are correlated – it is a nonlinear form of correlation, which is not diagnosed by the covariance.
Applying this idea to the problem at hand, we might ask how the velocites at two times are related to one-another by calculating \(C_{xx}(t,t^\prime)=\langle v_x(t) v_x(t^\prime)\rangle\), and similar expresions for other components. We don’t need to subtract off the mean, since we expect that \(\langle v_x(t)\rangle\)=0.
Because of reflection invarience, we expect \(C_{xy}=0\) and \(C_{xx}=C_{yy}\). Thus the natural quantity to calculate is \(C(t,t^\prime)=\langle \vec{v}(t) \cdot\vec{v}(t^\prime)\rangle=C_{xx}+C_{yy}=2 C_{xx}\). If we are in some sort of steady-state, then it makes sense to average over different times, writing
The most natural way to calculate this is to first make a time series of the velocities vlist=[v1,v2,v3,...vT] where vj=v(tj) is a vector of length \(2N\) containing the components of the velocity at each time.
Our estimate of \(C(\tau=0)\) is
Our estimate of \(C(\tau=\delta t)\) is
and so on.
The naive implementation would take a time which scales as \(O(N T^2)\). There are smart ways of doing the convolution which reduce this to \(O(N T\log(T))\). At some point we may talk about that algorithm. Julia has an implementation of the smart algorithm, which is in the DSP package, and called conv. There is also an implementation in the StatsBase package, called autocov, which I believe is just the naive algorithm (written in a pretty nice way).
Lets just write our own.
The first thing we need is the time-series data of the velocities. We could either modify our simulation code to extract that, or we can extract it from the SimulationData that the current program returns. Lets do the latter (as that was the whole reason for returning SimulationData.
getv(s::State)=vcat(s.vxvalues,s.vyvalues)
function velocitytimeseries(s::Simulationdata,dt)
# set up data structures
states=s.statelist
t=states[1].t
velocities=[getv(states[1])]
# step through collisions
for j in 2:length(states)
tj=states[j].t
# if collision is more than dt in the future
# then we want to sample the velocities,
# and increment t.
#
# We use a while loop, as we the next
# collision might be multiple dt's in the future
#
# We want to "eat up" all those timesteps
# before incrementing to the next collision
#
# The point is that the velocities are piece-wise constant
#
while tj>t+dt
t+=dt
prev=states[j-1]
push!(velocities,getv(prev))
end
end
return velocities
end
velocitytimeseries (generic function with 1 method)
timestep=0.01
vs1=velocitytimeseries(sfull_N100_r01,timestep)
scatter([v[1] for v in vs1[1:100]],xlabel="timestep",ylabel="vx",label="vx1")
scatter!([v[4] for v in vs1[1:100]],xlabel="timestep",ylabel="vx",label="vx4")
length(vs1)
8674
Now, by hand we can look at the correlations. We can just change the indices in this scatter plot to look at correlations in different time slices
scatter(vs1[1],vs1[1],
xlabel="vx,vy at first time",
ylabel="vx, vy at second time",label="")
Initially the velocities are perfectly correlated, but for large time separations they appear random.
function autocorrelation(velocities,largestsep)
timesteps=length(velocities)
numv=length(velocities[1])
[(velocities[1+j:end]⋅velocities[1:end-j])*2/(numv*(timesteps-j)) for j in 0:largestsep]
end
autocorrelation (generic function with 1 method)
times=range(start=0.,length=501,step=timestep)
scatter(times,autocorrelation(vs1,500),label="",
xlabel=L"$\tau$",ylabel=L"$\langle v(t+\tau)v(t)\rangle$")
That is kind of fun. At short times the velocities are clearly correlated: These correlations fall off on a timescale of order the collision rate. Then there is a period during which the velocities are anti-correlated. This makes sense, after you bounce off of another particle, you tend to be moving in the opposite direction from your start.
The typical time between collsions is
states1=sfull_N100_r01.statelist
finaltime=states1[end].t
initialtime=states1[1].t
numberofcollisions=length(states1)
numberofparticles=length(states1[1].xvalues)
meancollisiontime=(finaltime-initialtime)/(numberofcollisions/numberofparticles)
0.17346736804378635
This is sometimes referred to as the mean-free time. Lets scale our correlation function by the collision time, to give a scale to things
scatter(times./meancollisiontime,autocorrelation(vs1,500),label="",
xlabel=L"$\tau/\tau_c$",ylabel=L"$\langle v(t+\tau)v(t)\rangle$")
Finally, lets look at how the collision time changes with the radius of the spheres
function collisiontime(;
nx=10,ny=10,radius=0.01,Lx=1,Ly=1,
equilibrationcollisions=1000,
samplingcollisions=5000)
b=randomgrid(nx,ny,radius,Lx,Ly)
simulate(b,equilibrationcollisions)
sim=simulate(b,samplingcollisions)
states=sim.statelist
finaltime=states[end].t
initialtime=states[1].t
numberofcollisions=length(states)
numberofparticles=length(states[1].xvalues)
meancollisiontime=(finaltime-initialtime)/(numberofcollisions/numberofparticles)
end
collisiontime (generic function with 1 method)
@time collisiontime(radius=0.01)
0.497674 seconds (17.95 M allocations: 855.100 MiB, 14.50% gc time)
0.17525286988172276
If we do 50 radii it will take about 30 seconds
radlist=collect(0.001:0.001:0.049)
@time taulist=[collisiontime(radius=r) for r in radlist]
scatter(radlist,taulist,xlabel="radius",ylabel="τ",label="")
31.627615 seconds (895.04 M allocations: 41.504 GiB, 11.92% gc time, 0.04% compilation time)
Can we model this? Well each particle is moving with an average speed \(|v|=1\) in the units we are using. If we think of a ball randomly moving, in time \(dt\) will collide with anything other balls whose centers are within a strip that is \(v\,dt\) long and \(4 r\) wide. Thus the number of collisions you expect in that time is \(N_c= (v\,dt)(4 r)n\) where \(n\) is the density of balls. In this case \(n=100\). We find the collision time by setting \(N_c=1\),
tau_naive(r;density=100,vave=1)=1/(4*density*vave*r)
plot(tau_naive,0.001,0.05,label="")
scatter!(radlist,taulist,xlabel="radius",ylabel="τ",label="")
Plots.xlims!(0,0.05)
That is not bad – but it seems to miss things at both small and large radii. At small radii the discrepency is that we did not include collisions with the walls:
tau_better(r)=1/(1+1/tau_naive(r)) #r->1/(1+400*r)
plot(tau_better,0.001,0.05,label="")
scatter!(radlist,taulist,xlabel="radius",ylabel="τ",label="")
Another nice trick here is to plot the rate instead of the time
plot(x->1/tau_better(x),0.001,0.05,label="")
scatter!(radlist,1 ./taulist,xlabel="radius",ylabel="1/τ",label="")
Plots.ylims!(0,30)
Indeed, we seem to have captured what is happening at small radius.
The large radius discrepency has to do with the fact that in a time \(dt\) the ball will actually be able to hit other balls which are in a region of area \(A\sim(v dt+r)(4r)\). So the average time to get a collision with another ball is when
tau_naive_improved(r;density=100,vave=1)=1/vave*(1/(4*density*vave*r)-r)
#r->1/(1+400*r/(1-400*r^2))
tau_improved(r)=1/(1+1/tau_naive_improved(r))
plot(tau_improved,0.001,0.05,label="")
scatter!(radlist,taulist,xlabel="radius",ylabel="τ",label="")
Plots.ylims!(0,0.2)
Interestingly, it appears that there is some radius at which \(\tau\to0\). This is the ``Jamming’’ transition.
Pressure#
Another fun thing to look at is the pressure. The most physical way to calculate pressure is to look at the collisions of the balls with the walls. The pressure is the time averaged force exerted by these balls, normal to the wall, divided by the length of the wall. The force is exerted in little impulses at each collsion. So lets modify our simulation so that whenever there is a collision with the wall it stores the absolute value of the impulse. The pressure is proportional to the sum these impulses divided by the time interval.
function calculatepressure(billiard::Billiard,n)
setevents!(billiard)
startt=billiard.t
momentum=0.
for j in 2:n+1
setevents!(billiard,billiard.nextevent)
update!(billiard,billiard.nextevent)
next=billiard.nextevent
momentum+=wallmomentum(next.object1,next.object2)
end
endt=billiard.t
return momentum/(endt-startt)
end
wallmomentum(b1::Ball,b2::Ball)=0
wallmomentum(hw::Hwall,b::Ball)=abs(2*b.vy)
wallmomentum(vw::Vwall,b::Ball)=abs(2*b.vx)
wallmomentum(b::Ball,w::Wall)=wallmomentum(w,b)
wallmomentum (generic function with 4 methods)
p1=randomgrid(10,10,0.01,1,1)
@time calculatepressure(p1,1000) #equilibration run -- throw away
@time pressure1=calculatepressure(p1,5000)
0.118003 seconds (3.05 M allocations: 142.543 MiB, 8.86% compilation time)
0.651521 seconds (14.83 M allocations: 685.820 MiB, 26.82% gc time)
212.47473046608152
p2=randomgrid(10,10,0.048,1,1)
@time calculatepressure(p2,1000) #equilibration run -- throw away
@time pressure1=calculatepressure(p2,5000)
0.089052 seconds (3.18 M allocations: 147.409 MiB)
0.416132 seconds (15.72 M allocations: 720.496 MiB, 12.20% gc time)
4477.61131615306
Apparently the pressure quite strongly depends on the ball radius!
function pressure(rad)
b=randomgrid(10,10,rad,1,1)
calculatepressure(b,1000) # equilibrate
calculatepressure(b,5000) # calculate pressure
end
pressure (generic function with 1 method)
radlst=collect(0.01:0.001:0.049)
@time plst=[pressure(rad) for rad in radlst]
scatter(radlst,plst,xlabel="radius",ylabel="pressure",label="")
19.856300 seconds (745.31 M allocations: 33.524 GiB, 13.24% gc time, 0.11% compilation time)
scatter(radlst,1 ./plst,xlabel="radius",ylabel="1/P",label="")
Unfortunately, near the jamming transition the configurations can get stuck, and the result depends on the initial conditions. For example, the rad=0.049 data point is not indicative of a typical configuration
p3=randomgrid(10,10,0.049,1,1)
s3=simulate(p3,5000)
visualize(p3)
when the radius is 0.048, however, it can rearange itself
p4=randomgrid(10,10,0.048,1,1)
s4=simulate(p4,50000)
visualize(p4)
eventanimate(s4)
Task (runnable) @0x00000003eda0e270
There are tricks to avoid getting stuck, but fundamentally the system becomes “glassy” near the jamming transition – and one needs to do very careful calculations.
It turns out that somewhere before things jam up the hard disk gas even has a “hexatic” phase transition, which is a precursor to forming a liquid. We will not explore it.
Spatial Correlations#
If the Ergodic Hypothesis is correct, then every allowed configuration of hard spheres is equally likely. That then uniquely specified the velocity distribution. What does it tell us about the spatial distribution?
One challenge is how one even quantifies the spatial distribution. One convenient question to ask is about the distribution of the separations between the particles: \(g({\vec r})\) is the probability that 2 particles are separated by a vector \(\vec{r}\). This is referred to as the pair distribution function. If you have rotational symmetry, one expects that \(g\) will only depend on the magnitude of \(r\), and it makes sense to average it over angles
We can measure \(g(\vec{r})\) just like we did the speed distribution function – we just record all of the separation between particles, and make a histogram. As with the speeds, we can combine data from multiple times.
One slightly confusing feature is that we can consider a different function \(p(r)\) which is the probability that 2 particles are separated by distance \(r\).
If we bin the occurences of separation \(r\), we will have an estimate of \(p(r)\). Thus to get \(g(r)\) we need to divide by \(r\).
Rather than reporting back a big list of separations, it is better to bin the data as we go along. Lets make a Bin object, which can be called with data, storing it in a array
mutable struct Bin
min ::Float64
max ::Float64
step ::Float64
counts ::Array{Int64,1}
under ::Int64
over ::Int64
end
## Constructors
### Basic Constructor
function Bin(min,max,numbins)
step=(max-min)/numbins
counts=zeros(Int64,numbins)
Bin(min,max,step,counts,0,0)
end
### Named Constructor -- allowing one to specify either numbins or step
function Bin(;min,max,numbins=0,step=0)
if numbins==0 && step==0
throw("Bin must be constructed with either `numbin` or `step` arguments")
end
if step==0
return Bin(min,max,numbins)
end
numbins=Int64(round((max-min)/step))
MAX=min+numbins*step
return Bin(min,MAX,numbins)
end
function (b::Bin)(val)
if val<b.min
b.under+=1
return b
end
if val>b.max
b.over+=1
return b
end
bnum=1+Int64(floor((val-b.min)/b.step))
b.counts[bnum]+=1
return b
end
Base.push!(b::Bin,val)=b(val)
function hist(b::Bin,opts...;norm=false,namedopts...)
xvals=b.min:b.step:b.max
total = sum(b.counts)
scale=1
if norm
scale = total
end
bardata=[(xvals[j],b.counts[j]/scale) for j in 1:length(b.counts)]
bar(bardata,opts...;namedopts...)
end
hist (generic function with 3 methods)
tstbin=Bin(0.,10.,10)
Bin(0.0, 10.0, 1.0, [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 0, 0)
tstbin(6.6)
Bin(0.0, 10.0, 1.0, [0, 0, 0, 0, 0, 0, 3, 0, 0, 0], 0, 0)
tstbin(3.2)
Bin(0.0, 10.0, 1.0, [0, 0, 0, 1, 0, 0, 3, 0, 0, 0], 0, 0)
push!(tstbin,1.5)
Bin(0.0, 10.0, 1.0, [0, 2, 0, 1, 0, 0, 3, 0, 0, 0], 0, 0)
hist(tstbin)
Lets also make a N dimensional version.
mutable struct NBin{N,T}
min ::NTuple{N,T}
max ::NTuple{N,T}
step ::NTuple{N,T}
counts ::Array{Int64,N}
under ::Int64
over ::Int64
end
## Constructors
function NBin(min,max,numbins)
step= @. (max-min)/numbins
counts=zeros(Int64,Tuple(numbins))
NBin(Tuple(min),Tuple(max),Tuple(step),counts,0,0)
end
function NBin(;min,max,numbins=[],step=[])
if numbins==[] && step==[]
throw("NBin must be constructed with either `numbins` or `step` arguments")
end
if step==[]
return NBin(min,max,numbins)
end
numbins=@. Int64(round((max-min)/step))
MAX=@. min+numbins*step
NBin(Tuple(min),Tuple(MAX),Tuple(numbins))
end
function (b::NBin)(val)
if any(Tuple(val).<b.min)
b.under+=1
return b
end
if any(Tuple(val).>b.max)
b.over+=1
return b
end
bnum=@. 1+Int64(floor((val-b.min)/b.step))
b.counts[bnum...]+=1
return b
end
Base.push!(b::NBin,val)=b(val)
function hist(b::NBin{1,T},opts...;norm=false,rscale=false,namedopts...) where {T}
xvals=b.min[1]:b.step[1]:b.max[1]
total = sum(b.counts)
scale=1
if norm
scale = total
end
bardata=[(xvals[j],b.counts[j]/scale)
for j in 1:length(b.counts)]
if rscale
bardata=[(xvals[j],b.counts[j]/(scale*xvals[j]))
for j in 2:length(b.counts)]
end
bar(bardata,opts...;namedopts...)
end
function hist(b::NBin{2,T},opts...;norm=false,namedopts...) where {T}
total = sum(b.counts)
scale=1
if norm
scale = total
end
xvalues=collect(b.min[1]:b.step[1]:b.max[1])
yvalues=collect(b.min[2]:b.step[2]:b.max[2])
heatmap(xvalues,yvalues,b.counts./scale,opts...;namedopts...)
end
hist (generic function with 3 methods)
tstbin2=NBin([0.,0.],[1.,1.],[10,10])
data=[(rand(),rand()) for j in 1:100]
for d in data
push!(tstbin2,d)
end
tstbin2
NBin{2, Float64}((0.0, 0.0), (1.0, 1.0), (0.1, 0.1), [0 3 … 2 1; 0 0 … 3 0; … ; 0 0 … 2 0; 1 0 … 1 1], 0, 0)
hist(tstbin2,norm=true,aspect_ratio=:equal)
Now lets step through the simulation data, sampling the separations at each timestep:
function correlationtimeseries(s::Simulationdata;
dt,dx,dy,dr,xmin=-1/2,xmax=1/2,ymin=-1/2,ymax=1/2,
rmin=0.,rmax=1.)
# set up data structures
states=s.statelist
t=states[1].t
v_sep=NBin(min=(xmin,ymin),max=(xmax,ymax),step=(dx,dy))
s_sep=NBin(min=rmin,max=rmax,step=dr)
# step through collisions
for j in 2:length(states)
tj=states[j].t
# if collision is more than dt in the future
# then we want to sample the velocities,
# and increment t.
#
# We use a while loop, as we the next
# collision might be multiple dt's in the future
#
# We want to "eat up" all those timesteps
# before incrementing to the next collision
while tj>t+dt
t+=dt
prev=states[j-1]
nxt=states[j]
rbin!(v_sep,s_sep,prev,nxt,t)
end
end
return (v_sep=v_sep,s_sep=s_sep)
end
function rbin!(v_sep,s_sep,prev,nxt,t)
# Find positions at intermediate time
prevt=prev.t
nxtt=nxt.t
prevweight=(nxtt-t)/(nxtt-prevt)
nxtweight=(t-prevt)/(nxtt-prevt)
xvals=@. prev.xvalues*prevweight+nxt.xvalues*nxtweight
yvals=@. prev.yvalues*prevweight+nxt.yvalues*nxtweight
len=length(xvals)
dif=[0.,0.]
for j in 1:(len-1)
for i in j+1:len
dif .= (xvals[i]-xvals[j],yvals[i]-yvals[j])
v_sep(dif)
v_sep(-dif)
s_sep(norm(dif))
end
end
return nothing
end
rbin! (generic function with 1 method)
p1=randomgrid(10,10,0.01,1,1)
s1=simulate(p1,10000) #equilibrate
@time s1=simulate(p1,50000)
@time cs1=correlationtimeseries(s1,dt=0.5,dx=0.01,dy=0.01,dr=0.01,
xmin=-1.,ymin=-1.,xmax=1.,ymax=1.);
4.130292 seconds (149.00 M allocations: 6.898 GiB, 13.93% gc time)
3.771885 seconds (46.93 M allocations: 1.692 GiB, 4.49% gc time, 0.24% compilation time)
Since the x and y coordinates of each particle are bounded, the histogram of separations falls off with distance. This is solely a finite system size artifact. Also, rotational symmetry is broken by the walls.
hist(cs1.v_sep,norm=false,aspect_ratio=:equal,xlabel="x",ylabel="y")
Here is now a raw histogram of the separations between particles.
hist(cs1.s_sep,norm=false,rscale=false,
label="",xlabel="r",ylabel="Number of particles separated by r")
hist(cs1.s_sep,norm=true,rscale=false,label="",xlabel="r",ylabel="p(r)")
As already discussed, it makes sense to divide the vertical axis by \(r\), so we get the average number of particles at distance \(r\) in any one direction, \(g(r)\),
hist(cs1.s_sep,norm=false,rscale=true,label="",
xlabel="r",ylabel="g(r)")
In the thermodynamic limit that graph should be flat – but again, due to finite system size, there are fewer particles at large distances.
Now lets increase the size of the particles
p2=randomgrid(10,10,0.03,1,1)
s2=simulate(p2,10000) #equilibrate
@time s2=simulate(p2,50000)
@time cs2=correlationtimeseries(s2,dt=0.1,dx=0.01,dy=0.01,dr=0.01);
5.818289 seconds (155.39 M allocations: 7.147 GiB, 36.73% gc time)
5.155038 seconds (64.38 M allocations: 2.317 GiB, 3.84% gc time)
eventanimate(s2)
Task (runnable) @0x00000003770a84c0
hist(cs2.s_sep,norm=true,rscale=true,label="",xlabel="r",ylabel="g(r)")
The new thing that we are seeing is that there is a peak in \(g(r)\) at the particle radius, and a dip at roughly twice the particle radius. That means particles are more likely to have certain spacings. Lets see what that looks like in 2D
hist(cs2.v_sep,norm=false,aspect_ratio=:equal,xlim=[-0.2,0.2],ylim=[-0.2,0.2])
If we want more resolution we can average more frames, by making dt smaller, and take a smaller spacial grid
@time cs2hr=correlationtimeseries(s2,dt=0.01,dx=0.005,dy=0.005,dr=0.001);
51.633334 seconds (646.46 M allocations: 23.226 GiB, 3.83% gc time)
hist(cs2hr.v_sep,norm=false,aspect_ratio=:equal,xlim=[-0.2,0.2],ylim=[-0.2,0.2])
Things become even more interesting when we make the particle size larger yet
p3=randomgrid(10,10,0.04,1,1)
s3=simulate(p3,10000) #equilibrate
@time s3=simulate(p3,50000)
@time cs3=correlationtimeseries(s3,dt=0.1,dx=0.01,dy=0.01,dr=0.01);
4.395789 seconds (156.54 M allocations: 7.189 GiB, 14.58% gc time)
2.262194 seconds (27.84 M allocations: 1.002 GiB, 5.61% gc time)
hist(cs3.s_sep,norm=true,rscale=true,label="",xlabel="r",ylabel="g(r)")
We now have multiple bumps. This is characteristic of a liquid. The particles are roughly evenly spaced, so you are more likely to see particles at multiples of that spacing:
visualize(p3)
eventanimate(s3)
Task (runnable) @0x000000035378cb00
hist(cs3.v_sep,norm=false,aspect_ratio=:equal)
Another very clear feature is the horizontal and vertical bands in that correlation funciton plot. These are due to the boundaries. Just as particles are more likely to be found touching eachother, they are more likely to be found touching the walls.
Lets now make the radius even bigger
p4=randomgrid(10,10,0.048,1,1)
s4=simulate(p4,10000) #equilibrate
@time s4=simulate(p4,50000)
@time cs4=correlationtimeseries(s4,dt=0.01,dx=0.005,dy=0.005,dr=0.001);
4.294368 seconds (157.47 M allocations: 7.222 GiB, 13.72% gc time)
6.312658 seconds (79.63 M allocations: 2.860 GiB, 4.39% gc time)
hist(cs4.s_sep,norm=true,rscale=true)
eventanimate(s4)
Task (runnable) @0x0000000354bf0330
We now have a more complicated sequence of bumps. These are related to some sort of spatial ordering:
hist(cs4.v_sep,norm=true,aspect_ratio=:equal)
Here we have very distinct breaking of rotational symmetry. This is the sort of thing you tend to see in solids. If you look at the snapshots of the motion, you see definite “crystaline” regions.
visualize(p4)
Pair distribution function as a correlation function#
It is useful to connect the pair distribution function to a slightly more abstract quantity – the two particle correlation function. We make this connection to introduce the idea of correlation functions, and how they are probes of the physical system.
We first note that the probability of finding a particle at position \(\vec{r}\) can be written as an expectation value:
where \(\vec r_j\) is the position of the \(j\)’th particle. This expectation value can mean taking the average over many frames of our simulation. Note that \(P_1(\vec{r})\) can also be interpreted as the density of particles at \(\vec{r}\). We can write a single realization of the density as
and the ensemble average as
Similarly, our pair distribution function can be written as
It is telling us how correlated the density is at points separated by a displacement \(\vec{r}\).