In [1]:
import numpy
import tqdm
import wendy
%pylab inline
import matplotlib.animation as animation
from IPython.display import HTML
_SAVE_GIFS= False
rcParams.update({'axes.labelsize': 17.,
              'font.size': 12.,
              'legend.fontsize': 17.,
              'xtick.labelsize':15.,
              'ytick.labelsize':15.,
              'text.usetex': _SAVE_GIFS,
              'figure.figsize': [5,5],
              'xtick.major.size' : 4,
              'ytick.major.size' : 4,
              'xtick.minor.size' : 2,
              'ytick.minor.size' : 2,
              'legend.numpoints':1})
numpy.random.seed(2)
Populating the interactive namespace from numpy and matplotlib

Phase-mixing and violent relaxation in one dimension

We study gravitational collapse in one dimension, illustrating phase mixing and violent relaxation using the setup of Binney (2004) and Schulz et al. (2013). We start with a grid of particles on a uniform grid in $x$ with a small velocity perturbation

$x \sim [-\pi/2,\pi/2]\,,$

$v = -V_0\,\sin(x)+\epsilon\,V_0\,\mathcal{N}(0,1)\,,$

where we also add a very small $\epsilon \approx10^{-7}$ additional perturbation to the velocities to avoid numerical instabilities. We integrate this system with $N = 1001$ particles for $\approx132$ time steps:

In [2]:
N= 1001
dx= numpy.pi
V0= 0.001

x= dx*(numpy.arange(N)-N//2)/N
v= -V0*numpy.sin(x)+numpy.random.normal(size=N)*V0*10.**-7.
m= numpy.ones(N)/float(N)
In [3]:
g= wendy.nbody(x,v,m,0.05,maxcoll=100000000)
In [4]:
nt= 2700
xt= numpy.empty((N,nt+1))
vt= numpy.empty((N,nt+1))
Et= numpy.empty((nt+1))
xt[:,0]= x
vt[:,0]= v
Et[0]= wendy.energy(x,v,m)
for ii in tqdm.trange(nt):
    tx,tv= next(g)
    xt[:,ii+1]= tx
    vt[:,ii+1]= tv
    Et[ii+1]= wendy.energy(tx,tv,m)

First, we check that the total energy is conserved during the $N$-body integration:

In [5]:
plot(numpy.arange(len(Et))/20.,Et)
xlabel(r'$t$')
ylabel(r'$E_{\mathrm{tot}}$')
Out[5]:
<matplotlib.text.Text at 0x112e82210>

Then we make the equivalent of Figure 1 in Schulz et al. (2013), showing the sytem's phase space at 5 different times:

In [6]:
figsize(6,6)
subplot(2,2,1)
scatter(xt[:,0],vt[:,0],c=numpy.fabs(x),s=1.,edgecolors='None')
scatter(xt[:,360],vt[:,360],c=numpy.fabs(x),s=1.,edgecolors='None')
ylabel(r'$v$')
xlim(-3.49,3.49)
ylim(-2.49,2.49)
annotate(r'$t=0\,\mathrm{and}\ t=18$',(0.95,0.95),xycoords='axes fraction',
         horizontalalignment='right',verticalalignment='top',size=18.)
subplot(2,2,2)
scatter(xt[:,500],vt[:,500],c=numpy.fabs(x),s=1.,edgecolors='None')
xlim(-3.49,3.49)
ylim(-2.49,2.49)
annotate(r'$t=25$',(0.95,0.95),xycoords='axes fraction',
         horizontalalignment='right',verticalalignment='top',size=18.)
subplot(2,2,3)
scatter(xt[:,800],vt[:,800],c=numpy.fabs(x),s=1.,edgecolors='None')
xlabel(r'$x$')
ylabel(r'$v$')
xlim(-3.49,3.49)
ylim(-2.49,2.49)
annotate(r'$t=40$',(0.95,0.95),xycoords='axes fraction',
         horizontalalignment='right',verticalalignment='top',size=18.)
subplot(2,2,4)
scatter(xt[:,2640],vt[:,2640],c=numpy.fabs(x),s=1.,edgecolors='None')
xlabel(r'$x$')
xlim(-3.49,3.49)
ylim(-2.49,2.49)
annotate(r'$t=132$',(0.95,0.95),xycoords='axes fraction',
         horizontalalignment='right',verticalalignment='top',size=18.)
tight_layout()

We can also make a movie of the evolution of the system in phase space:

In [7]:
def init_anim_frame():
    line1= plot([],[])
    xlabel(r'$x$')
    ylabel(r'$v$')
    xlim(-3.49,3.49)
    ylim(-2.49,2.49)
    return (line1[0],)
figsize(6,4)
fig, ax= subplots()
line= ax.scatter(x,v,c=numpy.fabs(x),s=5.,edgecolors='None')
txt= ax.annotate(r'$t=%.0f$' % (0.),
                 (0.95,0.95),xycoords='axes fraction',
                 horizontalalignment='right',verticalalignment='top',size=18.)
subsamp= 4
def animate(ii):
    line.set_offsets(numpy.array([xt[:,ii*subsamp],vt[:,ii*subsamp]]).T)
    txt.set_text(r'$t=%.0f$' % (ii*subsamp/20.))
    return (line,)
anim = animation.FuncAnimation(fig,animate,init_func=init_anim_frame,
                               frames=nt//subsamp,interval=40,blit=True,repeat=True)
if _SAVE_GIFS:
    anim.save('schulz_sinx.gif',writer='imagemagick',dpi=80)
# The following is necessary to just get the movie, and not an additional initial frame
plt.close()
out= HTML(anim.to_html5_video())
plt.close()
out
Out[7]: