Exercise_8:Route to chaos: Period doubling

python physics 2014301020054 陈佩卓

Abstract

• Using Euler-Cromer method to solve the oscillatory motion.

• Also using class function,and take more factors into consideratioon.

• Find out how the value of $F_D$ influences the routes of the chaos.

Background

• We can use Euler method to solve some easy questions, but to solve the problem of the oscillatory motion and chaos with the conservation of energy, we cannot solve the problems easily by Euler method. we should use Euler-Cromer method to find the answer.

• we know the Euler-Cromer method. Now I will take the Euler-Cromer method into the oscillatory motion and chaos calculation:
  Newton’s second law for ideal pendulum with small angle:
  

  $$\frac{d^2\theta}{dx^2}=-\frac{l}{g}\theta$$

Write the second-order equations as two firest-order differential equations:

$$\frac{d\omega}{dt}=-\frac{g}{l}\omega$$

$$\frac{d\theta}{dt}=\omega$$

Finite difference form with Euler-Cromer method:

$$\omega_{i+1}=\omega_{i}-(\frac{g}{l})\theta_{i}\Delta{t}$$

$$\theta_{i+1}=\theta_{i}+\omega_{i+1}\Delta{t}$$

$$t_{i+1}=t_{i}+\Delta{t}$$

To make the question more realistic,we can take some factors into consideration.

1. We do not assume the small-angle approximation, and thus do not expand $sin\theta$ term in $\theta$

2. We include friction of the form

3. We add to our model a sinusoidal driving force $F_Dsin(\Omega_Dt)$

Putting all of these ingredients together, we have the equation of motion:

$$\frac{d^2\theta}{dt^2}=-\frac{g}{l}\theta-q\frac{d\theta}{dt}+F_Dsin(\Omega_Dt)$$
We can rewrite the two-order differential equation now as two first-order differential equations and obtian:
$$\frac{d\omega}{dt}=-\frac{g}{l}\theta-q\frac{d\theta}{dt}+F_Dsin(\Omega_Dt)$$

$$\frac{d\theta}{dt}=\omega$$
Finite difference form with Euler-Cromer method:
$$\omega_{i+1}=\omega_i-[(g/l)sin\theta_i-q\omega_i+F_Dsin(\Omega_Dt_i)]\Delta{t}$$

$$\theta_{i+1}=\theta_i+\omega_{i+1}\Delta{t}$$

$$t_{i+1}=t_i+\Delta{t}$$

If $\theta_{i+1}$ is out of the range $[-\pi,\pi]$, add or subtract $2\pi$ to keep it in this range.
And for the pioncare section, we chose the time that:

$$\omega_Dt=\theta+2n$$

Main Body

Question 3.18

Calculate Poincare sections for the pendulum as it undergoes the period-doubling route to chaos. Plot $\omega$ versus $\theta$, with one point plotted for each drive cycle, as in Figure 3.9. Do this for $F_D=1.4, 1.44$, and $1.465$, using the other parameters as given in connection with Figure 3.10. You should find that after removing the points corresponding to the initial transient the attractor in the period-1 regime will contain only a single point. Likewise, if the behavior is period $n$, the attractor will contain $n$ discrete points.

Codes

import pylab as pl
import math
class oscillatory:
    def __init__(self,g=9.8,l=9.8,q=0.5,F_D=1.2,O_D=2/3,time_step=3*math.pi/1000,\
    total_time=100,initial_theta=0.2,initial_omega=0):
        self.g=g
        self.l=l
        self.q=q
        self.F=F_D
        self.O=O_D
        self.t=[0]
        self.initial_theta=initial_theta
        self.initial_omega=initial_omega
        self.dt=time_step
        self.time= total_time
        self.omega= [initial_omega]
        self.theta= [initial_theta]
        self.nsteps=int(total_time//time_step+1)
        self.tmpo=[0]
        self.tmpt=[0]
    def run(self):
        for i in range(self.nsteps):
            tmpo=self.omega[i]+(-1*(self.g/self.l)*math.sin(self.theta[i])-\
            self.q*self.omega[i]+self.F*math.sin(self.O*self.t[i]))*self.dt
            tmpt=self.theta[i]+tmpo*self.dt
            while(tmpt<(-1*math.pi) or tmpt>math.pi):
                if tmpt<(-1*math.pi):
                   tmpt+=2*math.pi
                if tmpt>math.pi:
                   tmpt-=2*math.pi
            self.omega.append(tmpo)
            self.theta.append(tmpt)
            self.t.append(self.t[i]+self.dt)
        for i in range(len(self.t)):
            if self.t[i] // (3 * math.pi) > 300:
                if ((2 / 3 * self.t[i]) % (2 * math.pi)) <= (2 / 3 * self.dt * 0.5) or (2 * math.pi - ((2 / 3 * self.t[i]) % (2 * math.pi))) <= (2 / 3 * self.dt * 0.5):
                    self.tmpo.append(self.omega[i])
                    self.tmpt.append(self.theta[i])   
    def show_results(self):
        font = {'family': 'serif',
                'color':  'darkred',
                'weight': 'normal',
                'size': 16,}
        pl.plot(self.t,self.theta)
        #pl.scatter(self.tmpt, self.tmpo)
        pl.title(r'$\theta$ versus time', fontdict = font)
        pl.xlabel(r'$\theta$(radians)')
        pl.ylabel(r'$\omega$(rad/s)')
        #pl.xlim(0,3)
        #pl.ylim(-3,0)
        pl.legend((['$F_D$=1.2']))
        pl.show()
a = oscillatory()
a.run()
a.show_results()

• result:
  
  
  
  
  

If we chose the time :

$$\Omega_Dt=2n\pi+\theta$$

where n is integer and $\theta$ is a angle between $[-\pi,\pi]$.
We will have the result that:

Conclusion

1.We can solve the oscillatory motion by using Euler-Cromer method easily.We can see chaos when the problem become more complex,and from this problem,we know that the period of chaos is related to $F_D$.

2.we can easily find the period of the pendulum by choosing the time $\omega_Dt=\theta+2n$.

Acknowledgement

• Thank rf from hongyi.