【发布时间】:2019-03-26 23:18:20
【问题描述】:
我正在使用 Matlab 中未记录的 contours 函数来获取函数映射 R^2 到 R 的零轮廓。函数contours() 调用contourc() 以获取数组索引方面的零轮廓,然后contours() 执行线性插值以将(M,N) 数组索引转换为(Y,X) 数据坐标。这工作得很好,并给出了一个基本上精确到机器精度的零轮廓。然而,试图沿着这个零轮廓进行插值以获得额外的点失败了,显然是因为插值点在很小程度上偏离了真正的零轮廓;产生的错误对于预期的应用程序来说太大了。
作为一个简单的示例,可以计算peaks 函数的零轮廓,沿轮廓计算函数值,然后计算沿轮廓插值的新点处的函数值。 countours/countourc 找到的点的最大误差大约为 eps(max(Z(:)),而插值点的最大误差大约高出十个数量级。
% 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 是封闭源代码,因此无法将其用于此任务。
【问题讨论】:
标签: matlab interpolation