Блог

Я использую недокументированную 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 имеет закрытый исходный код, адаптировать его к этой задаче невозможно.

Я решил эту проблему, взяв интерполированные (X,Y) значения в качестве приближения к истинному нулевому контуру и повторно интерпретировав Y координату каждой пары относительно Z значений поверхности, примыкающей к приблизительному нулю, чтобы определить, при какой Y для каждой Z=0 пары (X,Y координате. При этом сохраняется равномерный интервал между интерполированными X значениями.

Для каждой (iX,iY) пары, где iX и iY обозначают (дробные) индексы массива в X и Y координатных сетках для точек, интерполированных вдоль нулевого контура, приблизительные iY значения повторно интерпретируются относительно значений Z в позициях, расположенных по бокам (в Y измерении) приблизительного нуля, т.е. в (iX,floor(iY)) и (iX,ceil(iY)) . Z Значения в этих граничных точках сами по себе получены путем выполнения билинейной интерполяции Z в этих точках.

После повторной интерпретации величина ошибки такая же, как у точек, первоначально возвращенных с помощью contourc / contours .