使用 Matlab 基于密度函数生成具有不等步长的网格

Question

我正在尝试生成一个步长不等的一维网格，并且具有固定数量的元素 [与初始网格相同]。 长度应与节点密度成正比。在示例中，该密度与函数的梯度成反比。 [例如，假设您在一维网格中具有温度分布，并且您希望在温度梯度较高的网格部分使网格更密集]

这是我目前编写的代码：

% % % Initial fixed-step 1D mesh
X=(0:.01:2)';
% % % Values of a function at each mesh node [in this example,  f(x)=5*sin(2*pi*x)*x ]
Y=5*sin(2*pi*X).*X;

% % % Calculate density of mesh points based on the Y function gradient
rho=[1e-9; abs(diff(Y))];

% % % Calculate x-steps from the original mesh
h = diff(X);
% % % Rescale the steps based on the inverse of the density
F = cumsum([0; h]./rho);
% % % Make sure F is between 0 and 1
F = F/F(end);
% % % Calculate the new mesh with scaled steps
X2 = X(1) + F * (X(end)-X(1));
% % % Interpolate the function Y at the new positions
Y2 = interp1(X,Y,X2);

% % % Plot
figure
subplot(2,1,1)
hold on
plot(X,Y,'ko-')
plot(X2,Y2,'r.-')
xlabel('x')
ylabel('y')
subplot(2,1,2)
plot(X,rho,'ko-')
xlabel('x')
ylabel('rho')

这种方法存在一些问题： 1. 正如你从这个例子中看到的，当密度非常低（梯度几乎为零）时会有很大的跳跃。如何实现最小/最大时间步长？ 2. 节点密度计算正确，但在“积分”不等步之后，可能会发生在梯度较小时施加的大时间步长会导致跳过应该具有更精细时间步长的高梯度区域。 [例如，请看下面示例中的区域 1.8-1.9，它应该具有较小的时间步长，因为它具有高节点密度，但 ~1.75 的大步长会导致跳过很大一部分 X]

任何改进我的代码的建议将不胜感激！

Answer 1

计算rho 的累积和 (CDF)。从CDF 中抽取等距样本。从CDF 映射到X 以获取新的X3。从X3映射到Y得到Y3：

CDF = cumsum(rho);
eq_smpl = linspace(CDF(1), CDF(end), numel(CDF)+1).';
eq_smpl = eq_smpl(1:end-1) + diff(eq_smpl)/2; %use center points
X3 = interp1(CDF, X, eq_smpl);
Y3 = interp1(X, Y, X3);

plot(X3,Y3,'ro-')
hold on
plot(X,Y,'k.')

第三个subplot 显示结果。

Answer 2

rahnema1的回答给了我很大的帮助，但还有两个问题： 1-新网格的第一个元素与原始网格的第一个元素不相同 2-如果梯度在某个点为零，interp1函数将给出错误[“网格向量不是严格单调递增的。”]

对于第一点，我将定义 eq_smpl 的两行替换为以下行：

eq_smpl = linspace(CDF(1), CDF(end), numel(CDF))';

[取与CDF一样多的元素，而不是使点居中]

对于第二点，我在计算 rho 之后添加了一行，以用小的非零值替换最终的 0：

rho(rho==0)=1e-12;

完成我想要的最终代码如下：

% % % Initial fixed-step 1D mesh
X=(0:.01:2)';
% % % Values of a function at each mesh node [in this example,  f(x)=5*sin(2*pi*x)*x ]
Y=5*sin(2*pi*X).*X;

% % % Calculate density of mesh points based on the Y function gradient
rho=[0; abs(diff(Y)./abs(diff(X)))];
% % % Replace eventual 0 with small non-zero values
rho(rho==0)=1e-12;

CDF = cumsum(rho);
eq_smpl = linspace(CDF(1), CDF(end), numel(CDF))';

% % % Calculate new mesh
X3 = interp1(CDF, X, eq_smpl);
% % % Interpolate the function Y at the new positions
Y3 = interp1(X, Y, X3);

% % % Plot
figure
subplot(2,1,1)
hold on
plot(X,Y,'ko-')
plot(X3,Y3,'r.-')
xlabel('x')
ylabel('y')
subplot(2,1,2)
plot(X,rho,'ko-')
xlabel('x')
ylabel('rho')

再次感谢 rahnema1 提供了 90% 的答案[可能我没有很好地解释最初的请求]！