Pendulum¶
This notebook shows some analysis of one of the most used systems, a pendulum.
We are going to do some simple simulations of a pendulum under different conditions that will lead to very important concepts. We are going to focus on the chaotic behavior, so Lyapunov exponets, Poincaré sections and bifurcation maps are going to be builded.
The system that we are going to use, is described by the following equation of motion
$$ \frac{d^2\theta}{dt^2}-\gamma\frac{d\theta}{dt}+\omega_0^2\sin(\theta)=F\sin(\Omega t) $$We have choosen this particular form of the external force so that some interesting behaviour emerges from the frequency $\Omega$.
For this particular regard, we are not going to use a sofisticated integrator, just going to use the Euler Cromer method which is very similar to the Euler method but has the huge advantage to be symplectic which means that it conserves energy.
Code¶
The code here presented is aimed to be as simple as possible to be able to do the simulations so that the post analysis is donde in this notebook, so we may first import the libraries to be used,
#include <iostream> //To print
#include <cstdlib> //To Convert Datatypes
#include <cmath> //To Calculate Mathematical Functions
There are some initial conditions and constants that are going to be fixed during the simulation.
This is a first attempt to do the study of this system, and therefore the fixed quantities are going to change when looking different problems. So far we are just getting familiar with the problem, so
const double omega_f=2./3.;
const double g=9.8;
const double l=9.8;
const double gamma0=0.5;
const double dt=0.001;
const double q0=1;
const double p0=0.0;
We now define a timestep function which basically solves the equation of motion. The equations of motion are those defined bellow.
void timestep(double & q, double & p, double& t, double F);
double dqdt(double q, double p, double t,double F);
double dpdt(double q, double p, double t, double F);
The main function is actually quite easy, because it only recieves the external parameters, sets up the initial conditions and evolve in time while evolves the system.
int main(int argc, char *argv[]) {
double Force=atof(argv[1]);
double q1=q0,p1=p0;
double t1=0.0;
double t2=0.0;
int N=atoi(argv[2]);
for(int i=0; i<N;i++){
std::cout<<t1<<' '<<q1<<' '<<p1<<std::endl;
timestep(q1,p1,t1,Force);
}
return 0;
}
As mentioned before, we are using the Euler Cromer algorithm to do the integration, so
void timestep(double & q, double & p, double& t, double F){
p+=dpdt(q,p,t,F)*dt;
q+=dqdt(q,p,t,F)*dt;
t=t+dt;
}
Now, we have the equations of motion for the system,
double dpdt(double q, double p, double t, double F_f){
return -(g/l)*sin(q)-gamma0*dqdt(q,p,t,F_f)+F_f*sin(omega_f*t);
}
double dqdt(double q, double p, double t, double F_f){
return p;
}
Trajectories on Phase Space and Lyapunov Exponent¶
Let us run this code two times, the second one with a slightly different initial conditions, so that we want to evaluate the evolution of this difference with time, this will require a small change in our code, or simply running it twice changing the initial conditions,
import numpy as np
import matplotlib.pylab as plt
We already implemented this change, so that we use as an external force amplitud $F=0.5$ in the units established by the other constant of the system, and for $200000$ timesteps. The initial difference in the angle is set to be $\Delta\theta=0.1$,
%%bash
g++ pendulum.cpp
time ./a.out 0.5 200000 0.1 > data_pend.dat
If we choose one of the two pendulums and plot the phase space, it would look like,
data=np.genfromtxt("data_pend.dat")
plt.plot((np.pi+data[:,1])%(2*np.pi)-np.pi,data[:,2],'.k', markersize=0.1)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
plt.savefig("test.png")
So if we take the difference between the two angles in time we may get something like
plt.yscale("log")
plt.plot((data[:,0]),np.abs(data[:,1]-data[:,3]),'.k', markersize=0.1)
plt.xlabel("Time $t$")
plt.ylabel("$\log|\\theta_1-\\theta_2|$")
plt.savefig("test1.png")
Or similarly with the frequence
plt.yscale("log")
plt.plot((data[:,0]),np.abs(data[:,2]-data[:,4]),'k.', markersize=0.1)
plt.xlabel("Time $t$")
plt.ylabel("$\log|\omega_1-\omega_2|$")
plt.savefig("test2.png")
Look that the $y$ axis is in logarithmic scale, wich basically means that the two pendulums are getting exponentially closer as time passes.
Let us consider the same situation, but with a very strong external force $F=1.4$
%%bash
g++ pendulum.cpp
time ./a.out 1.4 200000 0.1 > data_pend14.dat
The phase space now looks different, but still very regular
data=np.genfromtxt("data_pend14.dat")
plt.plot((np.pi+data[:,1])%(2*np.pi)-np.pi,data[:,2],'.k', markersize=0.1)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
plt.show()
The difference evolvution in time now is not that regular as we had before, but still it decreases almost exponential until a stabilization
plt.yscale("log")
plt.plot((data[:,0]),np.abs(data[:,1]-data[:,3]),'.k', markersize=0.1)
plt.xlabel("Time $t$")
plt.ylabel("$\log|\\theta_1-\\theta_2|$")
plt.show()
plt.savefig("1-4.png")
Exactly the same process we just did, but now with $F=1.44$
%%bash
g++ pendulum.cpp
time ./a.out 1.44 200000 0.1 > data_pend144.dat
data=np.genfromtxt("data_pend144.dat")
plt.plot((np.pi+data[:,1])%(2*np.pi)-np.pi,data[:,2],'.k', markersize=0.1)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
plt.show()
Even when the phase space look very similar, we get that the behavior of the difference of the angles now changes a lot, as it gets stabilized sooner and with larger fluctiations
plt.yscale("log")
plt.plot((data[:,0]),np.abs(data[:,1]-data[:,3]),'.k', markersize=0.1)
plt.xlabel("Time $t$")
plt.ylabel("$\log|\\theta_1-\\theta_2|$")
plt.show()
So, now let us increase one more time the amplitud of the external force, so that
%%bash
g++ pendulum.cpp
time ./a.out 1.465 200000 0.1 > data_pend1465.dat
Where the phase space may look something like,
data=np.genfromtxt("data_pend1465.dat")
plt.plot((np.pi+data[:,1])%(2*np.pi)-np.pi,data[:,2],'.k', markersize=0.1)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
plt.show()
Look that the behavior is not very different, so we may expect to have also a similar form on the difference plot,
plt.yscale("log")
plt.plot((data[:,0]),np.abs(data[:,1]-data[:,3]),'.k', markersize=0.1)
plt.xlabel("Time $t$")
plt.ylabel("$\log|\\theta_1-\\theta_2|$")
plt.show()
Now the system looks to increase the difference until a stabilization. We will explore more this in a moment. This makes us wonder What happens if we increas even more the force? so, let us try it.
%%bash
g++ pendulum.cpp
time ./a.out 1.5 200000 0.1 > data_pend15.dat
The behavior now is very different on the phase space plot,
data=np.genfromtxt("data_pend15.dat")
plt.plot((np.pi+data[:,1])%(2*np.pi)-np.pi,data[:,2],'.k', markersize=0.1)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
plt.show()
plt.yscale("log")
plt.plot((data[:,0]),np.abs(data[:,1]-data[:,3]),'.k', markersize=0.1)
plt.xlabel("Time $t$")
plt.ylabel("$\log|\\theta_1-\\theta_2|$")
plt.show()
The plot of the logarithm of the difference of the initial conditions is in fact, very important because it is related with the Chaoticity of the system. The
$$ \Delta \theta \sim \exp(-\lambda t) $$Where $\lambda$ is the first Lyapunov exponent. The sign of this exponent tells if the system is or not chaotic.
This is a preliminary inspection on the system which has lead to some ideas that we can summarize now,
- In this system the external force can change the system from regular, to chaotic.
- The Lyapunov exponent tells easily if the system is or not chaotic.
- The phase space sometimes cannot be enough to know if the system is chaotic.
But some questions still are there, for instance, The harder the external force is, the more chaotic the system gets?.
The answer to this is no, and let us explore that.
Consider an amplitud for the external force as $F=1.2$
%%bash
g++ pendulum.cpp
time ./a.out 1.2 200000 0.001 > data_pend_chaos.dat
data=np.genfromtxt("data_pend_chaos.dat")
So that the phase space looks like
plt.plot((np.pi+data[:,1])%(2*np.pi),data[:,2],'k.', markersize=0.1)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
And the plot of the difference,
plt.yscale("log")
plt.plot((data[:,0]),np.abs(data[:,1]-data[:,3]),'k.', markersize=0.1)
plt.xlabel("Time $t$")
plt.ylabel("$\log|\\theta_1-\\theta_2|$")
And the same happen for the frequency
plt.yscale("log")
plt.plot((data[:,0]),np.abs(data[:,2]-data[:,4]),'k.', markersize=0.1)
plt.xlabel("Time $t$")
plt.ylabel("$\log|\omega_1-\omega_2|$")
So it is not necessary to have a huge force to make the system chaotic, and it looks like the system can go from chaotic to regular and then chaotic again when increasing a parameter, which in this case is the amplitud of the external force.
Let us repeat the phase space plots, but now with more points to see if more clear the differences on the dynamics,
%%bash
g++ pendulum1.cpp
time ./a.out 0.5 2000000 0.001 > data_poinc.dat
data=np.genfromtxt("data_poinc.dat")
plt.plot((np.pi+data[:,1])%(2*np.pi),data[:,2],'k.', markersize=0.01)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
plt.show()
%%bash
g++ pendulum1.cpp
time ./a.out 1.2 2000000 0.001 > data_poinc.dat
data=np.genfromtxt("data_poinc.dat")
plt.plot((np.pi+data[:,1])%(2*np.pi),data[:,2],'k.', markersize=0.01)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
plt.show()
%%bash
g++ pendulum1.cpp
time ./a.out 1.5 2000000 0.001 > data_poinc.dat
data=np.genfromtxt("data_poinc.dat")
plt.plot((np.pi+data[:,1])%(2*np.pi),data[:,2],'k.', markersize=0.01)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
plt.show()
%%bash
g++ pendulum1.cpp
time ./a.out 1.4 2000000 0.001 > data_poinc.dat
data=np.genfromtxt("data_poinc.dat")
plt.plot((np.pi+data[:,1])%(2*np.pi),data[:,2],'k.', markersize=0.01)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
plt.show()
%%bash
g++ pendulum1.cpp
time ./a.out 1.44 2000000 0.001 > data_poinc.dat
data=np.genfromtxt("data_poinc.dat")
plt.plot((np.pi+data[:,1])%(2*np.pi),data[:,2],'k.', markersize=0.01)
plt.xlabel("Angle $\\theta$")
plt.ylabel("Frequency $\omega$")
plt.show()
Here we can see that the system tends to go in two paths, so it looks like the single line of the previous combination of parateters divided into two. Let us explore this more detailed.
Bifurcation plots¶
So we want to know when the system starts being chaotic and when it finishes, to do so we are about to build the bifurcation plots which are basically ploting the final points of a long simulation in terms of the external parameters.
We say a long simulation because we want to remove the transient behavior and just use the steady solution of the system, so this long will change a lot deppending on the system of interest.
We are going to take a single point of each simulation, so we select randomly a Force on the interval desired and evolve it for a large time and then plot just one point! (In fact more than one single point can be taken, but this would generate vertical fringes that we want to avoid, so we are going to have a simulation way heavier to get a densed populated plot)
%%bash
g++ pendulum2.cpp
time ./a.out 1 30000 50000 > data_bif.dat
If we plot this, we get two different behaviors, first one or a small number of lines, which means regular behavior, and regions where the points looks desordered and noisy, those regions are the chaotic ones
data= np.genfromtxt("data_bif.dat").T
q=data[1]
F=data[-1]
plt.plot(F,np.mod(q+np.pi,2*np.pi)-np.pi,'.k',markersize=0.1)
plt.xlabel("Fuerza $F$")
plt.ylabel("Ángulo $\\theta$")
plt.savefig("bif1.png")
This can be also seen on the frequencies,
p=data[2]
F=data[-1]
plt.plot(F,p,'.k',markersize=0.1)
plt.xlabel("Fuerza $F$")
plt.ylabel("Frecuencia $\omega$")
plt.savefig("bif2.png")
So if we repeat this procedure but restricted to a smaller region
%%bash
g++ pendulum2a.cpp
time ./a.out 1 30000 50000 > data_bif.dat
data= np.genfromtxt("data_bif.dat").T
q=data[1]
F=data[-1]
plt.plot(F,np.mod(q+np.pi,2*np.pi)-np.pi,'.k',markersize=0.1)
plt.xlabel("Force $F$")
plt.ylabel("Angle $\\theta$")
plt.show()
Now we can see why we got the chaotic behavior for $F=1.2$ and not for $F=1.4$.
Again, this can also been seen on the frequencies,
p=data[2]
F=data[-1]
plt.plot(F,p,'k.',markersize=0.2)
plt.xlabel("Force $F$")
plt.ylabel("Frequency $\omega$")
plt.show()
The key point here is that chaotic systems are ergodic but this is going to be presented with more detail later.
Look that before going chaotic the plot divides into two, that is why these are called Bifurcation plots.
But the external force is not the only important parameter that can change the behavior of the system, for example if we consider that the friction can be so hard that the external force just get shuted down, there wont be chaotic behavior, It can be seen on the following animation
Viscosity¶
So we can also leave fixed the external force and change the viscocity in order to get chaotic or regular behavior,
%%bash
g++ pendulum4.cpp
time ./a.out 1 50000 50000 > data_bif_damp.dat
Now, we started with a chaotic behavior, that leads to regular one, as the effect of increasing the viscosity is opposite to increase the force,
data= np.genfromtxt("data_bif_damp.dat").T
q=data[1]
F=data[-1]
plt.plot(F,np.mod(q+np.pi,2*np.pi)-np.pi,'.k',markersize=0.1)
plt.xlabel("Viscosity $\gamma$")
plt.ylabel("Angle $\\theta$")
plt.show()
And as usual, the frequencies also show the chaotic behavior,
p=data[2]
F=data[-1]
plt.plot(F,p,'.k',markersize=0.1)
plt.xlabel("Viscosity $\gamma$")
plt.ylabel("Frequency $\omega$")
plt.show()
Frequency of external force¶
There is a last parameter to change on this model, which is the frequency of the external force, as resonances can occur this parameter can lead to a huge impact
%%bash
g++ pendulum5.cpp
time ./a.out 1 50000 50000 > data_bif_frec.dat
data= np.genfromtxt("data_bif_frec.dat").T
q=data[1]
F=data[-1]
plt.plot(F,np.mod(q+np.pi,2*np.pi)-np.pi,'.k',markersize=0.1)
plt.xlabel("Frquency External Force $\Omega$")
plt.ylabel("Angle $\\theta$")
plt.show()
p=data[2]
F=data[-1]
plt.ylim(top=1)
plt.plot(F,p,'.k',markersize=0.1)
plt.xlabel("Frquency External Force $\Omega$")
plt.ylabel("Frquency $\omega$")
plt.show()
So we shouldc construct the complete space of parameters and see in which parts the combination of parameteres lead to chaotic behavior. This is of course computationally very expensive and very impracticle as this space may have more than 3 dimensions. Even in this case of a one dimensional problem, this study can be as complicated as we want given the freedom to do many things, but as this is a preliminar and introductory work, we are going to move to a different strategy to study chaotic systems which is very common for periodic systems, the stroboscopic maps.