Limit cycles

Run using Octave GUI v6.1.0

x = -1; y = -1;
A = [[0, 1-3*y**2]; [-1, -2*y]];
[v, lam] = eig(A);
disp(lam)
disp(v)

Diagonal Matrix

  -0.73205         0
         0   2.73205
  -0.93907   0.59069
  -0.34372  -0.80690

• We now consider the Manta Ray System
• The system of equations is given by: $$ \dot{x} = y - y^3$$ $$ \dot{y} = -x-y^2$$
• The fixed points for this system are: (0,0), (-1,-1) and (-1,1)
• We now linearize the system near the fixed points to understand the behaviour. We generate the Jacobian near each fixed point
• This example highlights the aspect of time-reversal symmetry

x = linspace(-3, 3, 100); y = linspace(-3, 3, 100);
[X, Y] = meshgrid(x, y);
U = Y-Y.^3;
V = -X-Y.^2;
#seedpoints
x1 = linspace(min(x), max(x), 20);
y1 = linspace(min(y), max(y), 20);
[X1, Y1] = meshgrid(x1,y1);

streamline(X, Y, U, V, [X1], [Y1]);
daspect([1 1 1]);

• Let us now look at the case of a simple pendulum
• The governing equation of motion is: $$ \ddot{\theta} + \sin{\theta} = 0 $$
• Let $x =\theta$, $y = \dot{x}$. We now have the following system of equations: $$ \dot{x} = y $$ $$ \dot{y} = -\sin{x}$$
• The fixed points of this system are $y = 0$ and $x = n\pi$. The system is infinitely periodic in the x direction so we can analyze only one time period and understand the behaviour of the system.
• The center (0,0) is a robust center since the equations have time-reversal symmetry.
• A trajectory which connects two fixed points is known as heteroclinic trajectory. In this example the saddle points at $y=0$ and $x = (2n+1)\pi$ are connected by such trajectories.

x = linspace(-6, 6, 100); y = linspace(-3, 3, 100);
[X, Y] = meshgrid(x, y);
U = Y;
V = -sin(X);
x1 = linspace(min(x), max(x), 20);
y1 = linspace(min(y), max(y), 20);
[X1, Y1] = meshgrid(x1,y1);
streamline(X, Y, U, V, [X1], [Y1]);
daspect([1 1 1]);

• Let us now add a damper to the previous system
• Now the system is: $$ \dot{x} = y $$ $$ \dot{y} = -\sin{x} - by$$
• The origin now is no longer a closed orbit but an inward spiral.
• There are some trajectories which can skip past the origin but they eventually spiral in at some other fixed point with $y=0$ and $x=2n\pi$
• As we increase the damping by increasing b, we notice that trajectories at larger energies are now spiraling in at the origin.

clear; clc; close all;

h.ax = axes ("position", [0.05 0.4 0.5 0.5]);  #reduce plot windows size to accomodate sliders

function update_plot (obj, init = false)
  ## gcbo holds the handle of the control
  h = guidata (obj);
  replot = false;
  recalc = false;
  switch (gcbo)                       # If we make any change then we replot
    case {h.print_pushbutton}
      fn =  uiputfile ("*.png");
      print (fn);
    case {h.slider1}
      recalc = true;
    end
    
    if (recalc || init)  
      b = 0 + 1*get (h.slider1, "value");
      set (h.label1, "string", sprintf ("b: %1f", b));
      x = linspace(-12, 12, 30); y = linspace(-3, 3, 30);
      [X, Y] = meshgrid(x, y);
      U = Y;
      V = -sin(X)-b*Y;
      x1 = linspace(min(x), max(x), 20);
      y1 = linspace(min(y), max(y), 20);
      [X1, Y1] = meshgrid(x1,y1);
      streamline(X, Y, U, V, [X1], [Y1]);
      daspect([1 1 1]);
      
      
      
    else
      set (h.plot, "ydata", y);
    end
  end
  
  # specify the GUI plot properties like sliders and label formatting
  
  ## print figure
  h.print_pushbutton = uicontrol ("style", "pushbutton",
  "units", "normalized",
  "string", "print plot\n(pushbutton)",
  "callback", @update_plot,
  "position", [0.6 0.45 0.35 0.09]);
  ## guess
  h.label1 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Guess:",
  "horizontalalignment", "left",
  "position", [0.05 0.25 0.35 0.08]);
  
  h.slider1 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.1,
  "position", [0.05 0.20 0.35 0.06]);
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);

• We can plot the time series

# Set default parameters for the plots
set(0, "defaultlinelinewidth", 5);
set (0, "defaulttextfontname", "TimesNewRoman")
set (0, "defaulttextfontsize", 20)
set (0, "DefaultAxesFontName", "TimesNewRoman")
set(0, 'DefaultAxesFontSize', 20)

function y=mysys(x,t) # returns the RHS
    y=([x(2), -0.25*x(2)-sin(x(1))]);
end

tspan = [0,10];
x0 = 0; 
y0 = 5;
ics = [x0, y0];
tspan=linspace(0,10,100)';
f = @(x,t) mysys(x, t);
sol = lsode(f, ics, tspan);

tout = linspace(0, max(tspan), 100);
xout = sol(:,1);
yout = sol(:,2);
E = 1/2*yout.^2 - cos(xout);
plot(tout, xout, tout, yout);

We now focus our discussion to Limit Cycles

• Consider the following set of equations: $$ \dot{r} = r(1-r^2)$$ $$ \dot{\theta} = 1$$
• The fixed point is $r=1$. The $1-r^2$ term prevents the trajectory from blowing up
• We solve the equation in polar coordinates and then cast the solution to cartesian coordinates using the transformation $x = r\cos(\theta)$ and $y = r\sin(\theta)$
• We take the initial condition away from the unit circle.
• The unit circle in this case is a limit cycle as all the trajectory ultimately converege onto the circle. A limit cycle is an isolated closed trajectory

function y=mysys(r, t) # returns the RHS
    y= ([r(1)*(1-r(1).^2), 1]);
end

tspan=linspace(0,10,100)';
x0 = 0; 
y0 = 2;
r0 = (x0.^2 + y0.^2).^0.5; theta0 = atan2(y0, x0);

ics = [r0, theta0];
f = @(r,t) mysys(r, t);
sol = lsode(f, ics, tspan);

tout = linspace(0, max(tspan), 100);
rout = sol(:,1);
thetaout = sol(:,2);

xout = rout.*cos(thetaout); yout = rout.*sin(thetaout);

subplot (2, 1, 1);
plot(xout, yout);
daspect([1 1 1]);
subplot(2,1,2);
plot(tout, xout, tout, yout);

• We can plot a group of trajectories by varying the initial point and see that all trajectories indeed home in on the unit circle
• Limit cycles need not be always circular in nature.
• Let us now look at an example of Van der Pol oscillator whose governing equation is: $$\ddot{x} + \mu \dot{x}(x^2-1) + x =0$$
  • It has a nonlinear damping.
• Limit cycles are present due to feedback terms. A linear equation cannot have limit cycles.

clear; clc; close all;

h.ax = axes ("position", [0.05 0.4 0.5 0.5]);  #reduce plot windows size to accomodate sliders
function y=mysys(x, t, a) # returns the RHS
  y= ([x(2), -x(1)-a*(x(1).^2-1)*x(2)]);
end


function update_plot (obj, init = false)
  ## gcbo holds the handle of the control
  h = guidata (obj);
  replot = false;
  recalc = false;
  switch (gcbo)                       # If we make any change then we replot
    case {h.print_pushbutton}
      fn =  uiputfile ("*.png");
      print (fn);
    case {h.slider1}
      recalc = true;
    case {h.slider2}
      recalc = true;
     case {h.slider3}
      recalc = true;
    end
    
    if (recalc || init)  
      x0 = -1.5+ 3*get (h.slider1, "value");
      y0 = -1.5 + 3*get (h.slider2, "value");
      a =  -1.5 + 3*get (h.slider3, "value");
      set (h.label1, "string", sprintf ("x0: %.1f", x0));
      set (h.label2, "string", sprintf ("y0: %.1f", y0));
      set (h.label3, "string", sprintf ("a: %.1f", a)); 
      
      x = linspace(-4, 4, 14); y = linspace(-4, 4, 14);
      [X, Y] = meshgrid(x, y);
      U = Y;
      V = -X-a*(X.^2-1).*Y;
      quiver(X, Y, U./(U.^2+V.^2).^0.5, V./(U.^2+V.^2).^0.5);
      daspect([1 1 1]);
      hold on;
      
      tspan=linspace(0,20,100)';    
      ics = [x0, y0];
      f = @(x,t) mysys(x, t, a);
      sol = lsode(f, ics, tspan);
      
      tout = linspace(0, max(tspan), 100);
      xout = sol(:,1);
      yout = sol(:,2);
      
      plot(xout,yout)
      hold off;
      
      
      
    else
      set (h.plot, "ydata", y);
    end
  end
  
  # specify the GUI plot properties like sliders and label formatting
  
  ## print figure
  h.print_pushbutton = uicontrol ("style", "pushbutton",
  "units", "normalized",
  "string", "print plot\n(pushbutton)",
  "callback", @update_plot,
  "position", [0.6 0.45 0.35 0.09]);
  ## guess
  h.label1 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Guess:",
  "horizontalalignment", "left",
  "position", [0.05 0.25 0.35 0.08]);
  
  h.slider1 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.1,
  "position", [0.05 0.20 0.35 0.06]);
  
   h.label2 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Iteration:",
  "horizontalalignment", "left",
  "position", [0.05 0.15 0.35 0.08]);
  
  h.slider2 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.5,
  "position", [0.05 0.10 0.35 0.06]);
  
    h.label3 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Guess:",
  "horizontalalignment", "left",
  "position", [0.6 0.25 0.35 0.08]);
  
  h.slider3 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.1,
  "position", [0.6 0.20 0.35 0.06]);
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);

• Let us also plot the trajectories

clear; clc; close all;

h.ax = axes ("position", [0.05 0.4 0.5 0.5]);  #reduce plot windows size to accomodate sliders
function y=mysys(x, t, a) # returns the RHS
  y= ([x(2), -x(1)-a*(x(1).^2-1)*x(2)]);
end


function update_plot (obj, init = false)
  ## gcbo holds the handle of the control
  h = guidata (obj);
  replot = false;
  recalc = false;
  switch (gcbo)                       # If we make any change then we replot
    case {h.print_pushbutton}
      fn =  uiputfile ("*.png");
      print (fn);
    case {h.slider1}
      recalc = true;
    case {h.slider2}
      recalc = true;
     case {h.slider3}
      recalc = true;
    end
    
    if (recalc || init)  
      x0 = -1.5+ 3*get (h.slider1, "value");
      y0 = -1.5 + 3*get (h.slider2, "value");
      a =  -1.5 + 3*get (h.slider3, "value");
      set (h.label1, "string", sprintf ("x0: %.1f", x0));
      set (h.label2, "string", sprintf ("y0: %.1f", y0));
      set (h.label3, "string", sprintf ("a: %.1f", a)); 
      
      x = linspace(-4, 4, 14); y = linspace(-4, 4, 14);
      [X, Y] = meshgrid(x, y);
      U = Y;
      V = -X-a*(X.^2-1).*Y;
      
      tspan=linspace(0,20,100)';    
      ics = [x0, y0];
      f = @(x,t) mysys(x, t, a);
      sol = lsode(f, ics, tspan);
      
      tout = linspace(0, max(tspan), 100);
      xout = sol(:,1);
      yout = sol(:,2);
      
      plot(tout, xout, tout, yout)
      
      
      
    else
      set (h.plot, "ydata", y);
    end
  end
  
  # specify the GUI plot properties like sliders and label formatting
  
  ## print figure
  h.print_pushbutton = uicontrol ("style", "pushbutton",
  "units", "normalized",
  "string", "print plot\n(pushbutton)",
  "callback", @update_plot,
  "position", [0.6 0.45 0.35 0.09]);
  ## guess
  h.label1 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Guess:",
  "horizontalalignment", "left",
  "position", [0.05 0.25 0.35 0.08]);
  
  h.slider1 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.1,
  "position", [0.05 0.20 0.35 0.06]);
  
   h.label2 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Iteration:",
  "horizontalalignment", "left",
  "position", [0.05 0.15 0.35 0.08]);
  
  h.slider2 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.5,
  "position", [0.05 0.10 0.35 0.06]);
  
    h.label3 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Guess:",
  "horizontalalignment", "left",
  "position", [0.6 0.25 0.35 0.08]);
  
  h.slider3 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.1,
  "position", [0.6 0.20 0.35 0.06]);
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);

• The trajectories has multiple time scales

• Let us now study the Glycolysis reaction
• It has been found that it can be modeled as 2 coupled nonlinear equations: $$ \dot{x} = -x + ay + x^2y$$ $$ \dot{y} = b -ay -x^2y$$ where, $x$ represents evolution of ATP and $y$ represents evolution of ADP
• The fixed point for this case by be found by first adding both and equations and putting $\dot{x}$ and $\dot{y}$ equal to 0. We obtain $x^* = b$ and $y^* = b/(a+b^2)$
• We plot the time series and the phase plot. Red circle denotes the fixed point and black circle denotes the initial point

• NOTE: using subplot with uicontrol is tricky. When subplot overlaps with another uicontrol object it deletes the axis configuration that we had specified. So that is why we have used a 4x4 subplot instead of 2x2 subplot for our two plots

clear; clc; close all;

h.ax = axes ("position", [0.05 0.4 0.5 0.5]);  #reduce plot windows size to accomodate sliders
function y=mysys(x, t, a, b) # returns the RHS
  y= ([-x(1) + a*x(2) + x(1).^2.*x(2), b - a*x(2) - x(1).^2.*x(2)]);
end


function update_plot (obj, init = false)
  ## gcbo holds the handle of the control
  h = guidata (obj);
  replot = false;
  recalc = false;
  switch (gcbo)                       # If we make any change then we replot
    case {h.slider1}
      recalc = true;
    case {h.slider2}
      recalc = true;
     case {h.slider3}
      recalc = true;
      case {h.slider4}
      recalc = true;
    end
    
    if (recalc || init)  
      x0 = 0+ 1.5*get (h.slider1, "value");
      y0 = 0 + 1.5*get (h.slider2, "value");
      a =  0 + 0.2*get (h.slider3, "value");
      b =  0 + 1.5*get (h.slider4, "value");
      set (h.label1, "string", sprintf ("x0: %.1f", x0));
      set (h.label2, "string", sprintf ("y0: %.1f", y0));
      set (h.label3, "string", sprintf ("a: %1f", a)); 
      set (h.label4, "string", sprintf ("b: %1f", b)); 

      x = linspace(0, 4, 14); y = linspace(0, 4, 14);
    [X, Y] = meshgrid(x, y);
    U = -X + a*Y + X.^2.*Y;
    V = b - a*Y - X.^2.*Y;
      
      tspan=linspace(0,40,100)';    
      ics = [x0, y0];
      f = @(x,t) mysys(x, t, a, b);
      sol = lsode(f, ics, tspan);
      
      tout = linspace(0, max(tspan), 100);
      xout = sol(:,1);
      yout = sol(:,2);
      
      subplot(2,2,1);
      quiver(X, Y, U./(U.^2+V.^2).^0.5, V./(U.^2+V.^2).^0.5);
      xlim([0 2.5]);
      daspect([1 1 1]);
      hold on; 
      plot(xout, yout, x0, y0, 'ok',b, b/(a+b^2), 'sr');
      hold off;
            
      subplot(2,2,2);
      plot(tout, xout, tout, yout);
      
    else
      set (h.plot, "ydata", y);
    end
  end
  
  # specify the GUI plot properties like sliders and label formatting
  ## guess
  h.label1 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Guess:",
  "horizontalalignment", "left",
  "position", [0.05 0.25 0.35 0.08]);
  
  h.slider1 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.1,
  "position", [0.05 0.20 0.35 0.06]);
  
   h.label2 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Iteration:",
  "horizontalalignment", "left",
  "position", [0.05 0.15 0.35 0.08]);
  
  h.slider2 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.5,
  "position", [0.05 0.10 0.35 0.06]);
  
    h.label3 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Guess:",
  "horizontalalignment", "left",
  "position", [0.6 0.25 0.35 0.08]);
  
  h.slider3 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.1,
  "position", [0.6 0.20 0.35 0.06]);
  
    h.label4 = uicontrol ("style", "text",
  "units", "normalized",
  "string", "Iteration:",
  "horizontalalignment", "left",
  "position", [0.6 0.15 0.35 0.08]);
  
  h.slider4 = uicontrol ("style", "slider",
  "units", "normalized",
  "string", "slider",
  "callback", @update_plot,
  "value", 0.5,
  "position", [0.6 0.10 0.35 0.06]);
  
  
  set (gcf, "color", get(0, "defaultuicontrolbackgroundcolor"))
  guidata (gcf, h)
  update_plot (gcf, true);