How can I use Interp1 here instead of ismembertol?

2 次查看（过去 30 天）
Charlie Rideout
Charlie Rideout 2025-3-6
评论： Stephen23 2025-3-8
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);

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2025-3-6
编辑：Star Strider 2025-3-6
Ran in:
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?
.
  2 个评论
Charlie Rideout
Charlie Rideout 2025-3-7
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.

请先登录，再进行评论。

更多回答（1 个）

Stephen23
Stephen23 2025-3-7
编辑：Stephen23 2025-3-7
Ran in:
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)
  1 个评论
Stephen23
Stephen23 2025-3-8
Ran in:
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)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Dates and Time 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by