使用 Matlab 基于密度函数生成步长不等的网格
Generate a mesh with unequal steps based on a density function, using Matlab
我正在尝试生成步长不等且元素数量固定的一维网格 [与初始网格相同]。
长度应与节点密度成正比。在示例中,此密度与函数的梯度成反比。
[例如,假设您在一维网格中有温度分布,并且您希望在温度梯度较高的网格部分使网格更密集]
这是我目前编写的代码:
% % % 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. 正如你从这个例子中看到的,当密度很低(梯度几乎为零)时会有很大的跳跃。我如何实现 minimum/maximum 时间步长?
2. 节点密度计算正确,但在 "integrating" 步长不等之后,可能会发生梯度较小时强加的大时间步长导致跳过本应具有更精细时间步长的高梯度区域。 [例如,请看下面示例中的区域 1.8-1.9,它应该有小的时间步长,因为它具有高节点密度,但是 ~1.75 的大步长导致跳过 X 的大部分]
任何改进我的代码的建议都将不胜感激!
计算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显示结果。
rahnema1的回答给了我很大的帮助,但还有两个问题:
1- 新网格的第一个元素与原始网格的第一个元素不相同
2- 如果梯度在某个点为零,interp1 函数将给出错误 ["The grid vectors are not strictly monotonic increasing."]
对于第一点,我将定义 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% 的答案[可能我没有很好地解释原始请求]!
我正在尝试生成步长不等且元素数量固定的一维网格 [与初始网格相同]。 长度应与节点密度成正比。在示例中,此密度与函数的梯度成反比。 [例如,假设您在一维网格中有温度分布,并且您希望在温度梯度较高的网格部分使网格更密集]
这是我目前编写的代码:
% % % 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. 正如你从这个例子中看到的,当密度很低(梯度几乎为零)时会有很大的跳跃。我如何实现 minimum/maximum 时间步长? 2. 节点密度计算正确,但在 "integrating" 步长不等之后,可能会发生梯度较小时强加的大时间步长导致跳过本应具有更精细时间步长的高梯度区域。 [例如,请看下面示例中的区域 1.8-1.9,它应该有小的时间步长,因为它具有高节点密度,但是 ~1.75 的大步长导致跳过 X 的大部分]
任何改进我的代码的建议都将不胜感激!
计算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显示结果。
rahnema1的回答给了我很大的帮助,但还有两个问题: 1- 新网格的第一个元素与原始网格的第一个元素不相同 2- 如果梯度在某个点为零,interp1 函数将给出错误 ["The grid vectors are not strictly monotonic increasing."]
对于第一点,我将定义 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% 的答案[可能我没有很好地解释原始请求]!