# Limit Cycle¶

### Example: Van der Pol Oscillator¶

Van der Pol Oscillator is described by the equation

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

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 [ ]: