Search code examples
matlabmathplotcycle

Limit cycle in Matlab for 2nd order autonomous system


I wish to create a limit cycle in Matlab. A limit cycle looks something like this:

enter image description here

I have no idea how to do it though, I've never done anything like this in Matlab.

The equations to describe the limit cycle are the following:

x_1d=x_2
x_2d=-x_1+x_2-2*(x_1+2*x_2)x_2^2

It is to be centered around the equilibrium which is (0,0)

Can any of you help me?


Solution

  • If you use the partial derivatives of your function to make a vector field, you can then use streamlines to visualize the cycle that you are describing.

    For example, the function f = x^2+y^2

    Gives me partial derivatives dx = 2x, dy=2y For the visualization, I sample from the partial derivatives over a grid.

    [x,y] = meshgrid(0:0.1:1,0:0.1:1);
    dx = 2*x;
    dy = 2*y;
    

    To visualize the vector field, I use quiver;

    figure; 
    quiver(x, y, dx, dy); 
    

    Using streamline, I can visualize the path a particle injected into the vector field would take. In my example, I inject the particle at (0.1, 0.1)

    streamline(x,y, dx, dy, 0.1, 0.1);
    

    This produces the following visualization

    Vector field streamline visualization example

    In your case, you can omit the quiver step to remove the hedgehog arrows at every grid point.

    Here's another example that shows the particle converging to an orbit.

    Vector field visualization example

    Edit: Your function specifically. So as knedlsepp points out, the function you are interested in is a bit ambiguously stated. In Matlab, * represents the matrix product while .* represents the element-wise multiplication between matrices. Similarly, '^2' represents MM for a matrix M, while .^2 represents taking the element-wise power.

    So,

    [x_1,x_2] = meshgrid(-4:0.1:4,-4:0.1:4);
    dx_1 = x_2;
    dx_2 = -x_1+x_2-2*(x_1+2*x_2)*(x_2)^2;
    figure; streamline(x_1,x_2, dx_1, dx_2, 0:0.1:4, 0:0.1:4);
    

    Looks like

    Vector field visualization

    This function will not show convergence because it doesn't converge.

    knedlsepp suggests that the function you are actually interested in is

    dx_1 = -1 * x_2;
    dx_2 = -1 * -x_1+x_2-2*(x_1+2*x_2).*(x_2).^2;    
    

    His post has a nice description of the rest.