Генерация me sh с неравными шагами на основе функции плотности с использованием Matlab

Question

Я пытаюсь сгенерировать 1D me sh с неравной длиной шага и с фиксированным числом элементов [такое же, как в начальной сетке]. Длина должна быть пропорциональна плотности узла. В этом примере эта плотность обратно пропорциональна градиенту функции. [представьте, например, что у вас есть распределение температуры в 1D me sh, и вы хотите сделать me sh более плотным в тех частях me sh, где градиент температуры выше]

Это то, что я до сих пор кодировал:

% % % 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

Рассчитать совокупную сумму (CDF) rho. Возьмите с равными интервалами образцы из 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.')

Третий субплот показывает результат.

Answer 2

Ответ rahnema1 оказал мне огромную помощь, но остались еще две проблемы: 1 - первый элемент нового me sh не идентичен первому элементу оригинальной сетки 2 - в случае, если градиент равен нулю в какой-то момент функция interp1 выдаст ошибку [«Векторы сетки не являются строго монотонными c увеличивается.»]

Для 1-й точки я заменил две строки, определяющие eq_smpl, следующей строкой:

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

[беря столько элементов, сколько в CDF, и не центрируя точки]

Для 2-й точки я добавил строку после вычисления 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% ответа [вероятно, я не сделал ' очень хорошо объясняет исходный запрос]!