How can I plot a phase portrait of x' vs x

6 views (last 30 days)
Liam Wheen
Liam Wheen on 23 Mar 2017
Answered: SK on 24 Mar 2017
I am conducting a project on friction looking at the stick/slip phenomena in a similar system to that linked below.
I have written the code in order to plot the state space (x against x') given some initial conditions for position and velocity. I would now like to plot a phase portrait to show the general direction of motion in a similar fashion to the one linked below.
However all of the help I have found online does not explain how to do this when one axis is just the derivative of the other.
The general pattern of behaviour is two ellipses that have centres at +-d/w^2 on the x axis where d=mu*g and w^2 = k/m where k is the spring stiffnes, mu is the kinetic coefficient of friction, g is gravity, and m is mass. The RHS ellipse only exists below the line x'=v where v is the constant belt velocity and the LHS ellipse only exists above this line.
I would like to have arrows showing what will happen to the values throughout the state space for some decided parameters.
Below is the code I made to calculate the x and x' values, if you run it with
X = Coulomb3([-2 4], 6, 0.005, 10, 2, 1, 0.3);
then you should get an idea of where the ellipses are placed.
function [X] = Coulomb3(Init, Tf, tstep, k, m, v, mu)
dmu = mu;
smu = 1.3*mu;
g = 9.81;
w = sqrt(k/m);
c = smu*m*g/k;
d = dmu*g;
tol = 0.05;
X(1, 1) = Init(1);
X(1, 2) = Init(2);
while length(X) <= Tf/tstep
sc = 1;
if X(end, 2) < v - tol || (abs(X(end, 2) - v) < tol && X(end, 1) >= c)
A = X(end, 2)/w;
if A<0; sc=-1; end
B = X(end, 1) - d/w^2;
R = sqrt(A^2 + B^2);
phi = atan(B/A);
for t = 1:Tf/tstep
if sc*cos(w*t*tstep + phi) > v/(R*w)
T = t; %finding where the line first touches the x'=v line
break;
end
T = Tf/tstep;
end
for t = 1:T
X(end+1, :) = [sc*R*sin(w*t*tstep + phi) + d/w^2 sc*R*w*cos(w*t*tstep + phi)];
end
end
if abs(X(end, 2) - v) < tol && abs(X(end, 1)) < c
for t = 1:Tf/tstep
if X(end, 1) + t*tstep*v > c
T = t; %finding where the x'=v line runs out (at c)
break;
end
T = Tf/tstep;
end
for t = 1:T
X(end+1, :) = [(X(end, 1) + tstep*v) v];
end
end
if X(end, 2) > v + tol || (abs(X(end, 2) - v) < tol && X(end, 1) <= -c)
A = X(end, 2)/w;
B = X(end, 1) + d/w^2;
R = sqrt(A^2 + B^2);
phi = atan(B/A);
for t = 1:Tf/tstep
if cos(w*t*tstep + phi) < v/(R*w)
T = t; %finding where the line first touches the x'=v line
break;
end
T = Tf/tstep;
end
for t = 1:T
X(end+1, :) = [R*sin(w*t*tstep + phi) - d/w^2 R*w*cos(w*t*tstep + phi)];
end
end
end
end
Any help is greatly appreciated as I would like to have this last figure to conclude my project.
Thanks

Answers (1)

SK
SK on 24 Mar 2017
It is a little tricky to get good looking arrows. As a first iteration you could do this - an example with x = t and x' = t^2, but it could be any other function like the one you calculate:
t = transpose(1:0.1:5);
x = t;
xprime = t.^2;
u = diff(x);
v = diff(xprime);
u = [u; u(end)];
v = [v; v(end)];
scale = 0.5;
quiver(x,xprime,u,v, scale);
This gives arrows whose size is proportional to the size of the derivative. Also the arrows are too closely spaced. To make them less closely spaced I did the following:
close all
w = sqrt(u.^2 + v.^2);
cw = cumsum(w);
num_arrows = 5;
for i = 1 : num_arrows
locs(i) = find(cw >= cw(end)*i/num_arrows, 1);
end
plot(x, xprime);
hold on
quiver(x(locs),xprime(locs),u(locs),v(locs), scale);
The only issue now is that the arrows vary in size depending on the size of the derivative. This probably can also be fixed - but you may have to put in a little more effort.

Categories

Find more on Programming in Help Center and File Exchange

Products

Community Treasure Hunt

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

Start Hunting!