## Collisions II

Continuation of the preceding post, in which I considered 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 conditions ${x_i(0)=p_i}$, where ${p_1,\dots,p_n}$ are distinct points in a Euclidean space. This can be seen as the gradient flow (in the space of all ${n}$-point subsets) of the function ${\sum_{i\ne j}|x_i-x_j|}$. Not surprisingly, minimization of the sum of distances leads to collisions. And they happen as quickly as they possibly could: here is a precise statement.

Lemma 1. Let ${t_c}$ be the time of the first collision. Then

${\displaystyle \frac{1}{2(n-1)}\,\delta\le t_c\le \frac{1}{2}\,\delta}$

where ${\displaystyle \delta=\min_{i\ne j}|p_i-p_j|}$.

The lower bound follows at once from ${|\dot{x}_i|\le n-1}$. To prove the upper bound on ${t_c}$, we need the concept of monotonicity. A map ${F\colon{\mathbb R}^d\rightarrow{\mathbb R}^d}$ is called monotone if ${\langle F(a)-F(b),a-b\rangle \ge 0}$ for all ${a,b\in{\mathbb R}^d}$. It is not hard to prove that the gradient of any convex function is monotone. In particular, ${F(x)=x/|x|}$ is a monotone map, being the gradient of convex function ${x\mapsto |x|}$ (the lack of smoothness at the origin need not concern us here).

Renumbering the points, we may assume ${|p_1-p_2|=\delta}$. Consider the function ${\varphi(t)=|x_1-x_2|}$, ${0. Differentiation yields

$\displaystyle \dot{\varphi}=|x_1-x_2|^{-1}\langle \dot{x}_1-\dot{x}_2, x_1-x_2\rangle.$

This inner product consists of the term

$\displaystyle \langle F(x_2-x_1)-F(x_1-x_2), x_1-x_2\rangle = -2|x_1-x_2|$

and the sum over ${j=3,\dots,n}$ of

$\displaystyle \langle F(x_j-x_1)-F(x_j-x_2), x_1-x_2\rangle \\ = - \langle F(x_j-x_1)-F(x_j-x_2),\ (x_j-x_1)-(x_j-x_2)\rangle \\ \le 0.$

Thus, ${\dot{\varphi} \le -2}$ for ${0, which implies ${t_c\le \varphi(0)/2}$. ${\Box}$

It is worth noting that if the process is stopped at the moment of first collision, it defines a map between symmetric products ${r\colon ({\mathbb R}^d)^{n}\rightarrow ({\mathbb R}^d)^{n-1}}$ which fixes ${({\mathbb R}^d)^{n-1}}$ pointwise. The latter statement can now be made quantitative: if a set ${A\in ({\mathbb R}^d)^{n}}$ is at distance ${\epsilon}$ from ${({\mathbb R}^d)^{n-1}}$, then ${d_H(r(A),A)\le (n-1)\,\epsilon}$. Indeed, in the notation of Lemma 1 we have ${\delta=2 \epsilon}$. Since the first collision occurs at a time ${t_c\le \epsilon}$, it follows that ${|x_i(t_c)-p_i|\le (n-1) t_c\le (n-1)\,\epsilon}$.

Time for some pictures. I changed the Scilab routine so that after each collision it drops one of the collided particles from the process, and restarts. As predicted, this did not change the qualitative picture. Main difference is that trajectories are now visibly non-smooth: angles are formed when a collision happens elsewhere.

In the above examples some trajectories cross in space, but not in space-time.

While the above plots may look like hierarchical clustering trees, solving a system of ${n}$ nonlinear ODE is hardly a practical algorithm for clustering ${n}$ points.

The updated scilab routine is below. This time I called it with random parameters:

 
 gravity2(grand(1,7,"unf",0,2),grand(1,7,"def"),300) 
function gravity2(vx,vy,Nsteps)
h = 0.002;
Npoints = length(vx);
moving = ones(vx);
tolerance = h*(Npoints-1)/3;
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
if moving(i)==0 then continue
end
for j = (i+1):Npoints
if moving(j)==0 then continue
end
v = [x(m,j),y(m,j)]-[x(m,i),y(m,i)];
if norm(v) < tolerance then
moving(j)=0;
else
x(m+1,i) = x(m+1,i) + h*v(1)/norm(v);
x(m+1,j) = x(m+1,j) - h*v(1)/norm(v);
y(m+1,i) = y(m+1,i) + h*v(2)/norm(v);
y(m+1,j) = y(m+1,j) - h*v(2)/norm(v);
end
end
end
end
clf();
plot(x(:,:),y(:,:),'o');
endfunction 

## 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