import numpy as np
import matplotlib.pylab as plt
The main reason we are exploring this relationship instead of looking for the solution of the simple pendulum is that we are very used to things without discuss them.
Here we want to see how well the small angle approximation works. (Specially for this approximation).
The easiest way to understand the approximation is to compare the plots of the two functions we are approximating, (This particular case is also gotten from taking the 1st order Taylor expansion on zero of the $\sin(\theta)$ function).
theta=np.linspace(0,np.pi/2.,1000)
fig=plt.figure(figsize=(15,6))
ax1=fig.add_subplot(121)
ax2=fig.add_subplot(122)
ax1.plot(theta,theta,label='$\\theta$')
ax1.plot(theta,np.sin(theta),label='$\sin(\\theta)$')
ax1.legend()
ax1.set_title('Two functions plotted')
ax1.set_xlabel('$\\theta$ [rad]')
ax1.set_ylabel('Functions')
ax2.plot(theta,np.sin(theta)-theta, color='C2',label='Difference')
ax2.set_title('Difference')
ax2.set_xlabel('$\\theta$ [rad]')
ax2.set_ylabel('Difference $\sin(\\theta)-\\theta$')
ax2.legend()
fig.suptitle('Comparison $\\theta$ and $\sin(\\theta)$')
plt.show()
Graphically, one can see that they are almost the same for $\theta\approx 0.5$ rad, which correspond $\approx 30 $deg. (which at least for me, it is not very small!).
This data, are the result from a simulation of a dumped pendulum, we expect to have a exponentially decreasing oscilating function.
url1="https://raw.githubusercontent.com/jmsevillam/Herramientas-Computacionales-UniAndes/master/Data/dumped.dat"
Dumped=np.genfromtxt(url1).T
Time,theta1,omega1,theta2,omega2=Dumped
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\\theta_1$ vs $t$')
ax.set_ylabel('$\\theta$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Time,theta1)
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\\theta_2$ vs $t$')
ax.set_ylabel('$\\theta_2$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Time,theta2)
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\omega_1$ vs $t$')
ax.set_ylabel('$\\omega_1$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Time,omega1)
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\omega_2$ vs $t$')
ax.set_ylabel('$\\omega_2$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Time,omega2)
plt.show()
On this case the phase spaces look very similar!
fig=plt.figure(figsize=(8,8))
ax=fig.add_subplot(111)
ax.set_title('$\omega_1$ vs $\\theta_1$')
ax.set_ylabel('$\\omega_1$ [rad/sec]')
ax.set_xlabel('$\\theta_1$ [rad]')
ax.plot(theta2,omega2)
plt.show()
fig=plt.figure(figsize=(8,8))
ax=fig.add_subplot(111)
ax.set_title('$\omega_2$ vs $\\theta_2$')
ax.set_ylabel('$\\omega_2$ [rad/sec]')
ax.set_xlabel('$\\theta_2$ [rad]')
ax.plot(theta2,omega2)
plt.show()
It is common to have two different similar initial conditions and evaluate their difference over time in order to evaluate the chaos.
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$|\\theta_2-\\theta_2|$ vs $t$')
ax.set_ylabel('$|\\theta_2-\\theta_2|$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Time,np.abs(theta2-theta1))
plt.show()
As this relationship looks exponential, one can take el logarithm.
If this difference decrease the system is not chaotic (On this regime).
One can calculate the exponent of this behaviour, this exponent recieves the name of Lyapunov Exponent
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$|\\theta_2-\\theta_2|$ vs $t$')
ax.set_ylabel('$|\\theta_2-\\theta_2|$ [rad]')
ax.set_xlabel('time $t$')
ax.set_yscale("log")
ax.plot(Time,np.abs(theta2-theta1))
plt.show()
The case of having a strong (also oscillating) external force, makes the system change dramatically.
The following plots should be calculated from $(0,2\pi)$ or $(-\pi,\pi)$, but we are not going to worry about that.
url2="https://raw.githubusercontent.com/jmsevillam/Herramientas-Computacionales-UniAndes/master/Data/Driven.dat"
Driven=np.genfromtxt(url2).T
Timeb,theta1b,omega1b,theta2b,omega2b=Driven
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\\theta_1$ vs $t$')
ax.set_ylabel('$\\theta$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Timeb,theta1b)
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\\theta_2$ vs $t$')
ax.set_ylabel('$\\theta_2$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Timeb,theta2b)
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\omega_1$ vs $t$')
ax.set_ylabel('$\\omega_1$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Timeb,omega1b)
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\omega_2$ vs $t$')
ax.set_ylabel('$\\omega_2$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Timeb,omega2b)
plt.show()
On this phase spaces the trajectories cross each other because of the external force.
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\omega_1$ vs $\\theta_1$')
ax.set_ylabel('$\\omega_1$ [rad/sec]')
ax.set_xlabel('$\\theta_1$ [rad]')
ax.plot(theta2b,omega2b)
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\omega_2$ vs $\\theta_2$')
ax.set_ylabel('$\\omega_2$ [rad/sec]')
ax.set_xlabel('$\\theta_2$ [rad]')
ax.plot(theta2b,omega2b)
plt.show()
And the differences look way different.
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$|\\theta_2-\\theta1|$')
ax.set_ylabel('$\\theta_2$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Timeb,np.abs(theta2b-theta1b))
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\\theta_2$ vs $t$')
ax.set_ylabel('$\\theta_2$ [rad]')
ax.set_xlabel('time $t$')
ax.set_yscale("log")
ax.plot(Timeb,np.abs(theta2b-theta1b))
plt.show()
fig=plt.figure(figsize=(15,6))
ax1=fig.add_subplot(121)
ax2=fig.add_subplot(122)
ax1.set_yscale("log")
ax1.set_title('Two functions plotted')
ax1.set_xlabel('$\\theta$ [rad]')
ax1.set_ylabel('Functions')
ax1.plot(Time,np.abs(theta1-theta2),label='Difference')
ax1.legend()
ax2.set_yscale("log")
ax2.set_title('Difference')
ax2.set_xlabel('$\\theta$ [rad]')
ax2.set_ylabel('Difference $\sin(\\theta)-\\theta$')
ax2.plot(Timeb,np.abs(theta1b-theta2b),label='Difference')
ax2.legend()
fig.suptitle('Comparison $\\theta$ and $\sin(\\theta)$')
plt.show()
As angles are limited (from $0$ to $2\pi$ or from $-\pi$ to $\pi$), one can correct these previous result by transforming the data.
url2="https://raw.githubusercontent.com/jmsevillam/Herramientas-Computacionales-UniAndes/master/Data/Driven.dat"
Driven=np.genfromtxt(url2).T
Timeb,theta1b,omega1b,theta2b,omega2b=Driven
theta1b=(theta1b+np.pi)%(2*np.pi)
theta2b=(theta2b+np.pi)%(2*np.pi)
And repeating everything we get,
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\\theta_1$ vs $t$')
ax.set_ylabel('$\\theta$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Timeb,theta1b)
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\\theta_2$ vs $t$')
ax.set_ylabel('$\\theta_2$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Timeb,theta2b)
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\omega_1$ vs $t$')
ax.set_ylabel('$\\omega_1$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Timeb,omega1b)
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\omega_2$ vs $t$')
ax.set_ylabel('$\\omega_2$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Timeb,omega2b)
plt.show()
New visualization of the phase spaces (They are cylinders)
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\omega_1$ vs $\\theta_1$')
ax.set_ylabel('$\\omega_1$ [rad/sec]')
ax.set_xlabel('$\\theta_1$ [rad]')
ax.plot(theta2b,omega2b,'.')
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\omega_2$ vs $\\theta_2$')
ax.set_ylabel('$\\omega_2$ [rad/sec]')
ax.set_xlabel('$\\theta_2$ [rad]')
ax.plot(theta2b,omega2b,'.')
plt.show()
And more drastic behaviour
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\\theta_2-\\theta_1$ vs $t$')
ax.set_ylabel('$\\theta_2$ [rad]')
ax.set_xlabel('time $t$')
ax.plot(Timeb,np.abs(theta2b-theta1b))
plt.show()
fig=plt.figure(figsize=(10,6))
ax=fig.add_subplot(111)
ax.set_title('$\\theta_2$ vs $t$')
ax.set_ylabel('$\\theta_2$ [rad]')
ax.set_xlabel('time $t$')
ax.set_yscale("log")
ax.plot(Timeb,np.abs(theta2b-theta1b))
plt.show()
fig=plt.figure(figsize=(15,6))
ax1=fig.add_subplot(121)
ax2=fig.add_subplot(122)
ax1.set_yscale("log")
ax1.set_title('Non Chaotic')
ax1.set_xlabel('Time $t$')
ax1.set_ylabel('Difference $|\\theta_2-\\theta_1|$')
ax1.plot(Time,np.abs(theta1-theta2),label='Difference')
ax1.legend()
ax2.set_yscale("log")
ax2.set_title('Chaotic')
ax1.set_xlabel('Time $t$')
ax2.set_ylabel('Difference $|\\theta_2-\\theta_1|$')
ax2.plot(Timeb,np.abs(theta1b-theta2b),label='Difference')
ax2.legend()
fig.suptitle('Comparison $\\theta$ and $\sin(\\theta)$')
plt.show()
This simulation shows the intensity pattern of a plane light which travels towards a spherical mirror, as the mirror has a focus depending on its curvature, the intensity will be higher on that point
url3="https://raw.githubusercontent.com/jmsevillam/Herramientas-Computacionales-UniAndes/master/Data/LBWaves.dat"
Mirror=np.genfromtxt(url3)
fig=plt.figure(figsize=(10,7))
ax=fig.add_subplot(111)
ax.set_aspect("equal")
ax.set_title('Lattice Boltzann - Waves on a Mirror')
ax.set_ylabel('$y$')
ax.set_xlabel('$x$')
pl=ax.imshow(Mirror)
plt.colorbar(pl)
plt.show()
The amplituds are higher at (x=115,y=100). And the real solution is!
np.unravel_index(np.argmax(Mirror, axis=None), Mirror.shape)
We are going to explore how Taylor expansions work in the particular case of the sine function.
def fac(n):
if n==0 or n==1:
return 1
else:
return n*fac(n-1)
def sin_tay(n,x):
suma=np.zeros(len(x))
for i in range(n):
suma+=(-1.)**i/fac(2.*i+1)*x**(2.*i+1.)
return suma
x=np.linspace(-np.pi,np.pi,1000)
fig=plt.figure(figsize=(10,10))
for i in range(0,9):
ax=fig.add_subplot(3,3,i+1)
ax.plot(x,np.sin(x),label='$\sin(\\theta)$')
y=sin_tay(i,x)
index=2*i-1
if index<0: index=0
ax.plot(x,y,label='Taylor order '+str(index))
ax.set_ylim(-5,5)
ax.legend()
plt.show()
x=np.linspace(-3*np.pi,3*np.pi,1000)
fig=plt.figure(figsize=(10,10))
for i in range(0,9):
ax=fig.add_subplot(3,3,i+1)
ax.plot(x,np.sin(x),label='$\sin(\\theta)$')
y=sin_tay(i,x)
index=2*i-1
if index<0: index=0
ax.plot(x,y,label='Taylor order '+str(index))
ax.set_ylim(-5,5)
ax.legend()
plt.show()