How can I use Interp1 here instead of ismembertol?

Hi! Essentially, I am trying to find the distance away from the origin such that the concentration is an exact match to a certain value, and thus requiring interpolation. I am currently using ismembertol as a work around because I've come into some issues with Interp1 (which I'm fairly sure is what I need to be using), but it simply isn't accurate enough for my purposes.

Using interp1 gives an error since the dataset does not consist of unique values, but I cannot just condense it to unique values because I have to use the mesh I've generated to find the distance from the origin (unless I'm wrong here of course).

I am wondering:

• Can I somehow doctor the dataset to retain its order and indicies while removing repeated datapoints?

I.e [1, 5, 7, 8, 8, 16, 35, 35, 35, 60, 120] might become [1, 5, 7, 8, NaN, 16, 35, NaN, NaN, 60, 120]

• Can I use a different method of indexing that would work better in this case?
• Is there another method entirely that would be better suited to this?

It may be important to note that the next step is to vary V and recieve an array of dc values to plot a dc-V curve - so the solution would need to not cause any problems here.

Any help is greatly appriciated, I'm very fair from the code wizard i need ot be for this...

Here's my code so far:

% process Parameters
V = 10 % S/L interface velocity in microns/sec

% Material Parameters
c0 = 0.2; % Alloy composition
k = 0.28; % Partition Coefficient (conc solute in solid/ liquid)
D = 10; % Diffusivity in microns^2/sec
dc_tol = 1.01*c0; % Defines how close to c0 is acceptable

% Mesh Parameters
n = 500; % number of nodes
L = 5; % Length of nodes in microns
dx = L / n; % defining node spacing
x = linspace(0, L, n); % Mesh

% Initialise concentration field
c_field = ones(1, n)*c0; % Creates a field vector of intial concentrration
c_field_update = c_field;

% Initialise dc Field
dc_field = zeros(1, 100);

% Boundary Conditions
c_field(1) = c0 / k;
c_field(n) = c0;

% Solver controls
tol = 1e-3; % Sets the tolerance for stopping the iteration
err = 1; % sets starting value to check against tol

F0 = 0.5*dx*D/V; % Finite Difference Coefficient

% Solver
while err > tol

    % Apply Boundary Conditions
    c_field_update(1) = c_field(1);
    c_field_update(n) = c_field(n);

    % Finite Difference Scheme 
    for i = 2:(n-1)
        c_field_update(i) = 0.5*(1+F0)*c_field(i+1)+0.5*(1-F0)*c_field(i-1);
    end 

    % Compute Error
    errVector = abs((c_field_update - c_field)./ c_field_update);
    err = max(errVector);

    % Update Concentration Field
    c_field = c_field_update;

end 

% Determine Solute Boundary Layer
% The ismembertol function gives number of nodes into the mesh of the first value that meets the dc_tol within a certain toleranmce, idx_tol)

idx_tol = 0.001*dc_tol;
[~, idx] = ismembertol(dc_tol, c_field, idx_tol);
dc = idx*dx % solute boundary layer thickness in microns
dcApprox = 2*D/V;


% Visualisation
figure;
plot (x, c_field, 'r-', 'lineWidth', 2);
xlabel('Position (microns)', "FontSize", 20);
ylabel('Concentration', "FontSize", 20);
title('Concentration Profile', "FontSize", 20);

Answer 1

[Star Strider]

That might be possible. The interp1 function requires two vectors of equal size, and then returns the interpolated value of the second argument vector that corresponds to the value in the first argument vector.

The ‘c_field’ variable is a vector. It needs a second vector and a value to be interpoalted (scalar or vector) to perform the interpolation.

I created something similar to what you may want, however I need to know if it does what you want —

% process Parameters
V = 10 % S/L interface velocity in microns/sec

V = 10

% Material Parameters
c0 = 0.2; % Alloy composition
k = 0.28; % Partition Coefficient (conc solute in solid/ liquid)
D = 10; % Diffusivity in microns^2/sec
dc_tol = 1.01*c0; % Defines how close to c0 is acceptable

% Mesh Parameters
n = 500; % number of nodes
L = 5; % Length of nodes in microns
dx = L / n; % defining node spacing
x = linspace(0, L, n); % Mesh

% Initialise concentration field
c_field = ones(1, n)*c0; % Creates a field vector of intial concentrration
c_field_update = c_field;

% Initialise dc Field
dc_field = zeros(1, 100);

% Boundary Conditions
c_field(1) = c0 / k;
c_field(n) = c0;

% Solver controls
tol = 1e-3; % Sets the tolerance for stopping the iteration
err = 1; % sets starting value to check against tol

F0 = 0.5*dx*D/V; % Finite Difference Coefficient

% Solver
while err > tol

    % Apply Boundary Conditions
    c_field_update(1) = c_field(1);
    c_field_update(n) = c_field(n);

    % Finite Difference Scheme 
    for i = 2:(n-1)
        c_field_update(i) = 0.5*(1+F0)*c_field(i+1)+0.5*(1-F0)*c_field(i-1);
    end 

    % Compute Error
    errVector = abs((c_field_update - c_field)./ c_field_update);
    err = max(errVector);

    % Update Concentration Field
    c_field = c_field_update;

end 

% Determine Solute Boundary Layer
% The ismembertol function gives number of nodes into the mesh of the first value that meets the dc_tol within a certain toleranmce, idx_tol)

idx_tol = 0.001*dc_tol;
[~, idx] = ismembertol(dc_tol, c_field, idx_tol);
dc = idx*dx % solute boundary layer thickness in microns

dc = 0.7300

dx

dx = 0.0100

dc_tol

dc_tol = 0.2020

c_field

c_field = 1×500

0.7143 0.6959 0.6778 0.6598 0.6422 0.6247 0.6076 0.5907 0.5743 0.5580 0.5423 0.5267 0.5117 0.4968 0.4826 0.4686 0.4552 0.4419 0.4294 0.4170 0.4053 0.3937 0.3829 0.3722 0.3622 0.3523 0.3431 0.3341 0.3257 0.3175

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[cfmin,cfmax] = bounds(c_field)

cfmin = 0.2000

cfmax = 0.7143

c_fieldu = unique(c_field, 'stable')                            % ‘interp1’ Requires Unique Values

c_fieldu = 1×206

0.7143 0.6959 0.6778 0.6598 0.6422 0.6247 0.6076 0.5907 0.5743 0.5580 0.5423 0.5267 0.5117 0.4968 0.4826 0.4686 0.4552 0.4419 0.4294 0.4170 0.4053 0.3937 0.3829 0.3722 0.3622 0.3523 0.3431 0.3341 0.3257 0.3175

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

idx = interp1(c_fieldu, (1:numel(c_fieldu)), dc_tol)            % ‘interp1’ Version

idx = 73.1043

dc = idx*dx

dc = 0.7310

dcApprox = 2*D/V;


% Visualisation
figure;
plot (x, c_field, 'r-', 'lineWidth', 2);
xlabel('Position (microns)', "FontSize", 20);
ylabel('Concentration', "FontSize", 20);
title('Concentration Profile', "FontSize", 20);

EDIT — (5 Apr 2025 at 16:06)

‘Can I somehow doctor the dataset to retain its order and indicies while removing repeated datapoints?’

I just used the unique function for this. Is there a reason to retain all the original data?

.

[Charlie Rideout]

Ahhhhh, this makes sense! It's exactly what I'm wanting!

I think I made two erroneous assumptions:

• I took the value at 4sf to be the exact value, and thought there were more duplicate data points than there are.
• Since the duplicates are only at the end of the dataset, the unique function does not change the position of the required values within in the vector - as I though it would.

[Star Strider]

As always, my pleasure!

Answer 2

[Stephen23]

Fiddling around with indices and ISMEMBERTOL are red-herrings and misdirections from your actual task. Essentially what you are trying to do is root-finding, so use FZERO (or similar). Using FZERO is conceptually neater and much more accurate. I also made a few other improvements (e.g. vectorizing the code inside the WHILE loop).

% Process parameters:
V   = 10;        % S/L interface velocity (microns/sec)

% Material parameters:
c0  = 0.2;       % Alloy composition
k   = 0.28;      % Partition coefficient (solute in solid/liquid)
D   = 10;        % Diffusivity (microns^2/sec)
dct = 1.01 * c0; % Threshold concentration value

% Mesh parameters:
n   = 500;       % Number of nodes
L   = 5;         % Length of domain (microns)
dx  = L / n;     % Node spacing
x   = linspace(0, L, n); % Spatial mesh

% Initialize concentration field:
cfd = ones(1, n) .* c0;  
cfd(1) = c0 / k;  
cfd(n) = c0;      

% Solver controls:
tol = 1e-3;
err = 1;
f0  = 0.5 * dx * D / V;

% Finite difference solver (partially vectorized):
while err > tol
    cfo = cfd;
    %
    cfd(2:n-1) = 0.5 * (1+f0) .* cfd(3:n) + 0.5 * (1-f0) .* cfd(1:n-2);
    %
    nzi = (cfd~=0);
    rel = zeros(size(cfd));  
    rel(nzi)  = abs((cfd(nzi) - cfo(nzi)) ./ cfd(nzi));
    rel(~nzi) = abs(cfd(~nzi) - cfo(~nzi));
    err = max(rel);
end

% Determine solute boundary layer:
cfu = @(xi) interp1(x, cfd, xi, 'pchip') - dct;
dc  = fzero(cfu, x([1,n]));  
dca = 2 * D / V;

% Display results:
fprintf('Computed solute boundary layer thickness: %.4f microns\n', dc);

Computed solute boundary layer thickness: 0.7225 microns

fprintf('Approximate boundary layer thickness: %.4f microns\n', dca);

Approximate boundary layer thickness: 2.0000 microns

% Visualization:
figure;
plot(x, cfd,'r-', dc,dct,'b*', 'LineWidth', 2);
xlabel('Position (microns)', "FontSize", 20);
ylabel('Concentration', "FontSize", 20);
title('Concentration profile', "FontSize", 20);
grid on;

●

Lets zoom in on that point to see how it matches the plotted data:

plot(x, cfd,'r-', dc,dct,'b*', 'LineWidth', 2);
xlim([0.999,1.001]*dc)
ylim([0.999,1.001]*dct)

[Stephen23]

Compared against your approach using ISMEMBERTOL:

% process Parameters
V = 10; % S/L interface velocity in microns/sec

% Material Parameters
c0 = 0.2; % Alloy composition
k = 0.28; % Partition Coefficient (conc solute in solid/ liquid)
D = 10; % Diffusivity in microns^2/sec
dc_tol = 1.01*c0 % Defines how close to c0 is acceptable

dc_tol = 0.2020

% Mesh Parameters
n = 500; % number of nodes
L = 5; % Length of nodes in microns
dx = L / n; % defining node spacing
x = linspace(0, L, n); % Mesh

% Initialise concentration field
c_field = ones(1, n)*c0; % Creates a field vector of intial concentrration
c_field_update = c_field;

% Initialise dc Field
dc_field = zeros(1, 100);

% Boundary Conditions
c_field(1) = c0 / k;
c_field(n) = c0;

% Solver controls
tol = 1e-3; % Sets the tolerance for stopping the iteration
err = 1; % sets starting value to check against tol

F0 = 0.5*dx*D/V; % Finite Difference Coefficient

% Solver
while err > tol
    
    % Apply Boundary Conditions
    c_field_update(1) = c_field(1);
    c_field_update(n) = c_field(n);
    
    % Finite Difference Scheme 
    for i = 2:(n-1)
        c_field_update(i) = 0.5*(1+F0)*c_field(i+1)+0.5*(1-F0)*c_field(i-1);
    end 
    
    % Compute Error
    errVector = abs((c_field_update - c_field)./ c_field_update);
    err = max(errVector);
    
    % Update Concentration Field
    c_field = c_field_update;
    
end 

% Determine Solute Boundary Layer
c_fieldu = unique(c_field, 'stable');
idx = interp1(c_fieldu, (1:numel(c_fieldu)), dc_tol);
dc = idx*dx

dc = 0.7310

plot(x, c_field,'r-', dc,dc_tol,'b*', 'LineWidth', 2);
xlim([0.999,1.001]*dc)
ylim([0.999,1.001]*dc_tol)