I am using the undocumented contours function in Matlab to obtain the zero contour of a function mapping R^2 to R. The function contours() calls contourc() to obtain the zero contour in terms of array indices, and contours() then performs linear interpolation to convert from (M,N) array indices to (Y,X) data coordinates. This works well and gives a zero contour which is essentially accurate to machine precision. However, attempting to interpolate along this zero contour to obtain additional points fails, apparently because the interpolated points deviate to a small extent from the true zero contour; the resulting error is too large for the intended application.
As a simple example, one can compute a zero contour of the peaks function, compute the function values along the contour, and then compute the function values at new points interpolated along the contour. The maximum error is roughly eps(max(Z(:)) for the points found by countours/countourc and about ten orders of magnitude higher for the interpolated points.
% 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})');
The error is almost certainly because the interpolant does not truly follow the zero contour (i.e., having the same tangent, curvature, and indeed all higher derivatives at each point). Each interpolated point thus deviates to a small degree from the path of the zero contour. This is true for all tested methods of interpolation. If the functional form of the zero contour were known, one could fit a function to the zero contour, which would likely give satisfactory results. However, for the application motivating this question, an analytical solution of the zero contour is not possible.
Given the inadequacy of interpolation to this task, how can new points along the zero contour be obtained? The points should be such that the value of the function is as close to zero as for the set of points returned by contourc. Increasing the grid resolution globally is not an option since the memory and computational expense increase with the square of the resolution. Locally sampling at high resolution in an iterative manner over consecutive short intervals of the zero contour is an option, but would take time to implement, and I'm not aware of existing code (e.g., in the FEX) to perform such a task. Also, since contourc is closed-source, adapting it to this task is not an option.
