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<t<t_c}. 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<t<t_c}, 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.

4 points
4 points
5 points
5 points
6 points
6 points

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

7 points
7 points

8 points
8 points
9 points
9 points

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 

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s