4140.sagews

%auto
import numpy as np
import scipy as sp
from scipy import integrate

omega2 = 9.8/5; omega = sqrt(omega2)
tau = 2*np.pi/omega
phi0 = np.pi/9; Phi0 = np.array([phi0, 0])
dPhi_dt = lambda q,t: np.array([q[1], -omega2*sin(q[0])])
phi_app = lambda t: phi0*cos(omega*t)

t = np.arange(0, 3*tau, 0.025)
Phi = integrate.odeint(dPhi_dt, Phi0, t)
phi = Phi[:,0]

Lnum = line([(t[i],phi[i]) for i in range(len(t))], linestyle='-', marker='',rgbcolor=(0,0,0), legend_label="numerical")
Lapp = line([(t[i],phi_app(t[i])) for i in range(len(t))], linestyle='--', marker='',rgbcolor=(1,0,0),
    legend_label="approximate")
(Lnum+Lapp).show(axes_labels=("$t$", "$\\phi$"), legend_back_color=(1,1,1), legend_loc=4)

%auto
t = np.arange(0, 150, 0.025)
Phi = integrate.odeint(dPhi_dt, Phi0, t)
phi = Phi[:,0]
phi_app = lambda t, omega_eff: phi0*cos(omega_eff*t)

epsilon = 0.0076

indices = range(t.searchsorted(140), len(t))
Lnum = line([(t[i],phi[i]) for i in indices], linestyle='-', marker='',rgbcolor=(0,0,0), legend_label="numerical")
Lapp = line([(t[i],phi_app(t[i], (1.0-epsilon)*omega)) for i in indices], linestyle='--', marker='',
  rgbcolor=(1,0,0), legend_label="approximate")
html('$\epsilon\;=\;%s$'%latex(np.double(epsilon)))
(Lnum+Lapp).show(axes_labels=("$t$", "$\\phi$"), legend_back_color=(1,1,1), legend_loc=4)

\epsilon\;=\;0.0076