使用midpoint search method：
应用于函数 f(x,y) = x^3 + y^2，我将线段的点在 XY 平面 y = x 上从 x = -1 投影到 x = 1。
了解一下，一次迭代，XY平面上的直线上只有4个点，黑色球体是投影到曲面上的直线的这4个原始点，而红点是单次迭代的中点, 黄点是红点沿表面法线投影的结果：
使用 Matlab fmincon()，经过 5 次迭代，我们可以得到从 A 点到 B 点的测地线：
代码如下：
% Creating the surface
x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = x.^3 + y.^2;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap summer
% Number of points
n = 1000;
% Line to project on the surface with n values to get a feel for it...
t = linspace(-1,1,n);
height = t.^3 + t.^2;
P = [t;t;height];
% Plotting the projection of the line on the surface:
hold on
%plot3(P(1,:),P(2,:),P(3,:),'o')
for j=1:5
% First midpoint iteration updates P...
P = [P(:,1), (P(:,1:end-1) + P(:,2:end))/2, P(:,end)];
%plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize', 20)
A = zeros(3,size(P,2));
for i = 1:size(P,2)
% Starting point will be the vertical projection of the mid-points:
A(:,i) = [P(1,i), P(2,i), P(1,i)^3 + P(2,i)^2];
end
% Linear constraints:
nonlincon = @nlcon;
% Placing fmincon in a loop for all the points
for i = 1:(size(A,2))
% Objective function:
objective = @(x)(P(1,i) - x(1))^2 + (P(2,i) - x(2))^2 + (P(3,i)-x(3))^2;
A(:,i) = fmincon(objective, A(:,i), [], [], [], [], [], [], nonlincon);
end
P = A;
end
plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize', 5,'Color','y')
在一个名为nlcon.m的单独文件中：
function[c,ceq] = nlcon(x)
c = [];
ceq = x(3) - x(1)^3 - x(2)^2;
对于在 XY 上具有直线非对角线的非常凉爽的表面上的测地线也是如此：
% Creating the surface
x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = sin(3*(x.^2+y.^2))/10;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap summer
% Number of points
n = 1000;
% Line to project on the surface with n values to get a feel for it...
t = linspace(-1,1,n);
height = sin(3*((.5*ones(1,n)).^2+ t.^2))/10;
P = [(.5*ones(1,n));t;height];
% Plotting the line on the surface:
hold on
%plot3(P(1,:),P(2,:),P(3,:),'o')
for j=1:2
% First midpoint iteration updates P...
P = [P(:,1), (P(:,1:end-1) + P(:,2:end))/2, P(:,end)];
%plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize', 20)
A = zeros(3,size(P,2));
for i = 1:size(P,2)
% Starting point will be the vertical projection of the first mid-point:
A(:,i) = [P(1,i), P(2,i), sin(3*(P(1,i)^2+ P(2,i)^2))/10];
end
% Linear constraints:
nonlincon = @nonlincon;
% Placing fmincon in a loop for all the points
for i = 1:(size(A,2))
% Objective function:
objective = @(x)(P(1,i) - x(1))^2 + (P(2,i) - x(2))^2 + (P(3,i)-x(3))^2;
A(:,i) = fmincon(objective, A(:,i), [], [], [], [], [], [], nonlincon);
end
P = A;
end
plot3(P(1,:), P(2,:), P(3,:), '.', 'MarkerSize',5,'Color','r')
带有nonlincon.m中的非线性约束：
function[c,ceq] = nlcon(x)
c = [];
ceq = x(3) - sin(3*(x(1)^2+ x(2)^2))/10;
一个令人烦恼的问题是使用这种方法可能会过度拟合曲线，后一个图就是一个例子。所以我调整代码只选择一个起点和一个终点，并允许迭代过程找到曲线的其余部分，100 次迭代似乎朝着正确的方向前进：
上面的例子似乎是在XY平面上遵循线性投影，但幸运的是这不是一个固定的模式，这会进一步对方法产生怀疑。参见例如双曲抛物面 x^2 - y^2：
请注意，有一些算法可以沿表面 f(x,y) 推进或推动测地线，其增量由起点和表面的法线向量确定，如here。感谢 Alvise Vianello 在该模拟中研究 JS 的工作和他的sharing in GitHub，我能够将该算法转换为 Matlab 代码，为第一个示例生成此图，f(x,y) = x^3 + y^2：
这是 Matlab 代码：
x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = x.^3 + y.^2;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap('gray');
hold on
f = @(x,y) x.^3 + y.^2; % The actual surface
dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y
N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.
C = {'k','b','r','g','y','m','c',[.8 .2 .6],[.2,.8,.1],[0.3010 0.7450 0.9330],[0.9290 0.6940 0.1250],[0.8500 0.3250 0.0980]}; % Color scheme
for s = 1:11 % No. of lines to be plotted.
start = -5:5; % Distributing the starting points of the lines.
y0 = start(s)/5; % Fitting the starting pts between -1 and 1 along y axis.
x0 = 1; % Along x axis always starts at 1.
dx0 = 0; % Initial differential increment along x
dy0 = 0.05; % Initial differential increment along y
step_size = 0.000008; % Will determine the progression rate from pt to pt.
eta = step_size / sqrt(dx0^2 + dy0^2); % Normalization.
eps = 0.0001; % Epsilon
max_num_iter = 100000; % Number of dots in each line.
x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values
for i = 2:(max_num_iter - 1) % Creating the geodesic:
xt = x(i); % Values at point t of x, y and the function:
yt = y(i);
ft = f(xt,yt);
xtm1 = x(i - 1); % Values at t minus 1 (prior point) for x,y,f
ytm1 = y(i - 1);
ftm1 = f(xtm1,ytm1);
xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
ysymp = yt + (yt - ytm1);
fsymp = ft + (ft - ftm1);
df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
n = N(xt,yt); % Normal vector at point t
gamma = df * n(3); % Scalar x change f x z value of N
xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
ytp1 = ysymp - gamma * n(2);
x(i + 1) = xtp1;
y(i + 1) = ytp1;
end
P = [x; y; f(x,y)]; % Compiling results into a matrix.
indices = find(abs(P(1,:)) < 1); % Avoiding lines overshooting surface.
P = P(:,indices);
indices = find(abs(P(2,:)) < 1);
P = P(:,indices);
units = 15; % Deternines speed (smaller, faster)
packet = floor(size(P,2)/units);
P = P(:,1: packet * units);
for k = 1:packet:(packet * units)
hold on
plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
'.', 'MarkerSize', 3.5,'color',C{s})
drawnow
end
end
这是上面的一个较早示例，但现在计算方式不同，线并排开始，仅遵循测地线（无点对点轨迹）：
x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = sin(3*(x.^2+y.^2))/10;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap('gray');
hold on
f = @(x,y) sin(3*(x.^2+y.^2))/10; % The actual surface
dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y
N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.
C = {'k','r','g','y','m','c',[.8 .2 .6],[.2,.8,.1],[0.3010 0.7450 0.9330],[0.7890 0.5040 0.1250],[0.9290 0.6940 0.1250],[0.8500 0.3250 0.0980]}; % Color scheme
for s = 1:11 % No. of lines to be plotted.
start = -5:5; % Distributing the starting points of the lines.
x0 = -start(s)/5; % Fitting the starting pts between -1 and 1 along y axis.
y0 = -1; % Along x axis always starts at 1.
dx0 = 0; % Initial differential increment along x
dy0 = 0.05; % Initial differential increment along y
step_size = 0.00005; % Will determine the progression rate from pt to pt.
eta = step_size / sqrt(dx0^2 + dy0^2); % Normalization.
eps = 0.0001; % Epsilon
max_num_iter = 100000; % Number of dots in each line.
x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values
for i = 2:(max_num_iter - 1) % Creating the geodesic:
xt = x(i); % Values at point t of x, y and the function:
yt = y(i);
ft = f(xt,yt);
xtm1 = x(i - 1); % Values at t minus 1 (prior point) for x,y,f
ytm1 = y(i - 1);
ftm1 = f(xtm1,ytm1);
xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
ysymp = yt + (yt - ytm1);
fsymp = ft + (ft - ftm1);
df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
n = N(xt,yt); % Normal vector at point t
gamma = df * n(3); % Scalar x change f x z value of N
xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
ytp1 = ysymp - gamma * n(2);
x(i + 1) = xtp1;
y(i + 1) = ytp1;
end
P = [x; y; f(x,y)]; % Compiling results into a matrix.
indices = find(abs(P(1,:)) < 1); % Avoiding lines overshooting surface.
P = P(:,indices);
indices = find(abs(P(2,:)) < 1);
P = P(:,indices);
units = 35; % Deternines speed (smaller, faster)
packet = floor(size(P,2)/units);
P = P(:,1: packet * units);
for k = 1:packet:(packet * units)
hold on
plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))), '.', 'MarkerSize', 5,'color',C{s})
drawnow
end
end
更多示例：
x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = x.^2 - y.^2;
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap('gray');
f = @(x,y) x.^2 - y.^2; % The actual surface
dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y
N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.
C = {'b','w','r','g','y','m','c',[0.75, 0.75, 0],[0.9290, 0.6940, 0.1250],[0.3010 0.7450 0.9330],[0.1290 0.6940 0.1250],[0.8500 0.3250 0.0980]}; % Color scheme
for s = 1:11 % No. of lines to be plotted.
start = -5:5; % Distributing the starting points of the lines.
x0 = -start(s)/5; % Fitting the starting pts between -1 and 1 along y axis.
y0 = -1; % Along x axis always starts at 1.
dx0 = 0; % Initial differential increment along x
dy0 = 0.05; % Initial differential increment along y
step_size = 0.00005; % Will determine the progression rate from pt to pt.
eta = step_size / sqrt(dx0^2 + dy0^2); % Normalization.
eps = 0.0001; % Epsilon
max_num_iter = 100000; % Number of dots in each line.
x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values
for i = 2:(max_num_iter - 1) % Creating the geodesic:
xt = x(i); % Values at point t of x, y and the function:
yt = y(i);
ft = f(xt,yt);
xtm1 = x(i - 1); % Values at t minus 1 (prior point) for x,y,f
ytm1 = y(i - 1);
ftm1 = f(xtm1,ytm1);
xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
ysymp = yt + (yt - ytm1);
fsymp = ft + (ft - ftm1);
df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
n = N(xt,yt); % Normal vector at point t
gamma = df * n(3); % Scalar x change f x z value of N
xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
ytp1 = ysymp - gamma * n(2);
x(i + 1) = xtp1;
y(i + 1) = ytp1;
end
P = [x; y; f(x,y)]; % Compiling results into a matrix.
indices = find(abs(P(1,:)) < 1); % Avoiding lines overshooting surface.
P = P(:,indices);
indices = find(abs(P(2,:)) < 1);
P = P(:,indices);
units = 45; % Deternines speed (smaller, faster)
packet = floor(size(P,2)/units);
P = P(:,1: packet * units);
for k = 1:packet:(packet * units)
hold on
plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))), '.', 'MarkerSize', 5,'color',C{s})
drawnow
end
end
或者这个：
x = linspace(-1,1);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = .07 * (.1 + x.^2 + y.^2).^(-1);
S = [x;y;z];
h = surf(x,y,z)
zlim([0 8])
set(h,'edgecolor','none')
colormap('gray');
axis off
hold on
f = @(x,y) .07 * (.1 + x.^2 + y.^2).^(-1); % The actual surface
dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y
N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.
C = {'w',[0.8500, 0.3250, 0.0980],[0.9290, 0.6940, 0.1250],'g','y','m','c',[0.75, 0.75, 0],'r',...
[0.56,0,0.85],'m'}; % Color scheme
for s = 1:10 % No. of lines to be plotted.
start = -9:2:9;
x0 = -start(s)/10;
y0 = -1; % Along x axis always starts at 1.
dx0 = 0; % Initial differential increment along x
dy0 = 0.05; % Initial differential increment along y
step_size = 0.00005; % Will determine the progression rate from pt to pt.
eta = step_size / sqrt(dx0^2 + dy0^2); % Normalization.
eps = 0.0001; % EpsilonA
max_num_iter = 500000; % Number of dots in each line.
x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values
for i = 2:(max_num_iter - 1) % Creating the geodesic:
xt = x(i); % Values at point t of x, y and the function:
yt = y(i);
ft = f(xt,yt);
xtm1 = x(i - 1); % Values at t minus 1 (prior point) for x,y,f
ytm1 = y(i - 1);
ftm1 = f(xtm1,ytm1);
xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
ysymp = yt + (yt - ytm1);
fsymp = ft + (ft - ftm1);
df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
n = N(xt,yt); % Normal vector at point t
gamma = df * n(3); % Scalar x change f x z value of N
xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
ytp1 = ysymp - gamma * n(2);
x(i + 1) = xtp1;
y(i + 1) = ytp1;
end
P = [x; y; f(x,y)]; % Compiling results into a matrix.
indices = find(abs(P(1,:)) < 1.5); % Avoiding lines overshooting surface.
P = P(:,indices);
indices = find(abs(P(2,:)) < 1);
P = P(:,indices);
units = 15; % Deternines speed (smaller, faster)
packet = floor(size(P,2)/units);
P = P(:,1: packet * units);
for k = 1:packet:(packet * units)
hold on
plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
'.', 'MarkerSize', 3.5,'color',C{s})
drawnow
end
end
或者一个sinc函数：
x = linspace(-10, 10);
y = linspace(-10, 10);
[x,y] = meshgrid(x,y);
z = sin(1.3*sqrt (x.^ 2 + y.^ 2) + eps)./ (sqrt (x.^ 2 + y.^ 2) + eps);
S = [x;y;z];
h = surf(x,y,z)
set(h,'edgecolor','none')
colormap('gray');
axis off
hold on
f = @(x,y) sin(1.3*sqrt (x.^ 2 + y.^ 2) + eps)./ (sqrt (x.^ 2 + y.^ 2) + eps); % The actual surface
dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y
N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.
C = {'w',[0.8500, 0.3250, 0.0980],[0.9290, 0.6940, 0.1250],'g','y','r','c','m','w',...
[0.56,0,0.85],[0.8500, 0.7250, 0.0980],[0.2290, 0.1940, 0.6250],'w',...
[0.890, 0.1940, 0.4250],'y',[0.2290, 0.9940, 0.3250],'w',[0.1500, 0.7250, 0.0980],...
[0.8500, 0.3250, 0.0980],'m','w'}; % Color scheme
for s = 1:12 % No. of lines to be plotted.
x0 = 10;
y0 = 10; % Along x axis always starts at 1.
dx0 = -0.001*(cos(pi /2 *s/11)); % Initial differential increment along x
dy0 = -0.001*(sin(pi /2 *s/11)); % Initial differential increment along y
step_size = 0.0005; % Will determine the progression rate from pt to pt.
% Making it smaller increases the length of the curve.
eta = step_size / sqrt(dx0^2 + dy0^2); % Normalization.
eps = 0.0001; % EpsilonA
max_num_iter = 500000; % Number of dots in each line.
x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values
for i = 2:(max_num_iter - 1) % Creating the geodesic:
xt = x(i); % Values at point t of x, y and the function:
yt = y(i);
ft = f(xt,yt);
xtm1 = x(i - 1); % Values at t minus 1 (prior point) for x,y,f
ytm1 = y(i - 1);
ftm1 = f(xtm1,ytm1);
xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
ysymp = yt + (yt - ytm1);
fsymp = ft + (ft - ftm1);
df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
n = N(xt,yt); % Normal vector at point t
gamma = df * n(3); % Scalar x change f x z value of N
xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
ytp1 = ysymp - gamma * n(2);
x(i + 1) = xtp1;
y(i + 1) = ytp1;
end
P = [x; y; f(x,y)]; % Compiling results into a matrix.
indices = find(abs(P(1,:)) < 10); % Avoiding lines overshooting surface.
P = P(:,indices);
indices = find(abs(P(2,:)) < 10);
P = P(:,indices);
units = 15; % Deternines speed (smaller, faster)
packet = floor(size(P,2)/units);
P = P(:,1: packet * units);
for k = 1:packet:(packet * units)
hold on
plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
'.', 'MarkerSize', 3.5,'color',C{s})
drawnow
end
end
最后一个：
x = linspace(-1.5,1.5);
y = linspace(-1,1);
[x,y] = meshgrid(x,y);
z = 0.5 *y.*sin(5 * x) - 0.5 * x.*cos(5 * y)+1.5;
S = [x;y;z];
h = surf(x,y,z)
zlim([0 8])
set(h,'edgecolor','none')
colormap('gray');
axis off
hold on
f = @(x,y) 0.5 *y.* sin(5 * x) - 0.5 * x.*cos(5 * y)+1.5; % The actual surface
dfdx = @(x,y) (f(x + eps, y) - f(x - eps, y))/(2 * eps); % ~ partial f wrt x
dfdy = @(x,y) (f(x, y + eps) - f(x, y - eps))/(2 * eps); % ~ partial f wrt y
N = @(x,y) [- dfdx(x,y), - dfdy(x,y), 1]; % Normal vec to surface @ any pt.
C = {'w',[0.8500, 0.3250, 0.0980],[0.9290, 0.6940, 0.1250],'g','y','k','c',[0.75, 0.75, 0],'r',...
[0.56,0,0.85],'m'}; % Color scheme
for s = 1:11 % No. of lines to be plotted.
start = [0, 0.7835, -0.7835, 0.5877, -0.5877, 0.3918, -0.3918, 0.1959, -0.1959, 0.9794, -0.9794];
x0 = start(s);
y0 = -1; % Along x axis always starts at 1.
dx0 = 0; % Initial differential increment along x
dy0 = 0.05; % Initial differential increment along y
step_size = 0.00005; % Will determine the progression rate from pt to pt.
% Making it smaller increases the length of the curve.
eta = step_size / sqrt(dx0^2 + dy0^2); % Normalization.
eps = 0.0001; % EpsilonA
max_num_iter = 500000; % Number of dots in each line.
x = [[x0, x0 + eta * dx0], zeros(1,max_num_iter - 2)]; % Vec of x values
y = [[y0, y0 + eta * dy0], zeros(1,max_num_iter - 2)]; % Vec of y values
for i = 2:(max_num_iter - 1) % Creating the geodesic:
xt = x(i); % Values at point t of x, y and the function:
yt = y(i);
ft = f(xt,yt);
xtm1 = x(i - 1); % Values at t minus 1 (prior point) for x,y,f
ytm1 = y(i - 1);
ftm1 = f(xtm1,ytm1);
xsymp = xt + (xt - xtm1); % Adding the prior difference forward:
ysymp = yt + (yt - ytm1);
fsymp = ft + (ft - ftm1);
df = fsymp - f(xsymp,ysymp); % Is the surface changing? How much?
n = N(xt,yt); % Normal vector at point t
gamma = df * n(3); % Scalar x change f x z value of N
xtp1 = xsymp - gamma * n(1); % Gamma to modulate incre. x & y.
ytp1 = ysymp - gamma * n(2);
x(i + 1) = xtp1;
y(i + 1) = ytp1;
end
P = [x; y; f(x,y)]; % Compiling results into a matrix.
indices = find(abs(P(1,:)) < 1.5); % Avoiding lines overshooting surface.
P = P(:,indices);
indices = find(abs(P(2,:)) < 1);
P = P(:,indices);
units = 15; % Deternines speed (smaller, faster)
packet = floor(size(P,2)/units);
P = P(:,1: packet * units);
for k = 1:packet:(packet * units)
hold on
plot3(P(1, k:(k+packet-1)), P(2,(k:(k+packet-1))), P(3,(k:(k+packet-1))),...
'.', 'MarkerSize', 3.5,'color',C{s})
drawnow
end
end