Computational Physics Chapter 3 problem3.12

by defeng kong(2014301020081)

未分类

摘要：

本文先得出F=1.2时混沌系统角速度随角度变化的图像即figure3.8右图，然后比较受力在最小值，和1/8周期处的庞加莱截面。

背景

3.12

In constructing the Poincare section in figure 3.9 we plotted points only at times that were in phase with the dirve force, that is,at time t=2pi n/D,where n is an interger. at these values of t the driving force passed through zero . however,we could just as easily have chosen to make the plot at times correaponding to a maximum of the drive force,or the time pi/8out-of-phase with this force,etc.construct the poincare sections for these cases and compare them with figure3.9

正文

代码为

import pylab as pl
from math import sin,pi
g=9.8
l=9.8
class pendulums:
    def _init_(self,F,q,time_step=0.04,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=[init_w]
        self.rad1=[init_rad]
    def calculate(self):
        while(self.t[-1]<400):
            self.w.append(self.w[-1]-((g/l)*sin(self.rad[-1])\
            +self.q*self.w[-1]-self.F*sin(self.dri*self.t[-1]))*self.dt)
            self.rad.append(self.rad[-1]+self.w[-1]*self.dt)
            if self.rad[-1]>pi:
                self.rad[-1]=self.rad[-1]-2*pi
            elif self.rad[-1]<-pi:
                self.rad[-1]=self.rad[-1]+2*pi
            self.t.append(self.t[-1]+self.dt)
            n=self.t[-1]*self.dri/(2*pi)-0.5/2
            if abs(n-round(n))<0.001:
                self.w1.append(self.w[-1]-((g/l)*sin(self.rad[-1])\
            +self.q*self.w[-1]-self.F*sin(self.dri*self.t[-1]))*self.dt)
                self.rad1.append(self.rad[-1]+self.w[-1]*self.dt)
    def show_results(self):
        pl.plot(self.rad, self.w)
        #pl.axis([0,60,-3,3])
        pl.axis([-4,4,-3,3])
        pl.show()
    def show_results(self):
        pl.plot(self.rad1, self.w1)
        #pl.axis([0,60,-3,3])
        pl.axis([-4,4,-3,3])
        pl.show()
a= pendulums()
a._init_(1.2,0.5)
a.calculate()
a.show_results()
a= pendulums()
b._init_(1.2,0.5)
b.calculate()
b.show_results1()

F=1.2时角速度随角度变化图像（运行时间为400s的效果图）：


上图中是受力在1/8周期和受力最小处的比较（绿线：1/8，蓝线：受力最小）