第七次作业

本次作业包括习题3.12 3.13 3.14
3道题代码见后方

习题3.12

对于一个单摆，我们给予它一个周期性正弦驱动力。在考虑阻尼的情况下，做出了它的角速度-摆角（相空间）的图像。我们只选取其中驱动力位于一个周期的$\pi/4$位置处的相空间的点。最终结果如下图

习题3.13

此题研究的系统与上图相同。在此观察两个几乎相同的单摆，其唯一不同是初始摆角相差千分之一弧度。如下做出其摆角差（取对数）与时间的图像，会发现摆角差随时间大致沿指数增长。由此可以观察混沌现象

习题3.14

本题与上一题类似，也是研究两个几乎相同的单摆。但这次是将阻尼系数设置为一个微小的差值。同样做出其摆角差（取对数）与时间的图像。观察发现，一段时间后摆角差依然会大致沿指数增长。

代码如下

"""
by Cao Xinyu,2016/10/31
for exercises 3.12 3.13 3.14
"""
import pylab
from math import *
def resetangle(angle):
    angle+=pi
    angle%=pi*2
    angle-=pi
    return angle
class pendulum:
    def __init__(self,thet0,omega0,FD,length=9.8,q=0.5,omegaD=2/3,timestep=0.04,tmax=100):
        self.thet=[thet0]
        self.omega=[omega0]
        self.times=[0]
        self._FD=FD
        self._l=length
        self._dt=timestep
        self._q=q
        self._tmax=tmax
        self.g=9.8
        self._omegaD=omegaD
        self.calculate()
    def calculate(self):
        while self.times[-1]<self._tmax:
            self.omega.append(self.omega[-1]+self._dt*(-self.g/self._l*sin(self.thet[-1])-self._q*self.omega[-1]+self._FD*sin(self._omegaD*self.times[-1])))
            self.thet.append(self.thet[-1]+self.omega[-1]*self._dt)
            self.times.append(self.times[-1]+self._dt)
#3.12
a=pendulum(0.2,0,1.2,tmax=10000)
size=len(a.times)
lomega=[]
lthet=[]
for i in range(size):
    if abs((a.times[i]*2/3)%(2*pi)-pi/4)<=0.04*2/3 :
        lomega.append(a.omega[i])
        lthet.append(resetangle(a.thet[i]))
pylab.plot(lomega,lthet,'.')
pylab.xlabel('$\omega$ (rad/s)')
pylab.ylabel('$\Theta$ (rad)')# cannot theta?
pylab.xlim(min(lomega),max(lomega))
pylab.ylim(-pi,pi)
pylab.show()
#3.13
b=pendulum(0.201,0,1.2)
size=len(b.times)
llogdiffthet=[]
for i in range(size):
    llogdiffthet.append(log(abs(resetangle(a.thet[i]-b.thet[i]))))
pylab.plot(b.times,llogdiffthet,'.')
pylab.xlabel('time (s)')
pylab.ylabel('log($\Delta$ $\Theta$) (rad)')
pylab.show()
#3.14
c=pendulum(0.2,0,1.2,q=0.501)
ldiffthet2=[]
for i in range(size):
    ldiffthet2.append(abs(a.thet[i]-c.thet[i]))
pylab.plot(c.times,ldiffthet2,'.')
pylab.xlabel('time (s)')
pylab.ylabel('log($\Delta$ $\Theta$) (rad)')
pylab.show()
pass