Continuation of the preceding post, in which I considered the ODE system

with the initial conditions , where are distinct points in a Euclidean space. This can be seen as the gradient flow (in the space of all -point subsets) of the function . 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 be the time of the first collision. Then

where .

The lower bound follows at once from . To prove the upper bound on , we need the concept of monotonicity. A map is called *monotone* if for all . It is not hard to prove that the gradient of any convex function is monotone. In particular, is a monotone map, being the gradient of convex function (the lack of smoothness at the origin need not concern us here).

Renumbering the points, we may assume . Consider the function , . Differentiation yields

This inner product consists of the term

and the sum over of

Thus, for , which implies .

It is worth noting that if the process is stopped at the moment of first collision, it defines a map between symmetric products which fixes pointwise. The latter statement can now be made quantitative: if a set is at distance from , then . Indeed, in the notation of Lemma 1 we have . Since the first collision occurs at a time , it follows that .

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 nonlinear ODE is hardly a practical algorithm for clustering 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
```