Placing points depending on density distribution

Hi,

Some context for my problem:

I am calculation diffusion (concentration profiles) in 1D (radial diffusion in a sphere). Because the numeric gets unstable at some time, I have to increase the lateral resolution.

I do not want to just decrease my overall dr and calculate mith that small dx and equidistant points, I just want to decrease dr, where it is needed ands keep a coarse dr where its just fine.

So I know that I would expect the largest changes during a time increment dt, when the concentration gradient is large. So for a constant dt, in this regions dx should be small, otherwise large.

I have a solution, but it lags some properties that I would benefit from.

So I just calculate thet gradient of the concentration, smooth it a little bit, normalize it to get a probability function, multipli the function with my desired number of points, and then divide my largest allowed dr by these number of points, and then take the cummulative some of all dr, to get the positions of my points.

My code looks like this:

C= some concentration profile

for iterate over time in increments dt

C(t+dt)=f(C,r);    

r_new=linspace(10^-12,r(end),10); % base increments/ largest allowed dr

grad2=abs(gradient(C)./gradient(r)); % first derivative

grad2=interp1(r,grad2,r_new); % interpolate first derivative at base increments

grad2=smooth(grad2,10); % smooth out gradient
grad2=1+round(40*(grad2/sum(grad2))); % seed 40 points weighted depending on gradient plus 1 aditional every 1/10 r ( in case gradient is zero)
grad2(isnan(grad2))=1; % in case there is no gradient at all in the previous step

diff_rnew=abs(diff(r_new)); % distance between new base points
r_new=cumsum(repelem(diff_rnew'./grad2(2:end),grad2(2:end))); % subdivide base increments depending on how many subpoints based on gradient should be placed in between
r_new(1)=[]; 
r_new(1)=10^-12;

C_new=interp1(r,C,r_new);
C=C_new;


end

So it works fine, but there are some problems that I cannot handle.

• I want a not changing number of points ( numel(r_new) ), that hat always the same value. Right now thats not the case because depending on the gradient and the round function, the number of points is just close to this value, but it always changes slightly
• Because I need to do this for multiple concentration profiles, I hope that there is a benefit when they all have the same number of points when it comes to the interpolation

So basicly:

Is there a better approach instead of my piecewise linspace point placing? Something more smooth?

Is there a way to do the interpolation more time efficient, when I have to do it like 10k times, for 10k C profiles, for 10k different r- arrays, but all C and r arrays having the same number of points?

many thanks in advance and best regards

Answer 1

[Torsten]

Use Kardos for 1d partial differential equations.

It is adaptive in space and time and available for download.

https://www.mathematik.tu-darmstadt.de/numerik-und-wissenschaftliches-rechnen/personen_numerik/details_num/jens_lang.en.jsp

[Marc Laub]

its just a side problem, that has to be implemented and run in the background of a larger matlab code