Limit cycles

import numpy as np;
import matplotlib.pyplot as plt;
plt.rcParams.update({"text.usetex":True});
%config InlineBackend.figure_format = "svg"
from ipywidgets import interactive
from scipy.integrate import solve_ivp
from mayavi import mlab

x = -1; y = -1;
A = np.array([[0, 1-3*y**2],[-1, -2*y]])
l, v = np.linalg.eig(A)
print(l)
print(v)

[-0.73205081  2.73205081]
[[-0.9390708   0.59069049]
 [-0.34372377 -0.80689822]]

• 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 = np.linspace(-3, 3, 100); y = np.linspace(-3, 3, 100);
X, Y = np.meshgrid(x, y);
U = Y-Y**3;
V = -X-Y**2;
#plt.quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
plt.streamplot(X, Y, U, V, integration_direction='forward', density=1.5)
ax = plt.gca(); ax.set_aspect(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 = np.linspace(-6, 6, 100); y = np.linspace(-3, 3, 100);
X, Y = np.meshgrid(x, y);
U = Y;
V = -np.sin(X);
#plt.quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
plt.streamplot(X, Y, U, V, integration_direction='forward', density=1.5)
ax = plt.gca(); ax.set_aspect(1);

• We vary $y$ and plot the velocity along the z xis

def mysys(t, x): # returns the RHS
    return ([x[1], -np.sin(x[0])]);

tspan = [0,40]
x0 = 0; 
for y0 in np.linspace(-2.5, 2.5,100):
    ics = [x0, y0]
    sol = solve_ivp(mysys, tspan, ics, dense_output=True, rtol=1e-8);

    tout = np.linspace(0, np.max(tspan), 200);
    xout = sol.sol(tout)[0];
    yout = sol.sol(tout)[1];
    E = 1/2*yout**2 - np.cos(xout);
    mlab.plot3d(np.cos(xout), np.sin(xout), yout )
    
mlab.show();

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

def show_2dphase(b = 0.1):
    x = np.linspace(-12, 12, 30); y = np.linspace(-3, 3, 30);
    X, Y = np.meshgrid(x, y);
    U = Y;
    V = -b*Y-np.sin(X);
    #plt.quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
    plt.streamplot(X, Y, U, V, integration_direction='forward', density=1.5)
    ax = plt.gca(); ax.set_aspect(1);

w = interactive(show_2dphase, b=(0, 1, 0.05))
w

• We can plot the time series

def mysys(t, x): # returns the RHS
    return ([x[1], -0.25*x[1]-np.sin(x[0])]);

tspan = [0,10]
x0 = 0; 
y0 = 5;
ics = [x0, y0]
sol = solve_ivp(mysys, tspan, ics, dense_output=True, rtol=1e-8);

tout = np.linspace(0, np.max(tspan), 100);
xout = sol.sol(tout)[0];
yout = sol.sol(tout)[1];
E = 1/2*yout**2 - np.cos(xout);
plt.plot(tout, xout, tout, yout)
#mlab.plot3d(np.cos(xout), np.sin(xout), yout )

[<matplotlib.lines.Line2D at 0x25275d28fd0>,
 <matplotlib.lines.Line2D at 0x25275d320d0>]

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

def mysys(t, r): # returns the RHS
    return ([r[0]*(1-r[0]**2), 1]);

tspan = [0,10]


x0 = 0; 
y0 = 2;
r0 = (x0**2 + y0**2)**0.5; t0 = np.arctan2(y0, x0);

ics = [r0, t0]
sol = solve_ivp(mysys, tspan, ics, dense_output=True, rtol=1e-8);

tout = np.linspace(0, np.max(tspan), 100);
rout = sol.sol(tout)[0];
tout = sol.sol(tout)[1];

xout = rout*np.cos(tout); yout = rout*np.sin(tout);

fig, ax = plt.subplots(2,1)
ax[0].plot(xout, yout);
ax[0].set_aspect(1);

ax[1].plot(tout, xout, tout, yout);

#mlab.plot3d(np.cos(xout), np.sin(xout), 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.

def show_traj(x0=0.9,y0=1.0, a = 0.0):
    x = np.linspace(-4, 4, 14); y = np.linspace(-4, 4, 14);
    X, Y = np.meshgrid(x, y);
    U = Y;
    V = -X-a*(X**2-1)*Y;
    plt.quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
    ax = plt.gca(); ax.set_aspect(1);

    def mysys(t, x): # returns the RHS
        return ([x[1], -x[0]-a*(x[0]**2-1)*x[1]]);

    tspan = [0,20]
    #x0 = -1.1; y0 = 2.5
    ics = [x0, y0]
    sol = solve_ivp(mysys, tspan, ics, dense_output=True, rtol=1e-8);

    tout = np.linspace(0, np.max(tspan), 200);
    xout = sol.sol(tout)[0];
    yout = sol.sol(tout)[1];
    
    plt.plot(xout,yout)
    #plt.plot(tout, xout, tout, yout)
w = interactive(show_traj, x0 = (-1.5, 1.5, 0.1), y0 = (-1.5, 1.5, 0.1), a = (-1.5, 1.5, 0.05))
w

• Let us also plot the trajectories

def show_traj(x0=0.9,y0=1.0, a = 0.0):
    x = np.linspace(-4, 4, 14); y = np.linspace(-4, 4, 14);
    X, Y = np.meshgrid(x, y);
    U = Y;
    V = -X-a*(X**2-1)*Y;
    #plt.quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
    #ax = plt.gca(); ax.set_aspect(1);

    def mysys(t, x): # returns the RHS
        return ([x[1], -x[0]-a*(x[0]**2-1)*x[1]]);

    tspan = [0,20]
    #x0 = -1.1; y0 = 2.5
    ics = [x0, y0]
    sol = solve_ivp(mysys, tspan, ics, dense_output=True, rtol=1e-8);

    tout = np.linspace(0, np.max(tspan), 200);
    xout = sol.sol(tout)[0];
    yout = sol.sol(tout)[1];
    
    #plt.plot(xout,yout)
    plt.plot(tout, xout, tout, yout)
w = interactive(show_traj, x0 = (-1.5, 1.5, 0.1), y0 = (-1.5, 1.5, 0.1), a = (-1.5, 1.5, 0.05))
w

• 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^2+b^2)$
• We plot the time series and the phase plot. Red circle denotes the fixed point and black circle denotes the initial point

def show_traj(x0=0.9, y0=1.0, a = 0.0, b = 0.1):
    x = np.linspace(0, 4, 14); y = np.linspace(0, 4, 14);
    X, Y = np.meshgrid(x, y);
    U = -X + a*Y + X**2*Y;
    V = b - a*Y - X**2*Y;

    def mysys(t, x): # returns the RHS
        return ([-x[0] + a*x[1] + x[0]**2*x[1], b - a*x[1] - x[0]**2*x[1]]);

    tspan = [0,40]
    #x0 = -1.1; y0 = 2.5
    ics = [x0, y0]
    sol = solve_ivp(mysys, tspan, ics, dense_output=True, rtol=1e-8);

    tout = np.linspace(0, np.max(tspan), 200);
    xout = sol.sol(tout)[0];
    yout = sol.sol(tout)[1];
    
    fig, ax = plt.subplots(1,2)
    ax[0].quiver(X, Y, U/(U**2+V**2)**0.5, V/(U**2+V**2)**0.5);
    ax[0].set_aspect(1)
    ax[0].plot(xout, yout)
    ax[0].plot(x0, y0, 'ok');
    ax[0].plot(b, b/(a+b**2), 'sr')
    ax[0].set_xlim(0,2.5)
    fig.set_size_inches(10, 5)
    
    ax[1].plot(tout, xout, tout, yout)
    
w = interactive(show_traj, x0 = (0, 1.5, 0.1), y0 = (0, 1.5, 0.1), a = (0, 0.2, 0.01), b = (0, 1.5, 0.01))
w