Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Values of intersection points of plot. Print results.

Asked by Pauline on 22 Apr 2013

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

Cedric Wannaz on 22 Apr 2013

I am not sure to fully understand your question.. you get a vector x and you want to print elements x(i-1) for all i such as x(i) equals 1?

Pauline on 23 Apr 2013

Hi, Cedric. Not quite, perhaps my edit to the question would help clarify my problem?

Walter Roberson on 23 Apr 2013

Saying x(:,1) and x(:,2) would make the question clearer.

Pauline

Tags

Products

No products are associated with this question.

2 Answers

Answer by Teja Muppirala on 23 Apr 2013
Accepted answer

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;

1 Comment

Pauline on 23 Apr 2013

Thank you, Teja!

Teja Muppirala
Answer by Cedric Wannaz on 23 Apr 2013
Edited by Cedric Wannaz 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 Wannaz 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.

Pauline on 23 Apr 2013

Thanks, Cedric, Walter.

Cedric Wannaz

Contact us