Error using matlab.int​ernal.math​.interp1, Sample points must be unique.

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.

[Torsten]

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]

Firstly, thank you very much for your answer. As you said, I am trying to use some duplicate values, is it not possible?

[Torsten]

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]

Another option is to use polyfit and polyval. This is much less efficient, however it may be an acceptable alternative.

[dpb]

Not with interp1, no. The doc says explicitly for input x

x — Sample points

vector

Sample points, specified as a row or column vector of real numbers. The values in x must be distinct. ...

As @Torsten says, you didn't attach the data unless it's in the function iteslf...

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]

@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]

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.

[Torsten]

@Furkan Sencer Kaçar

No, they are not distinct (see above).

[dpb]

Carrying on from @Torsten and @Star Strider on the specifics of the input data set,

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]

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]

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.

[Steven Lord]

It was introduced in release R2016b.

Answer 1

[Rik]

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

Answer 2

[dpb]

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.