Hw 8. Matplotlib

In [1]:
import numpy as np
import matplotlib.pylab as plt

Simple Pendulum!

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).

In [2]:
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!).

Dumped Pendulum

This data, are the result from a simulation of a dumped pendulum, we expect to have a exponentially decreasing oscilating function.

In [3]:
url1="https://raw.githubusercontent.com/jmsevillam/Herramientas-Computacionales-UniAndes/master/Data/dumped.dat"
In [4]:
Dumped=np.genfromtxt(url1).T
In [5]:
Time,theta1,omega1,theta2,omega2=Dumped
In [6]:
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()
In [7]:
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()
In [8]:
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()
In [9]:
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()

Phase spaces

On this case the phase spaces look very similar!

In [10]:
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()
In [11]:
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.

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

In [13]:
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()

Driven Pendulum

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.

In [14]:
url2="https://raw.githubusercontent.com/jmsevillam/Herramientas-Computacionales-UniAndes/master/Data/Driven.dat"
In [15]:
Driven=np.genfromtxt(url2).T
In [16]:
Timeb,theta1b,omega1b,theta2b,omega2b=Driven
In [17]:
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()
In [18]:
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()
In [19]:
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()
In [20]:
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()

Phase Spaces

On this phase spaces the trajectories cross each other because of the external force.

In [21]:
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()
In [22]:
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.

In [23]:
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()
In [24]:
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()
In [25]:
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()

Driven Pendulum***

As angles are limited (from $0$ to $2\pi$ or from $-\pi$ to $\pi$), one can correct these previous result by transforming the data.

In [26]:
url2="https://raw.githubusercontent.com/jmsevillam/Herramientas-Computacionales-UniAndes/master/Data/Driven.dat"
In [27]:
Driven=np.genfromtxt(url2).T
In [28]:
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,

In [29]:
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()
In [30]:
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()
In [31]:
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()
In [32]:
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()

Phase Spaces

New visualization of the phase spaces (They are cylinders)

In [33]:
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()
In [34]:
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

In [35]:
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()
In [36]:
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()
In [37]:
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()

Light on a Mirror

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

In [38]:
url3="https://raw.githubusercontent.com/jmsevillam/Herramientas-Computacionales-UniAndes/master/Data/LBWaves.dat"
In [39]:
Mirror=np.genfromtxt(url3)
In [40]:
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!

In [41]:
np.unravel_index(np.argmax(Mirror, axis=None), Mirror.shape)
Out[41]:
(100, 109)

Taylor Works

We are going to explore how Taylor expansions work in the particular case of the sine function.

In [42]:
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
In [43]:
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() 
In [44]:
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()