Values of intersection points of plot. Print results.
2 views (last 30 days)
Show older comments
Suppose I have a code that numerically integrates a set of ODEs, where my dependent variable is t and the vector of variables is x. I wish to print and plot values of x(1) for whenever x(2)=1. How might I achieve that? Some sort of "If" function?
In other words, graphically speaking, I can plot x(1) against x(2) and draw a vertical line at x(2)=1 and I want to know the x(1) values at the intersections of the plot and the vertical line.
Thanks.
3 Comments
Accepted Answer
Teja Muppirala
on 23 Apr 2013
If your x is the output of an ODE solver, then it might not hit x2 = 1 exactly and you will need to interpolate.
You could use a line intersection finder from the file exchange or the Mapping Toolbox's POLYXPOLY function, or you could kind of do it manually:
% Just making some test data
A = [-0.1 -2; 2 -0.2];
F = @(t,x) A*x;
[T,X] = ode45(F,[0 10],[5; 2]);
plot(X(:,2),X(:,1)); xlabel('X(2)'),ylabel('X(1)');
hold on;
x2 = X(:,2); % Get the second state
x2goal = 1; % Position of vertical line
% Find zero crossings by changes in sign
E = x2-x2goal; % Find when E == 0
dE = diff(sign(E));
t = find(abs(dE)==2);
t = [t + E(t)./(E(t)-E(t+1)); find(E==0)]; %Linear Interpolation
t = sort(t); %Not necessary
crossings = interp1(X(:,1),t)
% Plot the result
plot(x2goal,crossings,'r.','markersize',16)
grid on;
More Answers (1)
Cedric
on 23 Apr 2013
Edited: Cedric
on 23 Apr 2013
EDIT: this answer is wrong, I answered too quickly without paying attention to the fact that your x is the output of an ODE solver. Please jump directly to Teja's answer.
Ok, then you can proceed as follows:
id = x(:,2) == 1 ;
x(id,1) % This will provide you with all x(:,1) that
% correspond to a x(:,2) equal to 1.
id is a vector of logicals (outcome of the test of equality on x(:,2)), that we use for indexing x(:,1). If you want to have the positions numerically, you can use
find(id)
3 Comments
Cedric
on 23 Apr 2013
Edited: Cedric
on 23 Apr 2013
Ah yes, if x gets out of an ODE solver, Walter's comment and Teja's answer are really important; I should have thought a little further.
EDIT: and my answer is even wrong if you don't interpolate, because it is absolutely not trivial to quantify how far from 1 points that are from each side of 1 can be. Setting a large tolerance would be wrong too as it could lead to taking points that are not around any 1 crossing.
See Also
Categories
Find more on Historical Contests in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!