Why are the particles getting off the Torus

Question

David Miller 2016-7-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/293445-why-are-the-particles-getting-off-the-torus

编辑： José-Luis 2016-7-4

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

José-Luis 2016-7-4

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/293445-why-are-the-particles-getting-off-the-torus#answer_227570

编辑：José-Luis 2016-7-4

在 MATLAB Online 中打开

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.