Why are the particles getting off the Torus

1 次查看（过去 30 天）
David Miller
David Miller 2016-7-2
编辑： 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

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

回答（1 个）

José-Luis
José-Luis 2016-7-4
编辑：José-Luis 2016-7-4
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.

类别

在 Help Center 和 File Exchange 中查找有关 Particle Swarm 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by