In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
%matplotlib inline
plt.rcParams.update({'font.size': 18})

In [2]:
# The linearized ODE
def ode1(X, t, gamma=0.1, omega=1.0):
x, y = X
dxdt = y
dydt = -2*gamma*y - (omega**2)*x
return [dxdt, dydt]

# The original (non-linear) ODE
def ode2(X, t, gamma=0.1, omega=1.0):
x, y = X
dxdt = y
dydt = -2*gamma*y - (omega**2)*np.sin(x)
return [dxdt, dydt]

om = 1.0
gam = 0.1
alpha = np.sqrt(om**2 - gam**2)
delta = np.arctan(-gam/alpha)
AA = (np.pi/2.0)/np.cos(delta)
print(alpha, delta, AA)

mypar = (gam, 1.0)
X0 = (np.pi/2.0, 0.0)
ts = np.linspace(0.0, 40.0, 500)
sol1 = odeint(ode1, X0, ts, args=mypar)
sol2 = odeint(ode2, X0, ts, args=mypar)
sol3 = AA*np.exp(-gam*ts)*np.cos(alpha*ts + delta)

plt.plot(ts, sol2[:,0], 'r-', label='numerical solution')
plt.plot(ts, sol3, 'b-.', label='linear solution')
#plt.tight_layout()
plt.xlabel('$t$')
plt.ylabel('$\\theta(t)$')
plt.gcf().set_size_inches(8,6)
plt.legend()
plt.savefig('ass2_1a.pdf')
plt.show()

plt.clf()
plt.gcf().set_size_inches(7,6)
plt.xlabel("$\\theta(t)$")
plt.ylabel("$\\eta(t)$")
xticks = [-np.pi/2.0,-np.pi/4.0,0,np.pi/4.0,np.pi/2.0]
xlabels = ['$-\\pi/2$', '$-\\pi/4$', '$0$', '$\\pi/4$', '$\\pi/2$']
plt.axes().set_aspect('equal')
plt.xticks(xticks, xlabels)
plt.yticks(xticks, xlabels)
plt.plot(sol1[:,0], sol1[:,1], 'b-.')
plt.plot(sol2[:,0], sol2[:,1], 'r-')
plt.tight_layout()
plt.savefig('ass2_1b.pdf')

0.994987437107 -0.100167421162 1.5787097085

In [3]:
def odesys(X, t, alpha=1.0, gamma=0.25, Omega=0.6667):
theta = X[0]
eta = X[1]
phi = X[2]
dtheta = eta
deta = -2*gamma*eta - np.sin(theta) + alpha*np.cos(phi)
dphi = Omega
return [dtheta, deta, dphi]

X0 = [0,0,0]
npts = 10000
tmax = 500.0
ts = np.linspace(0.0,1000.0,npts)

nmin = int(npts/2)

alphas = [0.9, 1.07, 1.15, 1.35]
fig, axes = plt.subplots(2,2)
fig.set_size_inches(10,7.5)
cols = ['r','b','g','m']
labels = ['(i)', '(ii)', '(iii)', '(iv)']
axs = axes.flatten()
for i in range(4):
Xs = odeint(odesys, X0, ts, args=(alphas[i],))
thetas = Xs[:,0]
thetas = np.mod(thetas + np.pi, 2*np.pi) - np.pi
etas = Xs[:,1]
ax = axs[i]
ax.cla()
ax.set_aspect(1)
ax.set_xlabel("$\\theta$")
ax.set_ylabel("$\\eta$")
ax.set_xticks([-np.pi,-np.pi/2,0,np.pi/2,np.pi])
ax.set_xticklabels(["$-\\pi$","$-\\pi/2$","$0$","$\\pi/2$","$\\pi$"])
ax.set_title(labels[i] + "   $\\alpha = " + repr(alphas[i]) + "$")
#xticks([-np.pi,-np.pi/2,0,np.pi/2,np.pi], )
ax.plot(thetas[nmin::], etas[nmin::], cols[i]+',')

fig.savefig("ass2_2a.pdf")

plt.clf()
Omega = 0.6667
T = np.pi / Omega
N = 1000
nmin = int(N/2)
t2s = np.arange(N)*T

cols = ['ro','b+','g.','mx']
sizes=[6,6,3,6]
for i in range(4):
Xs = odeint(odesys, X0, t2s, args=(alphas[i],))
thetas = np.mod(Xs[:,0] + np.pi, 2*np.pi) - np.pi
etas = Xs[:,1]
plt.plot(thetas[nmin::], etas[nmin::], cols[i], label="$\\alpha = " + repr(alphas[i]) + "$",ms=sizes[i])

plt.gcf().set_size_inches(8,6)
plt.xlabel("$\\theta$")
plt.ylabel("$\\eta$")
plt.xticks([-np.pi,-np.pi/2,0,np.pi/2,np.pi], ["$-\\pi$","$-\\pi/2$","$0$","$\\pi/2$","$\\pi$"])
#plt.set_aspect(1)
plt.legend()
fig.savefig("ass2_2b.pdf")

In [6]:
# Extensions 1: 3D plots
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D

X0 = [0,0,0]
npts = 10000
tmax = 500.0
ts = np.linspace(0.0,1000.0,npts)

nmin = int(npts/2)

alphas = [0.9, 1.07, 1.15, 1.35]
fig = plt.figure()
fig.set_size_inches(10,9)
cols = ['r','b','g','m']
labels = ['(i) $\\alpha = 0.9$', '(ii) $\\alpha = 1.07$', '(iii) $\\alpha = 1.15$', '(iv) $\\alpha = 1.35$']
for i in range(4):
Xs = odeint(odesys, X0, ts, args=(alphas[i],))
thetas = Xs[:,0]
thetas = np.mod(thetas + np.pi, 2*np.pi) - np.pi
etas = Xs[:,1]
phis = Xs[:,2]
ax.plot(thetas[nmin:], etas[nmin:], np.sin(phis[nmin:]), ',')
ax.set_xlabel("$\\theta$");
ax.set_ylabel("$\\eta$");
ax.set_zlabel("$\\sin\\varphi$");
white = (1.0, 1.0, 1.0, 1.0)
ax.w_xaxis.set_pane_color(white);
ax.w_yaxis.set_pane_color(white);
ax.w_zaxis.set_pane_color(white);
#ax.view_init(30,30)
ax.set_xticks((-np.pi/2,np.pi/2));
ax.set_yticks((-2,2));
ax.set_zticks((-1,0,1));
ax.set_title(labels[i])
plt.tight_layout()

plt.savefig("3d_phase.pdf")

In [7]:
# Extensions 1: 3D plots
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D

X0 = [0,0,0]
npts = 10000
tmax = 500.0
ts = np.linspace(0.0,1000.0,npts)

nmin = int(npts/2)

alp = 1.15
angs = [0,60,120,180]
fig = plt.figure()
fig.set_size_inches(10,9)
cols = ['r','b','g','m']
labels = ['(i) $\\alpha = 0.9$', '(ii) $\\alpha = 1.07$', '(iii) $\\alpha = 1.15$', '(iv) $\\alpha = 1.35$']
for i in range(4):
Xs = odeint(odesys, X0, ts, args=(alp,))
thetas = Xs[:,0]
thetas = np.mod(thetas + np.pi, 2*np.pi) - np.pi
etas = Xs[:,1]
phis = Xs[:,2]
ax.plot(thetas[nmin:], etas[nmin:], np.sin(phis[nmin:]), ',')
ax.set_xlabel("$\\theta$");
ax.set_ylabel("$\\eta$");
ax.set_zlabel("$\\sin\\varphi$");
white = (1.0, 1.0, 1.0, 1.0)
ax.w_xaxis.set_pane_color(white);
ax.w_yaxis.set_pane_color(white);
ax.w_zaxis.set_pane_color(white);
ax.view_init(30,angs[i])
ax.set_xticks((-np.pi/2,np.pi/2));
ax.set_yticks((-2,2));
ax.set_zticks((-1,0,1));
ax.set_title(labels[i])
plt.tight_layout()

plt.savefig("3d_phase.pdf")

In [5]:
# Extensions 2: Bifurcation diagram
t2s = np.arange(N)*T
nmin = N//2
apts = 200;  # apts = 100 takes one minute to run.
amin = 0.5; amax = 1.5;
alphas = np.linspace(amin, amax, apts)
ones = np.ones(N)
for i in range(apts):
Xs = odeint(odesys, X0, t2s, args=(alphas[i],))
thetas = np.mod(Xs[:,0] + np.pi, 2*np.pi) - np.pi
etas = Xs[:,1]
plt.plot(alphas[i]*ones[nmin:], thetas[nmin:], '.', ms=2)
plt.gcf().set_size_inches(10,8)
plt.xlabel("$\\alpha$")
plt.ylabel("$\\theta$")
plt.yticks([-np.pi,-np.pi/2,0,np.pi/2,np.pi], ["$-\\pi$","$-\\pi/2$","$0$","$\\pi/2$","$\\pi$"])
plt.savefig("bifurcation.jpg")
#plt.set_aspect(1)
#plt.legend()

In [41]:
# Extensions 3: Sensitivity to initial conditions
alp = 1.15
X0 = [0,0,0]
X1 = [1e-8,0,0]
ts = np.linspace(0,350,2500)
X0s = odeint(odesys, X0, ts, args=(alp,))
X1s = odeint(odesys, X1, ts, args=(alp,))
theta0s = np.mod(X0s[:,0] + np.pi, 2*np.pi) - np.pi
theta1s = np.mod(X1s[:,0] + np.pi, 2*np.pi) - np.pi
plt.plot(ts, theta0s, 'r-', ts,theta1s, 'b-')
plt.xlim(150, 350)
plt.xlabel("$t$")
plt.ylabel("$\\theta$")
plt.gcf().set_size_inches(8,6)
plt.savefig("determinism.pdf")

In [13]:
1e-8

Out[13]:
1e-08
In [ ]: