Для контекста, это оригинальный пост , который мотивировал этот вопрос.
Вот код (пожалуйста, обратите внимание что я включил код как есть, хотя графики генерируются путем выклинивания графического окна, когда симуляция запускается, и позволяя программе автоматически открывать новый windows без поверхности):
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
Я считаю, что с учетом того, насколько эти (геодезические c) линии моделирования изменяются при минимальных изменениях в начальной точке, численные приближения могут составлять конечный результат. Это сбивает с толку, что даже кажется, что есть расхождение в количестве точек: сгруппированные пурпурные и синие линии кажутся более длинными в Октаве (правый график), чем в Матлабе (левый график).
Мне пришло в голову, что eps может мешать точности по умолчанию, которая была одинаковой (как для Matlab, так и для Octave). Избавление от этой строки (и увеличение количества итераций):
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
Но я не думаю, что это так, потому что при некоторые очки, которые я получил (с Octave):
И с Matlab, когда я открываю приложение и запускаю приведенный выше код (без введенного eps в качестве переменной), я получаю это:
Тем не менее, если я запускаю тот же самый код сразу после этого также в Matlab, я получаю это:
Иногда требуется несколько прогонов, чтобы получить это последнее (красивее) изображение. И чтобы не было неопределенности, код, использованный для двух последних иллюстраций:
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