Why are the particles getting off the Torus

I have a program that randomly places N particles on a torus and each particle chases the next one (particle 1 chases particle 2, particle 2 chases particle 3..., particle N chases particle 1). I have three functions: eqn function: returns (x,y,z) position when (theta, phi) are passed in, projTangent: projects the direction particle m has (vector joining particle m and particle m+1) onto its tangent space, projTorus: projects the particle back onto the torus to its closest point. prams is a struct I created that has all the constants the user can change (such as the constants in the parameterization of the torus, number of particles, etc). After the first time step, some of the particles are being plotted off the torus and keep getting further away as there are more time steps. Can anyone shed some light on this situation?

%DRIVER 
prams.NBegin = 20;
prams.N = 25;
prams.Radius = 2;
prams.radius = 1;
prams.perturbations = 10;
prams.runs = 1;
prams.dt = 0.1;
prams.ntime = 2500;
prams.newtonSteps = 3;
options.iplot = true;

 if prams.Radius > prams.radius
   torusResults = zeros(prams.N-prams.NBegin+1, 1); 

for N = prams.NBegin:prams.N
  prams.N = N;
  disp(N);
  for m = 1:prams.runs
    [posFinal] = bugs_torus(options, prams);
  end
    torusResults(prams.N-prams.NBegin+1,1) = total/prams.runs
end
elseif prams.Radius == prams.radius
  fprintf('This is a Horn Torus! Program Terminated!\n');
else
  fprintf('This is a Spindle Torus! Program  Terminated!\n');  
end

The next segment of code will be the function called in the driver.

function [posFinal] = bugs_torus(options, prams)

function [x,y,z] = eqn(theta, phi)
  x = (prams.Radius + prams.radius*cos(theta)).*cos(phi);
  y = (prams.Radius + prams.radius*cos(theta)).*sin(phi);
  z = prams.radius*sin(theta);
end

function vel = projTangent(vector, theta, phi)
  x = -1*sin(phi);
  y = cos(phi);
  z = 0;
  X = -1*sin(theta).*cos(phi);
  Y = -1*sin(theta).*sin(phi);
  Z  = cos(theta);  
  t1 = [x; y; z];
  t2 = [X; Y; Z];
  normal = cross(t1, t2);
  vel = vector - normal*dot(vector, normal);
end

function [newTheta, newPhi] = projTorus(theta, phi, x, y, z)
  pair = [theta; phi];
  for i = 1:prams.newtonSteps
    theta = pair(1,1);
    phi = pair(2,1);
    dF_dTheta = 2*prams.radius*(x*sin(theta).*cos(phi) + y*sin(theta).*sin(phi) - ...
      (prams.Radius  + prams.radius*cos(theta)).*sin(theta) + ...
      cos(theta).*(prams.radius*sin(theta) - z));
    dF_dPhi = 2*(prams.Radius + prams.radius*cos(theta)).*(x*sin(phi) - y*cos(phi));
    d2F_dThetadPhi = -2*prams.radius*sin(theta).*(x*sin(phi) - y*cos(phi));
    d2f_d2Theta = 2*prams.radius*(x*cos(theta).*cos(phi) - y*sin(theta).*sin(phi) - ...
      prams.Radius*cos(theta) + z*sin(theta));
    d2F_d2Phi = (prams.Radius + prams.radius*cos(theta)).*(x*cos(phi) + y*sin(phi));
    F = [dF_dTheta; dF_dPhi];
    inverseJacobian = inv([d2f_d2Theta d2F_dThetadPhi; d2F_dThetadPhi d2F_d2Phi]);
    pair = pair - inverseJacobian*F;
  end
  newTheta = pair(1,1);
  newPhi = pair(2,1);
end

%Initialize the positions of the particles
%angles is a 2xN matrix  where the first row is the theta angle
%and the second row is the phi angle. Column m corresponds to particle m.
angles = zeros(2, prams.N);
for k = 1:prams.N 
  angles(1, k) = rand*2*pi;
  angles(2, k) = rand*2*pi;
  thetaInit = angles(1, k);
  phiInit = angles(2, k);
  [xInit,yInit,zInit] = eqn(thetaInit, phiInit);
  posInit(1, k) = xInit
  posInit(2, k) = yInit
  posInit(3, k) = zInit
end

pos = posInit;

vel = zeros(3, prams.N);

for j = 1:prams.ntime

    for k = 1:prams.N-1
      direction = pos(:, k+1) - pos(:, k);
      vel(:, k) = projTangent(direction, angles(1,k), angles(2,k));
    end
    direction = pos(:, 1) - pos(:, prams.N);
    vel(:, prams.N) = projTangent(direction, angles(1,prams.N), angles(2,prams.N));

    %Plotting the simulation

    if options.iplot
      hold off;
      [thetaGrid,phiGrid] = meshgrid(0:pi/50:2*pi);
      xAxis = (prams.Radius + prams.radius*cos(thetaGrid)).*cos(phiGrid);
      yAxis = (prams.Radius + prams.radius*cos(thetaGrid)).*sin(phiGrid);
      zAxis= prams.radius*sin(thetaGrid);
      surf(xAxis,yAxis,zAxis);
      hold on;
      plot3(pos(1,:), pos(2,:), pos(3,:), 'r.', 'markersize', 20);
      hold on;
      axis equal
      axis([-(prams.Radius + prams.radius + .1) (prams.Radius + prams.radius + .1) ...
        -(prams.Radius + prams.radius + .1) (prams.Radius + prams.radius + .1) ...
        -(prams.radius + .1) (prams.radius + .1)]);
      pause();
    end

    pos = pos + vel*prams.dt;

    for i = 1:prams.N
      thetaOld = angles(1,i);
      phiOld = angles(2,i);
      xVel = pos(1,i);
      yVel = pos(2,i);
      zVel = pos(3,i);
      [newTheta, newPhi] = projTorus(thetaOld, phiOld, xVel, yVel, zVel);
      angles(1, i) = newTheta;
      angles(2, i) = newPhi;
      [xNew,yNew,zNew] = eqn(newTheta, newPhi);
      pos(1,i) = xNew;
      pos(2,i) = yNew;
      pos(3,i) = zNew;
    end

end

posFinal = pos;

end

Answer 1

[José-Luis]

Without going through all that code, one thing pops to mind: numerical precision and stability. Computational geometry operations are particularly prone to that.

Just try

 cos(pi/2)

Should be zero right? It isn't. Such approximations may introduce perturbations that can become larger and larger, leading to your points merrily going away.