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

# 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.

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.

## Products

No products are associated with this question.

Answer by 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;

## 1 Comment

Pauline on 23 Apr 2013

Thank you, Teja!

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)

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.