excercise 8

摘要
我们看到，在低驱动力的阻尼，非线性摆表现出简单的振荡运动，而在高驱动，它可以是混乱的，这提出了一个明显的问题，究竟是如何从简单到混沌行为的过渡发生？事实证明，钟摆过渡到混沌行为的几驱动力。

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 for F_D =1.4, 1.44, 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.

Part 1

import math
import matplotlib.pyplot as plt
g=9.8
dt=0.01

adjust theta to keep it in range of [-pi,+pi]

def adjust(x):
while x>math.pi:
x=x-2*math.pi
while x<-math.pi:
x=x+2*math.pi
return x

2-order Runge-Kutta method

def EC(omega0,theta0,q,l,FD,omegaD,T):
motion=[[]for i in range(3)]
motion[0].append(omega0)
motion[1].append(theta0)
motion[2].append(0)
while motion[2][-1]<=T:
motion[0].append(motion[0][-1]+(-g*math.sin(motion[1][-1])/l-q*motion[0][-1]+FD*math.sin(omegaD*motion[2][-1]))*dt)
motion[1].append(motion[1][-1]+motion[0][-1]*dt)
motion[2].append(motion[2][-1]+dt)
return motion#omega,theta,t
omegaD,T=0.66,1000
def Poin(motion):
poi=[[]for i in range(3)]
for n in range(int(omegaD*T/2/math.pi)):
number=int(round(2*n*math.pi/omegaD/dt))
poi[0].append(motion[0][number])
poi[1].append(motion[1][number])
poi[2].append(motion[2][number])
return poi

fig.2.3 poincare section

d=EC(0,0.2,0.5,9.8,1.2,0.66,T)
d1=Poin(d)
plt.plot(map(adjust,d1[1]),d1[0],'.')
plt.title('Fig.2.3 Poincare Section of the Physical Pendulum')
plt.xlabel(r'$\theta$ (radius)')
plt.ylabel(r'$\omega$ (radius/second)')
plt.text(0,0.7,r'$\omega_0=0,\theta_0=0.2,q=0.5,l=9.8,\Omega_D=0.66$')
plt.text(1.5,0.5,r'$F_D=1.2,T=1000s$')
plt.show()

fig.3.1.theta V.S. time with slightly different FD

d1=EC(0,0.2,0.5,9.8,1.35,0.66,100)#(omega0,theta0,q,l,FD,omegaD,T)
d2=EC(0,0.2,0.5,9.8,1.42,0.66,100)
d3=EC(0,0.2,0.5,9.8,1.44,0.66,100)
plt.subplot(1,3,1)#
plt.plot(d1[2],map(adjust,d1[1]))
plt.title(r'Fig.3.1 $\theta$ versus time')
plt.text(40,3.7,r'$F_D=1.35$')
plt.xlim(0,100)
plt.subplot(1,3,2)#
plt.plot(d2[2],map(adjust,d2[1]))
plt.title(r'Fig.3.2 $\theta$ versus time')
plt.text(40,3.7,r'$F_D=1.42$')
plt.xlim(0,100)
plt.subplot(1,3,3)#
plt.plot(d3[2],map(adjust,d3[1]))
plt.title(r'Fig.3.2 $\theta$ versus time')
plt.text(40,3.7,r'$F_D=1.44$')
plt.xlim(0,100)
plt.show()

fig 3.4 bifurcation diagram

def bif(motion,FD):
m=[[]for i in range(2)]#FD,theta
for j in range(30,60):
num=int(round(2*j*math.pi/omegaD/dt))
m[0].append(FD)
m[1].append(motion[1][num])
return m
for i in range(170):
F_D=1.3+i*0.001
d=EC(0,0.2,0.5,9.8,F_D,0.66,600)
b=bif(d,F_D)
plt.plot(b[0],map(adjust,b[1]),'.',markersize=3.0,color='k')
plt.title('Fig.3.4 Bifurcation Diagram')
plt.xlim(1.3,1.47)
plt.text(1.32,3,r'$\omega_0=0,\theta_0=0.2,q=0.5,l=g=9.8,\Omega_D=0.66$'+'\n30-60 drive periods')
plt.show()

part 2

import math
import pylab as py
g=9.8
l=9.8
q=0.5
ou_D=2/3
class pendulum0:
def init(self,F_D,sita):
self.sita=[sita]
self.omiga=[0]
self.t=[0]
self.dt=(3*(math.pi))/1000
self.F_D=F_D
self.omiga2=[]
self.sita2=[]
def calculate(self):
while self.t[-1]<6000:
omiga_new=self.omiga[-1]-((g/l)*(math.sin(self.sita[-1]))+q*self.omiga[-1]-self.F_D*math.sin(ou_D*self.t[-1]))*self.dt
self.omiga.append(omiga_new)
sita_new=self.sita[-1]+self.omiga[-1]self.dt
t_new=self.t[-1]+self.dt
while sita_new > math.pi:
sita_new -=2(math.pi)
while sita_new < -math.pi:
sita_new +=2*(math.pi)
self.sita.append(sita_new)
self.t.append(t_new)
sub1=py.subplot(221)
A=pendulum0(1.2,0.2)
A.calculate()
for i in range(300,600):
A.sita2.append(A.sita[1000*i])
A.omiga2.append(A.omiga[1000*i])
sub1.plot(A.sita2,A.omiga2,'r.')
sub1.set_title('$\omega$ versus $\\theta$')
sub1.legend(loc="upper right",frameon=False)
py.xlabel('$\\theta$(radians)')
py.ylabel('$\omega$(radians/s)')
sub1=py.subplot(222)
A=pendulum0(1.4,0.2)
A.calculate()
for i in range(300,600):
A.sita2.append(A.sita[1000*i])
A.omiga2.append(A.omiga[1000*i])
sub1.plot(A.sita2,A.omiga2,'r.')
sub1.set_title('$\omega$ versus $\\theta$')
sub1.legend(loc="upper right",frameon=False)
py.xlabel('$\\theta$(radians)')
py.ylabel('$\omega$(radians/s)')
sub1=py.subplot(223)
A=pendulum0(1.44,0.2)
A.calculate()
for i in range(300,600):
A.sita2.append(A.sita[1000*i])
A.omiga2.append(A.omiga[1000*i])
sub1.plot(A.sita2,A.omiga2,'r.')
sub1.set_title('$\omega$ versus $\\theta$')
sub1.legend(loc="upper right",frameon=False)
py.xlabel('$\\theta$(radians)')
py.ylabel('$\omega$(radians/s)')
sub1=py.subplot(224)
A=pendulum0(1.465,0.2)
A.calculate()
for i in range(300,600):
A.sita2.append(A.sita[1000*i])
A.omiga2.append(A.omiga[1000*i])
sub1.plot(A.sita2,A.omiga2,'r.')
sub1.set_title('$\omega$ versus $\\theta$')
sub1.legend(loc="upper right",frameon=False)
py.xlabel('$\\theta$(radians)')
py.ylabel('$\omega$(radians/s)')

part 3

from future import division
from visual import *
from math import *

add some constants

q=0.5
l=9.8
g=9.8
f=2/3
dt=0.04
theta0=0.2
omega0=0

add three ceilings

ceil1=box(pos=(-20,5,0),size=(5,0.5,2),material=materials.bricks)
ceil2=box(pos=(0,5,0),size=(5,0.5,2),material=materials.bricks)
ceil3=box(pos=(20,5,0),size=(5,0.5,2),material=materials.bricks)
ball1=sphere(pos=(-20+l*sin(theta0),5-l*cos(theta0),0),radius=1,material=materials.emissive,color=color.red)
ball2=sphere(pos=(l*sin(theta0),5-l*cos(theta0),0),radius=1,material=materials.emissive,color=color.yellow)
ball3=sphere(pos=(20+l*sin(theta0),5-l*cos(theta0),0),radius=1,material=materials.emissive,color=color.green)
ball1.omega=omega0
ball2.omega=omega0
ball3.omega=omega0
ball1.theta=theta0
ball2.theta=theta0
ball3.theta=theta0
ball1.t=0
ball2.t=0
ball3.t=0
ball1.F=0
ball2.F=0.5
ball3.F=1.2
string1=curve(pos=[ceil1.pos,ball1.pos],color=color.red)
string2=curve(pos=[ceil2.pos,ball2.pos],color=color.yellow)
string3=curve(pos=[ceil3.pos,ball3.pos],color=color.green)
po1=arrow(pos=ball1.pos,axis=(10*ball1.omega*cos(ball1.theta),10*ball1.omega*sin(ball1.theta),0),color=color.red)
po2=arrow(pos=ball2.pos,axis=(10*ball2.omega*cos(ball2.theta),10*ball2.omega*sin(ball2.theta),0),color=color.yellow)
po3=arrow(pos=ball3.pos,axis=(10*ball3.omega*cos(ball3.theta),10*ball3.omega*sin(ball3.theta),0),color=color.green)
while True:
rate(100)
ball1.omega=ball1.omega-((g/l)*sin(ball1.theta)+q*ball1.omega-ball1.F*sin(f*ball1.t))*dt
ball2.omega=ball2.omega-((g/l)*sin(ball2.theta)+q*ball2.omega-ball2.F*sin(f*ball2.t))*dt
ball3.omega=ball3.omega-((g/l)*sin(ball3.theta)+q*ball3.omega-ball3.F*sin(f*ball3.t))*dt
angle1_new=ball1.theta+ball1.omega*dt
while angle1_new>pi:
angle1_new-=2*pi
while angle1_new<-pi:
angle1_new+=2*pi
angle2_new=ball2.theta+ball2.omega*dt
while angle2_new>pi:
angle2_new-=2*pi
while angle2_new<-pi:
angle2_new+=2*pi
angle3_new=ball3.theta+ball3.omega*dt
while angle3_new>pi:
angle3_new-=2*pi
while angle3_new<-pi:
angle3_new+=2*pi
ball1.theta=angle1_new
ball2.theta=angle2_new
ball3.theta=angle3_new
ball1.pos=(-20+l*sin(ball1.theta),5-l*cos(ball1.theta),0)
ball2.pos=(l*sin(ball2.theta),5-l*cos(ball2.theta),0)
ball3.pos=(20+l*sin(ball3.theta),5-l*cos(ball3.theta),0)
string1.pos=[ceil1.pos,ball1.pos]
string2.pos=[ceil2.pos,ball2.pos]
string3.pos=[ceil3.pos,ball3.pos]
po1.pos=ball1.pos
po2.pos=ball2.pos
po3.pos=ball3.pos
po1.axis=(10*ball1.omega*cos(ball1.theta),10*ball1.omega*sin(ball1.theta),0)
po2.axis=(10*ball2.omega*cos(ball2.theta),10*ball2.omega*sin(ball2.theta),0)
po3.axis=(10*ball3.omega*cos(ball3.theta),10*ball3.omega*sin(ball3.theta),0)
ball1.t=ball1.t+dt
ball2.t=ball2.t+dt
ball3.t=ball3.t+dt

感谢郭思嘉同学对本人作业完成的帮助