## Space-time (holiday) trees

No new math here, only new graphics. The following plot is from Collisions II, where it illustrates the process of contracting a 6-point planar set into a single point.

I remarked there that the trajectories cross in space, but not in space-time. The space-time plot of the contraction process is always a tree, but I never showed these trees to you.

Well, here they are. Time flows down, so that the trees have branches (initial points) above root.

The four points above were placed into the vertices of a rectangle. The plots below use points with random coordinates.

The Scilab command was

gravity4(grand(1,10,"def"),grand(1,10,"def"),250)

and the function itself is below.

function gravity4(vx,vy,Nsteps)
Npoints = length(vx);
h = 0.01/Npoints;
moving = ones(vx);
tolerance = h*(Npoints-1)/3;
x = zeros(Nsteps,Npoints);
y = zeros(Nsteps,Npoints);
z = zeros(Nsteps,Npoints);
x(1,:) = vx;
y(1,:) = vy;
for m = 1:(Nsteps-1)
x(m+1,:)= x(m,:);
y(m+1,:)= y(m,:);
z(m+1,:) = z(m,:);
for i = 1:Npoints
if moving(i)==0 then continue
end
z(m+1,i)=z(m,i)-h*Npoints;
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();
param3d1(x(:,:),y(:,:),z(:,:));
a = get("current_axes");
curves = a.children.children;
for i=1:Npoints
curves(i).foreground = 2+modulo(i,5);
curves(i).thickness = 4;
end
endfunction 

I tried to replace a.children.children with a.grandchildren, but that did not work…

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