I want to find particular values of t at which y becomes 2 (by writing some commands in code itself). Can I extract such values of t as a vector?

clear
clc;
N=2000;  % Number of time stpes
dt=0.05;  %Time step
%Initialize
 y=zeros(N,1); 
 t=zeros(N,1);
%Initial values
t(1)=0;
y(1)=1;

      for j=2:N
            t(j)=t(j-1)+dt;
      y(j)=y(j-1)+dt*sin(t(j-1));
      end
      plot(t,y)

Answer 1

[Star Strider]

Try this:

N=2000;  % Number of time stpes
dt=0.05;  %Time step
%Initialize
y=zeros(N,1);
t=zeros(N,1);
%Initial values
t(1)=0;
y(1)=1;

for j=2:N
    t(j)=t(j-1)+dt;
    y(j)=y(j-1)+dt*sin(t(j-1));
end

zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                                    % Returns Approximate Zero-Crossing Indices Of Argument Vector
idx_2 = zci(y - 2);                                                                     % Vector Of Approximate Indices Where ‘y=2’
for k = 1:numel(idx_2)
    t2(k) = interp1(y(idx_2(k)+[-1 1]), t(idx_2(k)+[-1 1]), 2, 'linear','extrap');      % Interopolate: ‘t’ Where ‘y=2’
end

figure
plot(t,y)
hold on
plot(t2, 2*ones(size(t2)), 'pg')
hold off

This uses the ‘zci’ utility function to find the approximate indices where ‘t=2’, then uses the ‘idx_2’ index vector to loop through your data to find the closest values of ‘t’ that the interp1 function can calculate where ‘y=2’.

[Dharma]

Thank you so much Star. In my actual problem y goes from 0 to 1. How do I modify above equation for finding t2 values at which y=0.4?

[Star Strider]

My pleasure.

Change the ‘idx_2’ calculation to:

idx_2 = zci(y - 0.4);                                                                     % Vector Of Approximate Indices Where ‘y=0.4’

[Dharma]

It is working now. Thank you again Star!

[Star Strider]

My pleasure.

If my Answer helped you solve your problem, please Accept it!

[Dharma]

 I am requesting you for one more favor.
 C=1;
  for j=2:N
      t(j)=t(j-1)+dt;
      y(j)=y(j-1)+dt*sin(t(j-1))+C;
  end
  zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);  
  idx_2 = zci(y - 2);                                                                   
  for k = 1:numel(idx_2)
      t2(k) = interp1(y(idx_2(k)+[-1 1]), t(idx_2(k)+[-1 1]), 2, 'linear','extrap');   
  end
  % I want to use C=2 at t=t2 otherwise C=1 (continue with upper loop of j).
  % Which probably means I can terminate upper j loop at the end of code only. 
% If I repeat j loop like below then there is no calculation for C=1 
% for t not equal to t2. Please suggest me.
  for j=2:N
      t(j)=t(j-1)+dt;
      for k=1:numel(idx_2)
          if t(j)==t2(k)
              C=2;
            y(j)=y(j-1)+dt*sin(t(j-1))+C;
          end
      end
  end

[Star Strider]

I have no idea what you want to do. I cannot run the code you posted in your Comment, even using your previous value for the constants.

As best as I can understand, you want to add 1 to every ‘y’ value where ‘y=2’, and 2 to every ‘y’ value where ‘y~=2’. I would just do that in the vector of results.

Note that ‘t2’ may not be an element of ‘t’, since it is an interpolated value.

[Dharma]

% Sorry for confusion. I am trying to use t2 values for applying
 % if condition on C. If t=t2 then C=2 otherwise C=1. This is my
 % running code. Main problem here is for values of t other than 
% t2, y(j) is calculated using second (last) j loop not using
% first (uppermost) j loop. Any hints?

N=2000;  % Number of time stpes
dt=0.05;  %Time step
%Initialize
 y=zeros(N,1); 
 t=zeros(N,1);
%Initial values
t(1)=0;
y(1)=1;
 C=1;
    for j=2:N
          t(j)=t(j-1)+dt;
    y(j)=y(j-1)+dt*C*sin(t(j-1));
    end
    zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                                  
    % Returns Approximate Zero-Crossing Indices Of Argument Vector
    idx_2 = zci(y - 2);                                                                    
% Vector Of Approximate Indices Where ‘y=2’
for k = 1:numel(idx_2)
    t2(k) = interp1(y(idx_2(k)+[-1 1]), t(idx_2(k)+[-1 1]), 2, 'linear','extrap');      
    % Interopolate: ‘t’ Where ‘y=2’
end
 for j=2:N
     t(j)=t(j-1)+dt;
    for k=1:numel(idx_2)
        if t(j)==t2(k)
            C=2;
          y(j)=y(j-1)+dt*C*sin(t(j-1));
        end
    end
 end
plot(t,y)

[Star Strider]

Since ‘t2’ is an interpolated value, ‘t’ will never equal ‘t2’. You can determine this with:

Lia = ismember(t, t2);                                  % Logocal Indices Where ‘t==t2’
Lnz = nnz(Lia);                                         % Number Of Times That Occurs

The ‘Lnz’ result is 0, indicating that there are no values of ‘t’ and ‘t2’ that are the same.

Introducing a tolerance:

Liat = ismembertol(t, t2, 0.00025);
Lnzt = nnz(Liat);

reveals that there are 32 of them are close (within 0.00025).

[Dharma]

Thanks again. Now I have realized my error here.

   Final question: How do I write 
        if t(j)==t2(k)+ some tolerance (say 0.5)
      I am trying to select all value of t  which are within 0.5 (or another number) of t2 value.

[Star Strider]

My pleasure.

I would use the ismembertol function, as I illustrated, with the tolerance I used.

[Dharma]

Hi Star,
You have already answered my question correctly.
I just have one curiosity.
Regarding your answer (following two lines), what I am thinking is 
Liat vector will return 1 for those value of t lying around t2 value by 
±0.00025. For example I am setting tolerance equal to 0.02 
and first element of t2 is 10.0366. Now my expectation is that for the 
value of t within ±0.02 from t2=10.0366, logical matrix Liat should return 1. 
But in my case it has returned value 1 for t values starting from around 8.
Please give me suggestion what basic concept I am missing here.
Thanks

    Liat = ismembertol(t, t2, 0.00025);
    Lnzt = nnz(Liat);

[Star Strider]

The ‘tolerance’ value has the effect of returning elements

t = t2 ± tolerance

Setting a very small tolerance (for example 0.00025) will result in a very small number of returned values. Setting a much larger tolerance (for example 0.02) will result in a much larger number of returned values.

[Dharma]

Hi Star,

I am requesting you for your suggestion on the following code which uses your previous suggestion for extrapolation. My question is why does code give error message (in the line of t2 calculation i.e. inside for loop of m) when I change number of oscillators (N) from 50 to 60? Lets just ignore other parameters and loop. Thank you!

Dharma

N=60;   % number of oscillators
dt=0.05;  %time step
J=2000;  %number of time steps
%Initialize matrices
t=zeros(1,J);
theta=zeros(N,J);
r=zeros(1,J);
sigN=0.2; %standard deviation of natural frequencies distribution
rng default
 theta(:,1)=2*pi*rand(N,1);%random initial values of theta
 W = random('Normal',0,sigN,N,1);% random values drawn from Gaussian
 K=zeros(1,J);
 r(1)=1;
    K(1)=0.3;         
    for j=2:J
        t(j)=t(j-1)+dt;
        K(j)=0.3;
        kura = @(x,N) W+(K(j)/N)*...
        sum(sin(meshgrid(x)-meshgrid(x)'),2);
        %4th order Runge-Kutta method
        k1=kura(theta(:,j-1),N);
        k2=kura(theta(:,j-1)+0.5*dt*k1,N);
        k3=kura(theta(:,j-1)+0.5*dt*k2,N);
        k4=kura(theta(:,j-1)+dt*k3,N);

          theta(:, j)=theta(:,j-1)+(dt/6)*(k1+2*k2+2*k3+k4);
          r(j)=1/N*sqrt((sum(cos(theta(:,j))))^2+(sum(sin(theta(:,j))))^2);
      end
          zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);  
          idx_2 = zci(r - 0.2);                                                                   

      for m = 1:numel(idx_2)
          t2(m) = interp1(r(idx_2(m)+[-1 1]), t(idx_2(m)+[-1 1]), 0.2, 'linear','extrap'); 
       % Interpolate t where r=0.4
      end

[Star Strider]

You included the [-1 1] offsets as part of the subscripts for ‘r’ and ‘t’.

Try this instead:

t2(m) = interp1(r(idx_2(m))+[-1 1], t(idx_2(m))+[-1 1], 0.2, 'linear','extrap');

[Dharma]

Thanks. Actually I was using the same as you suggested (Except in comment where I wrote r=0.4). I found solution as well. If I choose value of r smaller than 0.2 then I will be able to increase N. Do you have any suggestion for this?

   zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);  
            idx_2 = zci(r - 0.2);                                                                   
        for m = 1:numel(idx_2)
            t2(m) = interp1(r(idx_2(m)+[-1 1]), t(idx_2(m)+[-1 1]), 0.2, 'linear','extrap'); 
         % Interpolate t where r=0.4
        end
Also in the line interp1(r(idx_2(m))+[-1 1], what does [-1 1] represent?
Thank you again!

[Star Strider]

My pleasure.

‘If I choose value of r smaller than 0.2 then I will be able to increase N. Do you have any suggestion for this?’

I have no idea what you are doing, so I cannot answer that.

‘Also in the line interp1(r(idx_2(m))+[-1 1], what does [-1 1] represent?’

That vector creates a region ±1 indices around the chosen point, creating a region over which interp1 can interpolate. It is designed to create a strictly monotonic region in a function that may not otherwise be monotonic.

Answer 2

[Benjamin Großmann]

You could rewrite your function so that your wanted data points are maxima and search them with findpeaks.

 [~,locs] = findpeaks(-abs(y-2));
 plot(t,y,t(locs),y(locs),'or')

The pretty simple and elegant t2 = t(y==2) does not work here, since your data is not exaxtly two at any given entry.