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

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

Answer 1

[SK]

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.