Exercise_08 : Period Doubling

Problem 3.18 & 3.19 & 3.20

冉峰 弘毅班 2014301020064

---

Abstract

• From the last work, we have used the method of the Euler-Cromer to solve the problem of the oscillatory motion and chaos. Now I want to have deeply consideration of the routes of the chaos. And we know that the routes of chaos have a relationship with the period doubling from the book. So I would like to use the bifurcation diagram and poincare section to do the consideration. And I will try the biger $F_D$ to do the work to understand the relationship between chaos and period doubling.

Background

• Now, 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}{dt^2}=-\frac{g}{l}\theta$$
  Write the second-order equations as two firest-order differential equations:
  $$\frac{d\omega}{dt}=-\frac{g}{l}\theta$$ $$\frac{d\theta}{dt}=\omega$$
  Finite difference form with Euler-Cromer method:
  $$\omega_{i+1}=\omega_i-(g/l)\theta_i\Delta t$$ $$\theta_{i+1}=\theta_i+\omega_{i+1}\Delta t$$ $$t_{i+1}=t_i+\Delta t$$
  Also to make the question a bit more realistic and more interesting we can take some factors into consideration.
  First, we do not assume the small-angle approximation, and thus do not expand $sin\theta$ term in $\theta$. Second, we include friction of the form $-q(d\theta/dt)$. Third, 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}sin\theta-q(d\theta/dt)+F_Dsin(\Omega_Dt)$$
  We can again rewrite the two-order differential equation as two first-order differential equations and otian:
  $$\frac{d\omega}{dt}=-\frac{g}{l}sin\theta-q(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.
  Also for the pioncare section, we chose the time that : $$\Omega_Dt=\theta+2n$$

---

The Main Body

Question 3.18 : Calculate Poincare section 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 work with F_D =1.35,1.4, 1.44, 1.465, using the other parameters as given in connection with Figure 3.10.

• code

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

It's easy to see that when $F_D=1.2$, the system is in chaotic state. When $F_D=1.35,F_D=1.4$, the system become well-aligned again, its period is the same as the drive period. When $F_D=1.44$, its period is the twice as the drive period. When $F_D=1.465$, its period is the fourth times as the drive period.

---

If we chose the time :

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

It's easy to see that when $F_D=1.35,F_D=1.4$, there is just one point,when $F_D=1.44$, there is two points, when $F_D=1.465$, there is four points.It's easy to understand this.

---

Question 3.19 : just like 3.18, we chose the time $$\Omega't=\theta+2n\pi$$ where $\Omega'$ is irrationally related to the $\Omega_D$.

• code

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] // (5 * math.pi) > 300:
                if ((2 / 5 * self.t[i]) % (2 * math.pi)) <= (2 / 3 * self.dt * 0.5) or (2 * math.pi - ((2/ 5 * self.t[i]) % (2 * math.pi))) <= (2 / 5 * 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

I chose the $\Omega_D'=1/5$, but unfortunately, I didn't get the right answer. Because I did'n have enough time to solve the bug, so it's not completely.

---

Vpython run

• code

from __future__ import division
from visual import *
from math import *
q=0.5
l=9.8
g=9.8
f=2/3
dt=0.04
theta0=0.2
omega0=0
ceil=box(pos=(0,5,0),size=(5,0.5,2),material=materials.bricks)
ball=sphere(pos=(l*sin(theta0),5-l*cos(theta0),0),radius=1,material=materials.emissive,color=color.green)
ball.omega=omega0
ball.theta=theta0
ball.t=0
ball.F=0.5
string=curve(pos=[ceil.pos,ball.pos],color=color.blue)
po=arrow(pos=ball.pos,axis=(9.8*ball.omega*cos(ball.theta),9.8*ball.omega*sin(ball.theta),0),color=color.blue)
while True:
    rate(100)
    ball.omega=ball.omega-((g/l)*sin(ball.theta)+q*ball.omega-ball.F*sin(f*ball.t))*dt
    angle_new=ball.theta+ball.omega*dt
    while angle_new>pi:
        angle_new-=2*pi
    while angle_new<-pi:
        angle_new+=2*pi
    ball.theta=angle_new
    ball.pos=(l*sin(ball.theta),5-l*cos(ball.theta),0)
    string.pos=[ceil.pos,ball.pos]
    po.pos=ball.pos
    po.axis=(9.8*ball.omega*cos(ball.theta),9.8*ball.omega*sin(ball.theta),0)
    ball.t=ball.t+dt

Result

$$F_D=0.5$$

$$F_D=1.2$$

$$F_D=1.35$$

$$F_D=1.465$$

---

Question 3.20 : Calculate the bifurcation diagrams for the pendulum in the vicinity of F_D=1.35 to 1.5. Make a magnified plot of the diagram and obtain an estimate of the Feigenbaum $\delta$ parameter.

• code

import math
import matplotlib.pyplot as pl
pl.rc('text', usetex=True)
pl.rc('font', family='serif')
class bifurcation_diagram:
    '''docstring for assignment this week, using the Runga-Kutta method, with reseting process'''
    def __init__(self, total_time, time_interval,driving_force_interval):
        self.FD = [1.35]
        self.dFD = driving_force_interval
        self.dt = time_interval
        self.steps = int(total_time // time_interval) + 1 
        self.t = [0]
        self.omega = [0]
        self.theta = [0.2]
        self.theta_in_phase = []
        self.theta_fixed = []
        while self.FD[-1] < 1.5:
            self.FD.append(self.FD[-1] + self.dFD)
    def calculate(self):
        for j in range(len(self.FD)):
            self.t = [0]
            self.omega = [0]
            self.theta = [0.2]
            self.omega_in_phase = []
            self.theta_in_phase = []
            for i in range(self.steps):
                midpoint_omega = self.omega[i] + 0.5 * (-math.sin(self.theta[i]) - 0.5 * self.omega[i] + self.FD[j] * math.sin(0.66666666667 * self.t[i])) * self.dt
                midpoint_time = self.t[i] + 0.5 * self.dt
                midpoint_theta = self.theta[i] + 0.5 * self.dt
                temporary_theta = self.theta[i] + midpoint_omega * self.dt
                temporary_omega = self.omega[i] + (-math.sin(midpoint_theta) - 0.5 * midpoint_omega + self.FD[j] * math.sin(0.66666666667 * midpoint_time)) * self.dt
                if math.fabs(temporary_theta) > math.pi :
                    if temporary_theta > 0 :
                        while temporary_theta > math.pi :
                            temporary_theta = temporary_theta - 2 * math.pi
                    else :  
                        while temporary_theta < -math.pi :
                            temporary_theta = temporary_theta  + 2 * math.pi
                self.theta.append(temporary_theta)
                self.omega.append(temporary_omega)
                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.theta_in_phase.append(self.theta[i])
            pl.plot([self.FD[j]] * len(self.theta_in_phase), self.theta_in_phase,'b.', linewidth = 0, alpha = 0.1)
            pl.xlim(1.35,1.5)
            pl.ylim(1,3)
            pl.xlabel("$F_D$")
            pl.ylabel(r'$\theta$(radians)')
            pl.title(r'Bifurcation diagram $\theta$ versus $F_D$')
        pl.show()
a1 = bifurcation_diagram(total_time = 4000, time_interval = 0.01, driving_force_interval = 0.001)
a1.calculate()

Result

---

Conclusion

• We cannot easily use Euler method to solve the oscillatory motion, but we can use Euler-Cromer method. When we take many factors into consideration, the oscillatory motion will be chaos from the figure we can see. From solving the problem about oscillatory motion, we have a better known about chaos. Also, with the help of computer and Euler-Cromer method we can solve the qusetion of pendulum easily. To have a better understand of the chaos, to do work on the period doubling with the bifuraction is very usefull. Also we can have a better understan
  of the chaos from the routes of chaos

---

Thanks For

• The book of Computational Physics, chapter 3
• The Lesson plan Chapter 3 of Cai Hao
• Thanks for the code of the vpython from ZZT