You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Error using matlab.internal.math.interp1, Sample points must be unique.
215 views (last 30 days)
Show older comments
I am getting this error when I run the attached matlab code. I couldn't figure out how to solve it and I would be very happy if you can help me.
Thanks in advance.
12 Comments
Torsten
on 24 Oct 2023
Your input file is missing.
But the error seems quite well explained: in a call to interp1, you seem to use some dupliciate values for the independent variable (in the below case x):
yq = interp1(x,y,xq)
Furkan Sencer Kaçar
on 24 Oct 2023
Firstly, thank you very much for your answer. As you said, I am trying to use some duplicate values, is it not possible?
Torsten
on 24 Oct 2023
As you said, I am trying to use some duplicate values, is it not possible?
No. To use interp1, x must be strictly monotone and unique.
Maybe you could prescribe y for non-unique x values as the arithmetic mean of the corresponding y-values ?
Star Strider
on 24 Oct 2023
dpb
on 24 Oct 2023
x — Sample points
vector
Sample points, specified as a row or column vector of real numbers. The values in x must be distinct. ...
type ddt_MatRANS
%##############################################################################
% Solves the leading order and higher order RANS-equation
% combined with the Wilcox (2006) k-omega turbulence
% model and gradient diffusion equation for suspended
% sediment. The model uses various finite difference
% approximations for calculating vertical gradients.
% This function is intended to be used by any of Matlab's
% internal ODE solvers, e.g. ode15s.
%
% Release 2 also incorporates:
% 1) Transitional k-omega model (turb=2), see Williams & Fuhrman (2016)
% 2) Sediment mixtures, see Caliskan & Fuhrman (2017)
%
% Programmed by:
% David R. Fuhrman
% Signe Schloer
% Johanna Sterner
% Ugur Caliskan
%##############################################################################
function [out]= ddt_MatRANS(t,f)
GlobalVars; % Declare global variables
% Uncomment to ramp up streaming effects at the beginning
%nramp=8; streaming = tanh(omega.*t./(nramp*2*pi)*pi); % Ramps up over nramp periods
%nramp=16; iconv = tanh(omega.*t./(nramp*2*pi)*pi); % Ramps up over nramp periods
% Update RHS evaluation counter
nrhs = nrhs + 1;
% The unknowns f=f(u; W; k; C)
u = real(f(1:ny,1));
W = real(f(1+ny:2*ny,1));
K = real(f(1+2*ny:3*ny,1));
% Eliminate any negative values
K = max(K,0); W = max(W,0);
% Reinforce bottom boundary conditions
if ikbc == 0 % k=0 boundary condition
K(1) = 0;
elseif ikbc == 1 % dk/dy=0 boundary condition
K(1) = K(2);
elseif ikbc == 2 % Dirichlet condition
K(1) = kbc;
end
u(1) = 0; % No slip condition
%u(end) = u(end-1); % Forces du/dy=0 at top boundary
% The pressure gradient term, 1/rho*dp/dx
tp = t + t0; % Phase adjusted time, to ensure ufs=0 at t=0
%t=t-500
if iwave == 0 % No wave signal
u_fs = 0; dudt_fs = 0;
elseif iwave == 1 % Second-order Stokes acceleration
u_fs = U_1m*sin(omega*tp) - U_2m*cos(2*omega*tp); % Desired free stream velocity
dudt_fs = U_1m*omega*cos(omega*tp) + 2*U_2m*omega*sin(2*omega*tp); % Desired free stream acceleration
elseif iwave == 2 % Wave form acceleration of Abreau et al. (2010)
u_fs = U_1m*omega*sqrt(1-r^2)*(sin(omega*tp) + r*sin(phi)/(1+sqrt(1-r^2)));
dudt_fs = U_1m*omega*sqrt(1-r^2)*(cos(omega*tp)-r*cos(phi)- ...
r^2/(1+sqrt(1-r^2))*sin(phi)*sin(omega*tp+phi))/(1-r*cos(omega*tp+phi))^2;
elseif iwave == 3 % Solitary wave acceleration
u_fs = U_1m*sech(omega*(t-t0)).^2; % t0=T set in main_MatRANS.m
dudt_fs = -2*U_1m*omega*sech(omega*(t-t0)).^2*tanh(omega*(t-t0));
elseif iwave == 4 % N-wave acceleration
u_fs = U_1m*1.165.*( sech(omega*(t-t0)-3*pi/4).^2 - sech(omega*(t-t0-pi/(2*omega))-3*pi/4).^2 );
dudt_fs = 2*U_1m*1.165*omega*( ...
sech(omega*(t-t0)-3*pi/4).^2*tanh(3*pi/4 - omega*(t-t0)) ...
- sech(omega*(t-t0-pi/(2*omega))-3*pi/4).^2*tanh(3*pi/4 - omega*(t-t0-pi/(2*omega))) );
%t, u_fs, dudt_fs
elseif iwave == 5 % Superimposed sech^2 signals
u_fs = 0; dudt_fs = 0; % Initialize
for i = 1:length(Umvec)
u_fs = u_fs + Umvec(i)*sech( Omegavec(i).*(t - t0 - tnvec(i)) )^2;
dudt_fs = dudt_fs - 2*Umvec(i)*Omegavec(i)*sech(Omegavec(i)*(t-t0-tnvec(i))).^2*...
tanh(Omegavec(i)*(t-t0-tnvec(i)));
end
%t, u_fs, dudt_fs
elseif iwave == 10 % User-defined time series
u_fs = interp1(t_ud,U_ud,t); % Interpolate velocity values from user-given data
dudt_fs = interp1(t_ud,Ut_ud,t); % Similarly, interpolate the acceleration
%t, u_fs, dudt_fs
end
dpdx = (streaming*u(end)/c - 1)*dudt_fs; % Pressure gradient term, 1/rho*dp/dx
dpdx = dpdx + Px; % Add steady contribution
dpdx = dpdx - iconv*S*u(end)^2/depth; % Add contribution from converging-diverging effects
% Determine velocity gradient
dy = diff(y'); % Vector of grid spacings
dudy = diff(u)./dy; % Velocity gradients
dudy(ny) = 0; % Top boundary conditions (du/dy=0)
% Transition modifications (Wilcox 2006, pp. 200--206)
if turb == 2 % Transition model on
Re_T = K./(W*nu); % Turbulence Reynolds number
%if isinf(Re_T(1)); Re_T(1) = 0; end; % Force zero value at the wall
if ikbc == 0 % Force zero value at the wall
Re_T(1) = 0;
elseif ikbc == 1; % Zero gradient
Re_T(1) = Re_T(2);
elseif ikbc == 2; % Dirichlet condition
Re_T(1) = kbc./(W(1)*nu);
%if isnan(Re_T(1)); Re_T(1) = 0; end
%Re_T(1) = 1;
end
alpha_star = (alpha_star0 + Re_T/R_k)./(1 + Re_T/R_k);
alpha = 13/25.*(alpha0 + Re_T/R_omega)./(1 + Re_T/R_omega)./alpha_star; % Modified closure coefficients
beta_star = 9/100.*(100*betaW/27 + (Re_T/R_beta).^4)./(1 + (Re_T/R_beta).^4);
%alpha_star=1; beta_star0=9/100; beta_star=9/100; alpha=13/25; % Uncomment to force std. model values
else
alpha_star = 1; beta_star0 = beta_star;
end
% Calculate omega bottom wall boundary condition
W_t = max(W, C_lim*abs(dudy)./sqrt(beta_star0./alpha_star)); % Stress-limited Omega tilde
nu_T = alpha_star.*K./W_t; % Kinematic eddy viscosity
if ikbc == 0 | ikbc == 2
tau_b = rho.*(nu)*dudy(1); % Bed shear stress
elseif ikbc == 1 % dk/dy=0 boundary condition or generalized Dirichlet condition
tau_b = rho.*(nu_T(1)+nu)*dudy(1); % Bed shear stress
end
U_f = sqrt(abs(tau_b)./rho); % Friction velocity
Shields=(tau_b)./(rho.*g.*(s-1).*d); % Shields parameter
Shields = ((d./d_50).^hiding_coef).*Shields; % Hiding / exposure effects in mixtures (van Rijn 2007)
k_Np = max(U_f*k_N/nu,1e-8); % Wall roughness
if k_Np > 5 % Intermediate or hydraulically rough
S_r = (coef1/k_Np) + ((coef2/k_Np)^2 - coef1/k_Np)*exp(5-k_Np);%-0.27*10^(-3);
else % Hydraulically smooth
S_r = (coef2/k_Np)^2;
end
W(1,1) = U_f^2*S_r/nu; % Omega at the wall boundary
W_t = max(W, C_lim*abs(dudy)./sqrt(beta_star0./alpha_star)); % Stress-limited Omega tilde
nu_T = alpha_star.*K./W_t; % Kinematic eddy viscosity
% Re-calculate bed shear stress quantities
if ikbc == 0 | ikbc == 2
tau_b = rho.*(nu)*dudy(1); % Bed shear stress
elseif ikbc == 1 % dk/dy=0 boundary condition or generalized Dirichlet condition
tau_b = rho.*(nu_T(1)+nu)*dudy(1); % Bed shear stress
end
U_f = sqrt(abs(tau_b)./rho); % Friction velocity
Shields=(tau_b)./(rho.*g.*(s-1).*d); % Shields parameter
Shields = ((d./d_50).^hiding_coef).*Shields; % Hiding / exposure effects in mixtures (van Rijn 2007)
k_Np = max(U_f*k_N/nu,1e-8); % Wall roughness
% Leading order terms in Navier-Stokes equation
dudt0(:,1) = -dpdx + nu*gradient(dudy,y') + (turb>0)*gradient(nu_T.*dudy,y');
%dudt0 = dudt0./(1-streaming.*u./c); % Adds u*du/dx term here
dudt0(1,1) = 0; % Boundary condition at the seabed (forces u=0 at the wall)
% Approximate horizontal velocity gradient
dudx = streaming.*(-dudt0(:,1)/c) + iconv.*(u.*S/depth);
% Approximate vertical velocity from local continuity equation
v = -cumtrapz(y,dudx);
%v = -cumsimps(y,dudx); % This also works
% Derivatives in K and Omega equations
dKdy = [diff(K)./dy; 0]; % Gradients
if ikbc == 1 % dKdy=0 boundary condition
dKdy(1) = 0;
end
dWdy = [diff(W)./dy; 0];
sigma_d = sigma_do.*(dKdy.*dWdy >= 0); % Closure coefficient
dKdy2 = gradient(dKdy,y'); % Second derivatives
dWdy2 = gradient(dWdy,y');
% Leading order K equation
Kw = K./W;
dKwdy = gradient(Kw,y); % Original
Two_dudy = dudy.*dudy;
dKdt0 = (nu_T.*Two_dudy - beta_star.*K.*W + nu*dKdy2 + ...
gradient(sigma_star.*alpha_star.*Kw.*dKdy,y')).*(turb>0);
% k wall boundary condition
if ikbc == 0 | ikbc == 2 % k=0 boundary condition or generalized Dirichlet condition
dKdt0(1) = 0;
elseif ikbc == 1 % dk/dy=0 boundary condition
dKdt0(1) = dKdt0(2);
end
dKdx = streaming.*(-dKdt0/c); % Approximate x derivative from leading-order time derivative
dKdx = dKdx + iconv.*(K.*S/depth); % Also add slope contribution
% Leading order Omega equation
dWdt0 = (alphaW*alpha_star.*W./W_t.*Two_dudy - betaW.*W.^2 + ...
sigma_d./W.*dKdy.*dWdy + nu*dWdy2 + ...
gradient(sigmaW.*alpha_star.*Kw.*dWdy,y')).*(turb>0);
dWdx = streaming.*(-dWdt0/c); % Approximate x derivative from leading-order time derivative
dWdx = dWdx + iconv.*(W.*S/depth); % Also add slope contribution
% Final Navier-Stokes equation
dudt = dudt0 - (u.*dudx + v.*gradient(u,y')) - 0*2/3*dKdx*(turb>0); % du/dt
dudt(1,1) = 0; % Boundary condition at the seabed (keeps u=0)
dudt(end,1) = dudt(end-1,1); % Forces du/dy=0 at top boundary
% Final K equation
dKdt = (-(u.*dKdx + v.*gradient(K,y')) + dKdt0)*(turb>0); % dK/dt
if ikbc == 0 | ikbc == 2 % k=0 boundary condition or generalized Dirichlet condition
dKdt(1) = 0; % Also forces k=0 at the bed, with k=0 initial condition
elseif ikbc == 1
dKdt(1) = dKdt(2); % dk/dy=0 boundary condition
end
% Final Omega equation
dWdt = (-(u.*dWdx + v.*gradient(W,y')) + dWdt0)*(turb>0); % dW/dt
% The C equation(s) (suspended sediment concentration)
% Calculate the reference concentration
if susp == 1 % Do if suspended sediment calculation is desired
for i=1:length(d) % Loop over individual grain sizes
% Determine critical Shields parameter, modified for beach slope
if Shields(i) >= 0 % Positive
Shields_c(i) = Shields_cU; % Uphill critical Shields parameter
else % Negative
Shields_c(i) = Shields_cD; % Downhill critical Shields parameter
end
Shields_rel(i) = abs(Shields(i))-Shields_c(i); Shields_rel(i) = max(Shields_rel(i),0);
if icb == 1 % Zyserman & Fredsoe cb
%cb(i) = 0.331*Shields_rel(i)^1.75/(1+(0.331/0.46)*Shields_rel(i)^1.75);
cb(i) = w_f(i)*0.331*Shields_rel(i)^1.75/(1+(0.331/0.32)*Shields_rel(i)^1.75); % cb_max replaced with 0.32
elseif icb == 2 % O'Donoghue & Wright cb
cb(i) = w_f(i)*0.264;
elseif icb == 3 || icb == 12 % Engelund & Fredsoe cb and Pickup function 2
p_s(i) = (1+((pi*mu_d/6)/(Shields_rel(i)+1e-16))^4)^-0.25;
if abs(Shields(i)) > Shields_c(i)+pi*p_s(i)*mu_d/6;
lambda(i) = 0.4*2*((Shields_rel(i)-pi*p_s(i)*mu_d/6)/(0.013*abs(Shields(i))*s))^0.5;
cb(i) = w_f(i)*0.6/(1+1/lambda(i))^3;
else
lambda(i) = 0;
cb(i) = 0;
end
elseif icb == 4 % % Einstein's formula for the reference concentration
p_s(i) = (1+((pi*mu_d/6)./(Shields_rel(i)+1e-16)).^4).^-0.25;
cb(i) = w_f(i)*pi*p_s(i)/12;
elseif icb == 5 % Brors
cmax = 0.3; theta_crs = 0.25; theta_sheet = 0.75; % Brors values
cb(i) = w_f(i).*cmax.*(abs(Shields(i))-theta_crs)./(theta_sheet-theta_crs);
cb(i) = cb(i).*(abs(Shields(i))>theta_crs);
cb(i) = min(cb(i),cmax);
elseif icb == 6 % Tanh
theta_susp = 0.25; cmax = 0.264; theta_sf = 2.0;
cb(i) = w_f(i).*cmax.*tanh(pi./theta_sf.*(abs(Shields(i))-theta_susp));
cb(i) = cb(i).*(abs(Shields(i))>theta_susp);
elseif icb == 11 % Pickup function of van Rijn
if abs(Shields(i)) > Shields_c(i)
gammaS(i) = 0.00083*(abs(Shields(i))/Shields_c(i)-1)^1.5;
pickup(i) = w_s(i)*gammaS(i);
else
pickup(i) = 0;
end
end
cb(i) = cb(i)*(((2*d(i))/b)^(ref_conc_coef));
if icb == 12 % Pickup function 2
if abs(Shields(i)) > Shields_c(i)
gammaS(i) = cb(i);
pickup(i) = w_s0(i)*gammaS(i);
else
pickup(i) = 0;
end
end
end
C_total = zeros(nyc,1);
for i=1:length(d)
C = f(1+3*ny+nyc*(i-1):3*ny+nyc*i);
C(1) = cb(i);
C = max(C,0);
C_total = C_total + C;
end
for i=1:length(d)
C = f(1+3*ny+nyc*(i-1):3*ny+nyc*i);
C(1) = cb(i);
C = max(C,0);
if iextrap && icb < 10
cextrap = C(2) + (C(3)-C(2))/(yc(3)-yc(2))*(y(1)-y(2));
C(1) = max(cb(i),cextrap);
end
% Sediment diffusivity
if beta_s == 0 % Use van Rijn (1984) correction
eps_s = (1 + 2.0.*(w_s0/U_f)^2).*nu_T(end-nyc+1:end);
else
eps_s = beta_s.*nu_T(end-nyc+1:end); % Sediment diffusivity
end
eps_s = eps_s + nu; % Add moledular diffusion
dCdy = [diff(C)./diff(yc'); 0]; % dc/dy
if icb > 10 % Gradient condition
dCdy(1) = -pickup/eps_s(1);
end
% Settling term
if Hind_set == 1 % Hindered settling
%w_s(1:nyc,i) = w_s0(:,i).*(1-C_total).^nhs(i); % w_s varies with concentration
w_s = w_s0(i).*(1-C_total).^nhs(i); % w_s varies with concentration
settling = gradient(w_s.*C,yc');
else % Constant w_s
settling = w_s0(i).*dCdy; % w_s*dc/dy
end
% Leading order C equation
diffusion = gradient(eps_s.*dCdy,yc');
dCdt0 = settling + diffusion;
dCdx = streaming.*(-dCdt0./c);
% Final C equation
dCdt(1+(i-1)*nyc:i*nyc,1) = -(u(end-nyc+1:end).*dCdx + v(end-nyc+1:end).*gradient(C,yc')) + dCdt0; % Full dC/dt
end
else
dCdt(1:length(d)*nyc,1) = zeros(nyc,1); % Otherwise, just set vector to zero
end
% Turbulence supression due to density gradients
if turb && iturbsup && susp
if iturbsup == 1
rho_m = rho.*s.*C_total + rho.*(1-C_total); % Density of fluid-sediment mixture
rhomax = n*rho + (1-n)*rho*s; % Maximum fluid-sediment density
drdy_b = (rho_m(2)-rho_m(1))./(yc(2)-yc(1)); % d rho_m/dy @ y=b
rho_m = [rho_m(1) + drdy_b.*(y(1:ib-1)' - y(ib)); rho_m]; % Extrapolate fluid-sediment density
rho_m = min(rho_m,rhomax); rho_m = max(rho_m,rho); % Bound values
drho_mdy = gradient(rho_m,y');
N2 = -g./rho_m.*drho_mdy; % Square of Brunt-Vaisala frequency [1/s^2]
elseif iturbsup == 2
dCdy_b = (C_total(2)-C_total(1))./(yc(2)-yc(1)); % dCdy @ y=b
C_m = [C_total(1) + dCdy_b.*(y(1:ib-1)' - y(ib)); C_total]; % Extrapolate concentration
C_m = min(C_m,1-n); C_m = max(C_m,0); % Bound values
dC_mdz = gradient(C_m,y'); % dC/dy
N2 = -g.*(s-1).*dC_mdz; % Leading-order contribution
end
B = N2.*nu_T./sigma_p; % Buoyancy flux [m^2/s^3]
dKdt = dKdt - B; % Modify k-equation
c3eps = N2<0; % 1 for N^2<0, 0 for N^2>=0 (Reussink et al. 2009)
dWdt = dWdt - c3eps.*N2; % Modify omega-equation
% Enforce K boundary condition
if ikbc == 0 | ikbc == 2 % k=0 boundary condition or generalized Dirichlet condition
dKdt(1) = 0; % Also forces k=0 at the bed, with k=0 initial condition
elseif ikbc == 1
dKdt(1) = dKdt(2); % dk/dy=0 boundary condition
end
end
% Filter du/dt
if filt_U > 0
filt = [filt_U 1-2.*filt_U filt_U]; % Filter
dudt = filter(filt,1,dudt); % Initial filtering
dudt = [dudt(2:end); dudt(end)]; % Eliminate phase shift
dudt(1) = 0; % Enforce boundary condition
dudt(end) = dudt(end-1); % du/dy=0 at top
end
% Filter dK/dt
if filt_K > 0 && turb
filt = [filt_K 1-2.*filt_K filt_K]; % Filter
dKdt = filter(filt,1,dKdt); % Initial filtering
dKdt = [dKdt(2:end); dKdt(end)]; % Eliminate phase shift
% Enforce K boundary condition
if ikbc == 0 | ikbc == 2 % k=0 boundary condition or generalized Dirichlet condition
dKdt(1) = 0; % Also forces k=0 at the bed, with k=0 initial condition
elseif ikbc == 1
dKdt(1) = dKdt(2); % dk/dy=0 boundary condition
end
end
% Filter d(omega)/dt
if filt_W > 0 && turb
filt = [filt_W 1-2.*filt_W filt_W]; % Filter
dWdt = filter(filt,1,dWdt); % Initial filtering
dWdt = [dWdt(2:end); dWdt(end)]; % Eliminate phase shift
end
% Filter dc/dt
if filt_C > 0 && susp
filt = [filt_C 1-2.*filt_C filt_C]; % Filter
dCdt = filter(filt,1,dCdt); % Initial filtering
dCdt = [dCdt(2:end); dCdt(end)]; % Eliminate phase shift
dCdt(end) = 0; % Forces C=0 at top boundary condition
end
% Create output vector containing all time derivatives
out = [dudt; dWdt; dKdt; dCdt];
% Display time to screen every 500th stage evaluation
if mod(nrhs,noutput) == 0
disp([ 't = ' num2str(t) ' s, k_N^+ = ' num2str(k_Np) ...
', Shields = ' num2str(Shields)...
', CPU time = ' num2str(toc) ' s, cb = ' num2str(cb)]) % Display time to screen
if ioutplot % Also show plots, if desired (to turn off set to 'if 0')
figure(100);
ym = h_m;
% Plot physical variables
subplot(3,4,1); plot(u,y); xlabel('u (m/s)'); ylabel('y (m)'); ylim([0 ym]);
subplot(3,4,2); plot(K,y); xlabel('k (m^2/s^2)'); ylim([0 ym]);
subplot(3,4,3); plot(W,y); xlabel('\omega (1/s)'); ylim([0 ym]);
if susp; subplot(3,4,4); plot(C,yc); xlabel('C'); ylim([0 ym]); end
% Plot time derivatives
subplot(3,4,5); plot(dudt,y); xlabel('du/dt (m/s^2)'); ylabel('y (m)'); ylim([0 ym]);
subplot(3,4,6); plot(dKdt,y); xlabel('dk/dt (m^2/s^3)'); ylim([0 ym]);
subplot(3,4,7); plot(dWdt,y); xlabel('d\omega/dt (1/s^2)'); ylim([0 ym]);
if susp; subplot(3,4,8); plot(dCdt,yc); xlabel('dC/dt (1/s)'); ylim([0 ym]); end;
subplot(3,1,3); hold on; plot(t,u(end),'b.'); hold off;
%hold on; plot(t,u_fs,'r.'); hold off; % Add analytic value
xlabel('t (s)'); ylabel('u_0(t) (m/s)'); box on; grid on;
if turb == 2; % Transitional model
disp(['max(Re_T) = ' num2str(max(Re_T))]);
%figure(200); plot(alpha_star,y,'b-o'); xlim([0 1]);
%xlabel('\alpha^*'); ylabel('y (m)');
end;
% Update plot
drawnow;
end
end
which it wasn't but shows trying to interpolate over a time vector. You'll have to use distinct time points to use the function as it was written--not knowing what it is you're trying to do as to why you would have duplicated times, the closest thing you could do would be to introduce a slight difference between points that hopefully wouldn't otherwise effect the results by introducing huge gradients if, for example, the corresponding y values are grossly different as in a step change; the model is likely not designed to handle such transients anyway.
Furkan Sencer Kaçar
on 24 Oct 2023
Edited: Torsten
on 24 Oct 2023
@dpb These are the input values I've used, which seem to me distinct
data = readmatrix("Inputs.txt",'Delimiter',{',','='})
data = 4004×2
0 0.0500
0.1000 0.0500
0.2000 0.0500
0.3000 0.0500
0.4000 0.0500
0.5000 0.0500
0.6000 0.0500
0.7000 0.0500
0.8000 0.0500
0.9000 0.0500
size(data(:,1))
ans = 1×2
4004 1
size(unique(data(:,1)))
ans = 1×2
4001 1
idx = find(diff(data(:,1))==0)
idx = 3×1
401
2002
2403
data(idx,1)
ans = 3×1
40
200
240
Stephen23
on 24 Oct 2023
Edited: Stephen23
on 24 Oct 2023
M = readmatrix('Inputs.txt')
M = 4004×2
0 0.0500
0.1000 0.0500
0.2000 0.0500
0.3000 0.0500
0.4000 0.0500
0.5000 0.0500
0.6000 0.0500
0.7000 0.0500
0.8000 0.0500
0.9000 0.0500
numel(unique(M(:,1)))
ans = 4001
numel(unique(M(:,2)))
ans = 1602
The first column has 4001 unique values out of 4004 rows. So not distinct.
dpb
on 24 Oct 2023
Edited: dpb
on 24 Oct 2023
data=readmatrix('Inputs.txt');
[u,iu]=unique(data(:,1));
subplot(2,2,[1,2])
plot(u,data(iu,2),'k-')
hold on
ix=find(diff(data(:,1))==0);
plot(data(ix,1),data(ix,2),'r*')
ylim([0 0.6])
subplot(2,2,3)
plot(u,data(iu,2),'k-')
hold on
plot(data(ix(1),1),data(ix(1),2),'r*')
xlim(data(ix(1))+[-5 +5])
ylim([0.5455 0.5505])
subplot(2,2,4)
plot(u,data(iu,2),'k-')
hold on
plot(data(ix(3),1),data(ix(3),2),'r*')
xlim(data(ix(3))+[-5 +5])
ylim([0.5455 0.5505])
shows the trace is smooth and that the duplicated points could be safely discarded...looks like maybe somebody/some prior function augmented the original response with the peak locations and the end point of the first transient with a pretrigger interval before the second...I didn't check to see, but maybe the second is a replication of the first????
Steven Lord
on 24 Oct 2023
FYI if users want to mark the peak values, one way to do that without duplicating the peak values (assuming you ONLY want those peaks marked) is to use the MarkerIndices property of line objects.
x = 1:10;
y = x.^2;
plot(x, y, 'o-', 'MarkerIndices', 1:2:10) % Mark every other point
xline(1:2:10, ':')
You can use islocalmax to identify the maximum points.
figure
x = 0:720;
y = sind(x);
loc = islocalmax(y);
plot(x, y, 'o-', 'MarkerIndices', find(loc))
dpb
on 24 Oct 2023
Edited: dpb
on 25 Oct 2023
Good sidebar comment @Steven Lord -- I wasn't aware that 'MarkerIndices' had been added to line properties....or maybe it's always been there and I just didn't recall it; I know I've combed the list in its entirety before.
Answers (2)
Rik
on 24 Oct 2023
Edited: Rik
on 24 Oct 2023
interp1 is not intended to fit equations. There are other tools to do that. If you insist on using this function, you will either have to calculate some consensus value for each unique x (e.g. the mean, the median, the maximum, the minimum, ...), or you should adjust the x values slightly so each value becomes unique.
Note that the order of your x-values suddenly matters for the second option, because each element will have a different offset based on the position.
x = [1 1 1 2 3];
y = [0 1 2 3 4];
xq = 1:3;
[x2,y2]=AdjustWithMedian(x,y)
x2 = 1×3
1 2 3
y2 = 1×3
1 3 4
[x3,y3]=AdjustWithEps(x,y)
x3 = 1×5
1.0000 1.0000 1.0000 2.0000 3.0000
y3 = 1×5
0 1 2 3 4
[x4,y4]=AdjustWithPolyfit(x,y)
x4 = 1×3
1 2 3
y4 = 1×3
1.0625 2.6250 4.1875
plot(x,y,'b*','DisplayName','raw data'),axis([0.5 3.5 -0.5 4.5])
hold on
plot(xq,interp1(x2,y2,xq),'m--o','DisplayName','Adjusted with median')
plot(xq,interp1(x3,y3,xq),'k:o','DisplayName','Adjusted with eps')
plot(xq,interp1(x4,y4,xq),'g-.o','DisplayName','Adjusted with polyfit')
hold off
legend('Location','SouthEast')
function [x_adj,y]=AdjustWithEps(x,y)
% Add a different amount to each x-value to make them unique.
% The sum of the shift should be equal to 0 to avoid an overall drift.
shift = 1:numel(x);
shift = shift - mean(shift);
x_adj = x(:) + shift(:)*eps;
x_adj = reshape(x_adj,size(x));
end
function [new_x,new_y]=AdjustWithMedian(x,y)
% For each unique x, only keep the median y-value
[new_x,ignore,ind] = unique(x);
new_y = accumarray(ind,y,[numel(new_x) 1],@median);
if isrow(x),new_y = new_y.';end % Keep directions consistent
end
function [x,y]=AdjustWithPolyfit(x,y)
p=polyfit(x,y,1);
x=unique(x);
y=polyval(p,x);
end
dpb
on 24 Oct 2023
Edited: dpb
on 24 Oct 2023
As the follow on comment above shows, you can safely ignore the duplicate points and retain the original function -- use
[~,iu]=unique(data(:,1));
data=data(iu,:);
in place of the full data array to reduce it to the set of unique times and associated values; then pass the resulting time, response vectors to the function instead and all should work just fine...although there might be a question regarding there being two transients in the overall trace? Or is this an input concentration to be modelled its dispersion, maybe, in which case it could be intended and correct--we don't know anything about the problem trying to solve, just that there were duplicated times in this input trace that aren't allowed by the construction of the function.
See Also
Categories
Find more on Interactive Control and Callbacks in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)