For context, this is the original post that motivated this question.
Here is the code (please note that I am included the code as is, although the plots are generated by x-ing out of the graphics window when the simulation starts running, and letting the program open automatically new windows without the surface):
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.3010 0.7450 0.9330],'m',[0.8500 0.3250 0.0980]}; % 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.
eta = step_size / sqrt(dx0^2 + dy0^2); % Normalization.
eps = 0.0001; % EpsilonA
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.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', 4,'color',C{s})
drawnow
end
end
I am considering that, given how much these (geodesic) line simulations vary with minimal changes in the starting point, numerical approximations may account for the end result. It is confusing that there even seems to be a discrepancy in the number of points: the bunched up magenta and blue lines seem to run longer on Octave (right plot) than on Matlab (left plot).
It occurred to me that eps could be interfering with the default numerical precision, which was identical (both for Matlab and Octave). Getting rid of that line (and increasing the number of iterations):
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.3010 0.7450 0.9330],'m',[0.8500 0.3250 0.0980]}; % 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.
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,'color',C{s})
drawnow
end
end
But I don't think this is it, because at some points I got (with Octave):
And with Matlab, as I open the application and run the code above (without eps introduced as a variable), I get this:
Yet, if I run the very same code immediately afterwards also in Matlab, I get this:
Sometimes it takes multiple runs to get this last (prettier) image. And just to leave no uncertainty, the code used for the last two illustrations is:
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',...
'b','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.
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
