Limit Cycle

Example: Van der Pol Oscillator

Van der Pol Oscillator is described by the equation

\begin{equation} \displaystyle x'' +\mu(x^2-1)x'+x=0 \tag{1} \end{equation}

where the parameter $\mu \geq 0$.

$\mu(x^2-1)x' $ is the nonlinear damping term. This term causes large amplitude oscillations to decay, but it also pumps them back up if thet become too small.

We can write equation (1) in two dimensional form by defining: $~~~x'=y$

\begin{eqnarray} x' &=& y \\ y'&=& \mu(1-x^2)y-x \end{eqnarray}

Next we solve this equation numerically to obtain the limit cycle.

In [1]:
# import all the packages required 
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
# Solve the differential equations and plot the phase portrait  
from scipy.integrate import odeint


def fun(X,t):
    mu=1.5 
    x=X[0];y=X[1];
    dx=y; dy=mu*(1-x*x)*y-x
    return [dx,dy]


plt.figure(figsize=(10,6))

# Times at which the solution is to be computed.
t = np.linspace(0, 20, 1500)




#let's plot some orbits/trajectories  
    
Xinit=np.array([0.01,0.01])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=1.5)


Xinit=np.array([-0.01,-0.01])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=1.5)

Xinit=np.array([1,-4])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=1.5)


Xinit=np.array([-1,4])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=1.5)

Xinit=np.array([-2,4])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=1.5)


Xinit=np.array([2,-8])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=1.5)




# let's plot the limit cycle in different color
Xinit=np.array([1.68,2])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(Xsolu[:,0], Xsolu[:,1],lw=5,color='k')


# let's plot the critical point as well 
plt.plot(0.0,0.0,'or',ms=12)


# stuff required for plotting arrows/vectors on the plot    
xmax = plt.xlim(xmin=0)[1];ymax = plt.ylim(ymin=0)[1]       

xl = np.linspace(-6.0,xmax,20)
yl = np.linspace(-8.0,ymax,20)
x1 , y1  = np.meshgrid(xl, yl) 
[u,v]= fun([x1,y1],0); 


# normalize the lengths of the vectors 
uvn=np.sqrt(u**2 + v**2); uvn[uvn==0]=1.
u=u/uvn; v=v/uvn;


plt.quiver(x1,y1,u,v,uvn,scale=50,cmap=plt.cm.gray)
plt.xlim(-3,2.5); plt.ylim(-9,7)   
plt.xlabel('$x(t)$',fontsize=20); plt.ylabel('$y(t)$',fontsize=20)

plt.show;
In [3]:
# Solve the differential equations and plot the time series


plt.figure(figsize=(12,3))


Xinit=np.array([0.01,0.01])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(t,Xsolu[:,0],lw=1.5)


Xinit=np.array([-0.01,-0.01])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(t,Xsolu[:,0],lw=1.5)

Xinit=np.array([1,-4])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(t,Xsolu[:,0],lw=1.5)


Xinit=np.array([-1,4])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(t,Xsolu[:,0],lw=1.5)

Xinit=np.array([-2,4])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(t,Xsolu[:,0],lw=1.5)


Xinit=np.array([2,-8])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(t,Xsolu[:,0],lw=1.5)




# let's plot the limit cycle in different color
Xinit=np.array([1.68,2])
Xsolu = odeint(fun,Xinit, t)    
plt.plot(t,Xsolu[:,0],lw=5,color='k')

    
    
plt.xlabel('$t$',fontsize=20); plt.ylabel('$x(t)$',fontsize=20)
plt.title('Time Series')
plt.show;
In [ ]: