Generate a mesh with unequal steps based on a density function

I'm trying to generate a 1D mesh with unequal step length, and with a fixed number of elements [same as the initial mesh]. The length should be proportional to a node density. In the example, this density is inversely proportional to the gradient of a function. [imagine for example that you have a distribution of the temperature in a 1D mesh, and you want to make the mesh denser in the parts of the mesh where the temperature gradient is higher]

This is what I coded so far:

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

There are a few problems with this approach: 1. as you see from this example, there are big jumps when the density is very low (gradient almost zero). How could I implement a minimum/maximum time step size? 2. the node density is calculated correctly, but after "integrating" the unequal steps it can happen that the imposed large time step when the gradient is small causes to skip a high-gradient region that should have finer time-steps. [for example, please take a look at the region 1.8-1.9 in the example below, which should have small time step because it has high node density, but the large step size at ~1.75 causes to skip a large section of X]

Any suggestion to improve my code will be appreciated!

Answer 1

[ER2018]

After receiving some help from another forum [https://stackoverflow.com/questions/61248310/generate-a-mesh-with-unequal-steps-based-on-a-density-function-using-matlab], the final code that does what I want is the following:

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

Answer 2

[darova]

What about this?

x = (0:.02:2)';
y = 5*sin(2*pi*x).*x;

t = [0; cumsum(hypot(diff(x),diff(y)))];
t1 = linspace(t(1),t(end));
x1 = interp1(t,x,t1);
y1 = interp1(t,y,t1);

plot(x1,y1+1,'.-r')
hold on
plot(x,y,'.-g')
hold off

The idea: LINK