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?

6 views (last 30 days)
Dharma
Dharma on 23 Apr 2018
Commented: Star Strider on 30 May 2018
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)

Accepted Answer

Star Strider
Star Strider on 23 Apr 2018
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’.
  16 Comments
Star Strider
Star Strider on 30 May 2018
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.

Sign in to comment.

More Answers (1)

Benjamin Großmann
Benjamin Großmann on 23 Apr 2018
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.

Community Treasure Hunt

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

Start Hunting!

Translated by