# Clustering by collisions

This construction looks like it should answer the question raised at the end of the previous post, but instead of checking whether it really works, I decided to post some pictures.

Given ${n}$ distinct points ${p_1,\dots,p_n}$ in a Euclidean space ${{\mathbb R}^d}$, consider the ODE system

$\displaystyle \dot{x}_i = \sum_{j\ne i}\frac{x_j-x_i}{|x_j-x_i|},\quad i=1,\dots,n$

with the initial condition ${x_i(0)=p_i}$. The solution exists and is well-behaved until the first collision: ${x_i=x_j}$ for some ${i\ne j}$. We can stop the process then (thus obtaining a set of at most ${n-1}$ points) or remove the duplicates and restart. Let’s try the latter approach. To simplify things (at a cost), we can write the system as follows:

$\displaystyle \dot{x}_i = \sum_{x_j\ne x_j}\frac{x_j-x_i}{|x_j-x_i|},\quad i=1,\dots,n.$

Now the solution exists for all times, but becomes stationary when all points come together. The cost to which I alluded above is that keeping multiple copies of collided points distorts the post-collision behavior, but I hope that the qualitative picture survives.

I used Scilab to plot the trajectories of this system for several initial configurations. For purposes of numeric computation, the ODE was modified once again:

$\displaystyle \dot{x}_i = \sum_{|x_j- x_j|>\epsilon }\frac{x_j-x_i}{|x_j-x_i|},\quad i=1,\dots,n.$

where ${\epsilon}$ is tolerance, related to the step size used by the solver (Euler’s method works fine).

In each case the particles obviously want to collide: the points placed nearby move toward each other. Indeed, suppose that ${p_1}$ and ${p_2}$ are much close to each other than to the other points. The ODE system yields

$\displaystyle \frac{d}{dt}(x_1-x_2) = -2\,\frac{x_1-x_2}{|x_1-x_2|} + \text{(small terms)}$

where the smallness is due to the near cancellation of ${\displaystyle \frac{x_j-x_1}{|x_j-x_1|}}$ and ${\displaystyle \frac{x_j-x_2}{|x_j-x_2|}}$ for ${j\ne 1,2}$. As long as the sum of small terms is bounded by a constant less than ${2}$, the distance ${|x_1-x_2|}$ will decrease at linear rate, heading toward zero.

Here is the (hastily written) scilab function. It takes three parameters: vector of x-coordinates, vector of y-coordinates, and the number of steps. I use 100 for the latter, increasing it if the points do not have enough time to collide. For example, try gravity([0,1,2,3],[2,3,0,0.5],100).


function gravity(vx,vy,Nsteps)
h = 0.01;
Npoints = length(vx);
tolerance = h*(Npoints-1)/2;
x = zeros(Nsteps,Npoints);
y = zeros(Nsteps,Npoints);
x(1,:) = vx;
y(1,:) = vy;
for m = 1:(Nsteps-1)
x(m+1,:)= x(m,:);
y(m+1,:)= y(m,:);
for i = 1:Npoints
for j = 1:Npoints
v = [x(m,j),y(m,j)]-[x(m,i),y(m,i)];
if norm(v) > tolerance then
x(m+1,i) = x(m+1,i) + h*v(1)/norm(v);
y(m+1,i) = y(m+1,i) + h*v(2)/norm(v);
end
end
end
end
clf();
plot(x(:,:),y(:,:));
endfunction