Computational Physics Chapter 3 problem 3.18，20

by defeng kong(2014301020081)

未分类

Abstract

Using the way and the code of problem 3.12,we can work out the problem and draw a conclusion that if the behavior is period n,Poincare sections will contain n discreter points.
As for 3.20,I just make a simple work.(I have a trouble)

Background

3.18
Plot $\omega$ versus $\theta$,with one point plotted for each drive cycle,as in Figure3.9,Do ihin for $F_{D}=1.4,1.44,$and$1.465$,using the other parameters as given in connection with Figure 3.10.
3.20
Calculate the bifurcation diagram for the pendulum in the vicinty of $F_{D}=1.35$to$1.5$.Make amagnified plot of the diagram (as compared to Figure 3.11)and obtain an estimate of the Fegenbaum $\delta$ parameter.

The GIF of problem 3.18,3.20:
for period-1

for period-2

for period-4

Main body

For 3.18：

I have written the code from problem 3.12.In problem 3.12,the $F_{D}=1.2$,For 3.18,we should do something for $F_{D}=1.4,1.44,$and$1.465$.Just append it.
The code is:

import matplotlib.pyplot as plt 
from math import sin,pi
g=9.8
l=9.8
class pendulums:
    def _init_(self,F,q,time_step=0.01,init_w=0,init_rad=0.2,dri=2.0/3):
        self.w=init_w
        self.rad=init_rad
        self.t=0
        self.dt=time_step
        self.dri=dri
        self.F=F
        self.q=q
        self.w1=[]
        self.rad1=[]
    def calculate(self):
        while(self.t<60000):
            self.w=self.w-((g/l)*sin(self.rad)\
            +self.q*self.w-self.F*sin(self.dri*self.t))*self.dt
            self.rad=self.rad+self.w*self.dt
            if self.rad>pi:
                self.rad=self.rad-2*pi
            elif self.rad<-pi:
                self.rad=self.rad+2*pi
            self.t=self.t+self.dt
            n=self.t*self.dri/(2.0*pi)
            if abs(n-round(n))<10e-6:
                self.w1.append(self.w)
                self.rad1.append( self.rad)
    def show_results1(self):
        plt.scatter(self.rad1,self.w1,s=12)
        plt.legend()
        plt.xlabel(r'$\omega$')
        plt.ylabel(r'$\theta$')
        plt.title(r'$\theta$ versus $\omega$ $F_D$=1.35')
        plt.show()
a= pendulums()
a._init_(1.40,0.5)
a.calculate()
a.show_results1() 
b= pendulums()
b._init_(1.44,0.5)
b.calculate()
b.show_results1() 
c= pendulums()
c._init_(1.465,0.5)
c.calculate()
c.show_results1() 
d= pendulums()
d._init_(1.2,0.5)
d.calculate()
d.show_results1() 

For 3.20

The code is

import matplotlib.pyplot as plt 
from math import sin,pi
g=9.8
l=9.8
class pendulums:
    def _init_(self,F,q,time_step=0.01,init_w=0,init_rad=0.2,dri=2.0/3):
        self.w=init_w
        self.rad=init_rad
        self.t=0
        self.dt=time_step
        self.dri=dri
        self.F=F
        self.q=q
        self.rad1=[]
        self.FD=[]
    def circulation(self,F1=1.35,F2=1.50,dF=0.001,init_w=0,init_rad=0.2):
        for i in range(int((F2-F1)/dF)):
            self.F=self.F+i*dF
            self.calculate()
            self.w=init_w
            self.rad=init_rad
            self.t=0
    def calculate(self):
        while(self.t<400*3*pi):
            self.w=self.w-((g/l)*sin(self.rad)\
            +self.q*self.w-self.F*sin(self.dri*self.t))*self.dt
            self.rad=self.rad+self.w*self.dt
            if self.rad>pi:
                self.rad=self.rad-2*pi
            elif self.rad<-pi:
                self.rad=self.rad+2*pi
            self.t=self.t+self.dt
            if self.t/3.0*pi>300:
                n=self.t*self.dri/(2*pi)
                if abs(n-round(n))<0.001:
                #if self.t%(3.0*pi)<10e-4:
                    self.FD.append(self.F)
                    self.rad1.append( self.rad)
    def show_results(self):
        plt.scatter(self.FD, self.rad1,s=2)
        plt.axis([1.35,1.5,1,3])
        plt.xlabel('$F_D$')
        plt.ylabel('angle$\theta$')
        plt.title('$\theta$ versus $F_D$')
        plt.show()
a= pendulums()
a._init_(1.35,0.5)
a.circulation()
a.calculate()
a.show_results()   

Conclusion

For 3.18

Form the Figure of FD=1.40,1.44,1465.I can find that in the period-1,it contain only one point.Likewise,if the behavior is period n,it will contain n discreter points.

For 3.20

in my code :

def circulation(self,F1=1.35,F2=1.50,dF=0.001,init_w=0,init_rad=0.2):
        for i in range(int((F2-F1)/dF)):
            self.F=self.F+i*dF

I set dt =0.001,but figure .... about dt=0.02,and the time of running is very very long.I don't think my code is uncorrect,or I can't find the uncorrect position.so I want to wait ,I think next week I can debug.

Thanks

myself