Я использую недокументированную функцию contours в Matlab, чтобы получить нулевой контур функции, отображающей R^2 в R. Функция contours() вызывает contourc() для получения нулевого контура в терминах индексов массива, а contours() затем выполняет линейную интерполяцию для преобразования из (M,N) индексов массива в (Y,X) координаты данных. Это работает хорошо и дает нулевой контур, который по существу точен к точности станка. Однако попытка интерполировать вдоль этого нулевого контура для получения дополнительных точек не удалась, по-видимому, потому что интерполированные точки в незначительной степени отклоняются от истинного нулевого контура; Полученная ошибка слишком велика для предполагаемого приложения.
В качестве простого примера можно вычислить нулевой контур функции peaks, вычислить значения функции по контуру, а затем вычислить значения функции в новых точках, интерполированных по контуру. Максимальная ошибка составляет примерно eps(max(Z(:)) для точек, найденных на countours/countourc, и примерно на десять порядков выше для интерполированных точек.
% compute the zero contours; keep part of one contour
[X,Y,Z] = peaks(999);
XY0 = contours(X, Y, Z, [0 0]);
XY0 = XY0(:,140:(XY0(2,1)+1)); % keep only ascending values in the first contour for interpolation
% interpolate points along the zero contour
x2y = griddedInterpolant(XY0(1,:), XY0(2,:), 'linear','none');
X0i = linspace(XY0(1,1), XY0(1,end), 1e4);
Y0i = x2y(X0i);
% compute values of the function along the zero contour
Zi = interp2(X,Y,Z, XY0(1,:), XY0(2,:), 'linear', NaN);
Z0i = interp2(X,Y,Z, X0i, Y0i, 'linear', NaN);
% plot results
figure;
subplot(1,3,1); plot(XY0(1,:), XY0(2,:), '.'); hold on; plot(X0i, Y0i, 'Linewidth',1);
xlabel('X_0'); ylabel('Y_0'); title('(X_0,Y_0), (X_{0i},Y_{0i})');
subplot(1,3,2); plot(XY0(1,:), Zi, '.');
xlabel('X_0'); ylabel('Z_0'); title('f(X_0,Y_0)');
subplot(1,3,3); plot(XY0(1,:), Zi, '.'); hold on; plot(X0i, Z0i, '.');
xlabel('X_0'); ylabel('Z_0'); title('f(X_{0i},Y_{0i})');
Ошибка почти наверняка состоит в том, что интерполант действительно не следует нулевому контуру (то есть имеет одинаковую касательную, кривизну и, действительно, все более высокие производные в каждой точке). Таким образом, каждая интерполированная точка в небольшой степени отклоняется от пути нулевого контура. Это верно для всех проверенных методов интерполяции. Если бы была известна функциональная форма нулевого контура, можно было бы подогнать функцию к нулевому контуру, что, вероятно, дало бы удовлетворительные результаты. Однако для приложения, мотивирующего этот вопрос, аналитическое решение нулевого контура невозможно.
Учитывая неадекватность интерполяции для этой задачи, как можно получить новые точки вдоль нулевого контура? Точки должны быть такими, чтобы значение функции было так же близко к нулю, как и для набора точек, возвращаемого contourc. Глобальное увеличение разрешения сетки не вариант, так как затраты памяти и вычислений увеличиваются с квадратом разрешения. Локальная выборка с высоким разрешением итеративным образом на последовательных коротких интервалах нулевого контура является вариантом, но для ее реализации потребуется время, и я не знаю существующего кода (например, в FEX) для выполнения такой задачи. Кроме того, поскольку contourc является закрытым исходным кодом, адаптировать его к этой задаче не вариант.