第七次作业

1.摘要

本次作业主要是研究非线性单摆在加入阻尼和驱动力等因素下产生的混沌现象，并在相空间里对其分析与研究。

2.背景

混沌理论：混沌理论(Chaos theory)是关于非线性系统在一定参数条件下展现分岔（bifurcation）、周期运动与非周期运动相互纠缠，以至于通向某种非周期有序运动的理论。在耗散系统和保守系统中，混沌运动有不同表现，前者有吸引子，后者无（也称含混吸引子）。
对于一个单摆，考虑阻尼和驱动力，并且其不再做小角度运动，由此我们可以得到：

$$\frac{\mathrm{d^2} \theta }{\mathrm{d} x^2}=-\frac{g}{l} \sin \theta -q\frac{\mathrm{d} \theta }{\mathrm{d}t}+F_{D}\sin (\Omega _{D}t)$$
我们可以将上式花间为一阶常微分方程：
$$\frac{\mathrm{d} \omega }{\mathrm{d} t}=-\frac{g}{l}\sin \theta -q\frac{\mathrm{d} \theta }{\mathrm{d} t}+F_{D}\sin (\Omega _{D}t)$$
$$\frac{\mathrm{d} \theta }{\mathrm{d} t}=\omega$$
我们用Euler-Cromer method进行计算：
$$\omega _{i+1}=\omega _{i}-[(g/l)\sin \theta _{i}-q\omega _{i}+F_{D}\sin (\Omega _{D}t)]\Delta t$$
$$\theta _{i+1}=\theta _{i}+\omega _{i+1}\Delta t$$
并且如果$\theta _{i+1}$在【$-\pi ，\pi$】之外，则$\theta _{i+1}$需加上或减去2$\pi$

3.正文

import math
import pylab as pl
#初始化各个值，以及计算单摆的运功轨迹，展示结果曲线图像
class pendulums_basis:
    def __init__(self, F_D = 1.2, total_time = 60, init_theta = 0.2, q = 0.5, l = 9.8, g = 9.8, Omega_D = 2/3, time_step = 0.04):
        self.omega = [0]
        self.theta = [init_theta]
        self.t = [0]
        self.dt = time_step
        self.total_time = total_time
        self.q = q
        self.l = l
        self.g = g
        self.Omega_D = Omega_D
        self.F_D = F_D
    def calculate(self):
        i = 0
        loop_calculate = True
        while(loop_calculate):
            self.omega.append(self.omega[i] + ((-self.g / self.l) * math.sin(self.theta[i]) - self.q * self.omega[i] + self.F_D * math.sin(self.Omega_D * self.t[i])) * self.dt)
            self.theta.append(self.theta[i] + self.omega[i + 1] * self.dt)
            self.t.append(self.t[i] + self.dt)
            if self.theta[i + 1] < -math.pi:
                self.theta[i + 1] = self.theta[i + 1] + 2 * math.pi
            if self.theta[i + 1] > math.pi:
                self.theta[i + 1] = self.theta[i + 1] - 2 * math.pi
            i += 1
            if self.t[i] > self.total_time:
                loop_calculate = False
    def show_results(self):
        pl.plot(self.t, self.theta)
        pl.title("$\\theta$ versus $t$   $F_D=%.1f$"%self.F_D, fontsize=20)
        pl.xlabel('time(s)', fontsize=20)
        pl.ylabel('$\\theta$(radians)', fontsize=20)
        pl.show()
#计算两个几乎相同的摆，只是初始的$\theta$有微小的差异，每个时间点对应的$\Delta t=|\theta_1-\theta_2|$与时间$t$的关系曲线
class pendulums_comparison(pendulums_basis):
    def compare(self):
        pendulum_1 = pendulums_comparison(self.F_D, self.total_time, self.theta[0], 0.5)
        pendulum_1.calculate()
        pendulum_2 = pendulums_comparison(self.F_D, self.total_time, self.theta[0], 0.501)
        pendulum_2.calculate()
        dtheta = []
        loop_compare = True
        i = 0
        while(loop_compare):
            dtheta.append(abs(pendulum_1.theta[i] - pendulum_2.theta[i]))
            i += 1
            if i == len(pendulum_1.t):
                loop_compare = False
        pl.semilogy(pendulum_1.t, dtheta)
        pl.title("$\\Delta\\theta$ versus $t$   $F_D=%.1f$"%self.F_D, fontsize=20)
        pl.xlabel('time(s)', fontsize=20)
        pl.ylabel('$\\Delta\\theta$(radians)', fontsize=20)
        pl.show()
#展示$\omega$与$\theta$的关系曲线
class omega_versus_theta(pendulums_basis):
    def show_results(self):
        pl.plot(self.theta, self.omega, '.',)
        pl.title("$\\omega$ versus $\\theta$   $F_D=%.1f$"%self.F_D, fontsize=20)
        pl.xlabel('$\\theta$(radians)', fontsize=20)
        pl.ylabel('$\\omega$(radians/s)', fontsize=20)
        pl.show()
#展示在满足一定相位条件的$t$下$\omega$与$\theta$的关系曲线
class poincare_section(pendulums_basis):
    def calculate(self):
        i = 0
        n = 1
        self.ps_omega = []
        self.ps_theta = []
        loop_calculate = True
        while(loop_calculate):
            self.omega.append(self.omega[i] + ((-self.g / self.l) * math.sin(self.theta[i]) - self.q * self.omega[i] + self.F_D * math.sin(self.Omega_D * self.t[i])) * self.dt)
            self.theta.append(self.theta[i] + self.omega[i + 1] * self.dt)
            self.t.append(self.t[i] + self.dt)
            if self.theta[i + 1] < -math.pi:
                self.theta[i + 1] = self.theta[i + 1] + 2 * math.pi
            if self.theta[i + 1] > math.pi:
                self.theta[i + 1] = self.theta[i + 1] - 2 * math.pi
            #$\Omega_Dt=2n\pi$
            if  abs(self.t[i + 1] - 2 * n * math.pi / self.Omega_D) < self.dt / 2:
                self.ps_omega.append(self.omega[-1])
                self.ps_theta.append(self.theta[-1])
                n += 1
#            #$\Omega_Dt=n\pi-\pi/2$
#            if  abs(self.t[i + 1] - (n * math.pi - 0.5 * math.pi) / self.Omega_D) < self.dt / 2:
#                self.ps_omega.append(self.omega[-1])
#                self.ps_theta.append(self.theta[-1])
#                n += 1
#            #$\Omega_Dt=2n\pi+\pi/4$
#            if  abs(self.t[i + 1] - (2 * n * math.pi + math.pi / 4) / self.Omega_D) < self.dt / 2:
#                self.ps_omega.append(self.omega[-1])
#                self.ps_theta.append(self.theta[-1])
#                n += 1
            i += 1
            if self.t[i] > self.total_time:
                loop_calculate = False
    def show_results(self):
        pl.plot(self.ps_theta, self.ps_omega, '.')
        pl.xlabel('$\\theta$(radians)', fontsize=20)
        pl.ylabel('$\\omega$(radians/s)', fontsize=20)
        pl.title('$\\omega$ versus $\\theta$ $F_D=1.2$', fontsize=20)
        pl.show()
#按需求运行程序
start_1 = pendulums_basis(0, 60)
start_1.calculate()
start_1.show_results()
start_2 = pendulums_basis(0.5, 60)
start_2.calculate()
start_2.show_results()
start_3 = pendulums_basis(1.2, 60)
start_3.calculate()
start_3.show_results()
#start_comparison_1 = pendulums_comparison(0.5, 50)
#start_comparison_1.compare()
#start_comparison_2 = pendulums_comparison(1.2, 50)
#start_comparison_2.compare()
#start_ovt_1 = omega_versus_theta(0.5, 40)
#start_ovt_1.calculate()
#start_ovt_1.show_results()
#start_ovt_2 = omega_versus_theta(1.2, 200)
#start_ovt_2.calculate()
#start_ovt_2.show_results()
#start_ps = poincare_section(1.2, 7000)
#start_ps.calculate()
#start_ps.show_results()

3.12:

$\omega$ versus $\theta$,$F_{D}$=1.2,$\Omega _{D}t$=2n$\pi$
q=$\frac{1}{2}$,l=g=9.8,$\Omega _{D}=\frac{2}{3}$,$d_{t}=0.04$,$\theta (0)=0.2$,$\omega (0)=0$

$\omega$ versus $\theta$,$F_{D}$=1.2,$\Omega _{D}=n\pi -\frac{\pi }{2}$
q=$\frac{1}{2}$,l=g=9.8,$\Omega _{D}=\frac{2}{3}$,$d_{t}=0.04$,$\theta (0)=0.2$,$\omega (0)=0$

$\omega$ versus $\theta$,$F_{D}$=1.2,$\Omega _{D}=2n\pi +\frac{\pi }{4}$
q=$\frac{1}{2}$,l=g=9.8,$\Omega _{D}=\frac{2}{3}$,$d_{t}=0.04$,$\theta (0)=0.2$,$\omega (0)=0$

3.13:

$\Delta \theta$ versus t,$F_{D}$=0.5

$$q=1/2, l=g=9.8, \Omega_D=2/3, dt=0.04, \theta(0)_1=0.2, \theta(0)_2=0.201, \omega(0)=0$$

$\Delta \theta$ versus t,$F_{D}$=1.2

$$q=1/2, l=g=9.8, \Omega_D=2/3, dt=0.04, \theta(0)_1=0.2, \theta(0)_2=0.201, \omega(0)=0$$
此处本应有图，显不出来┌( ಠ_ಠ)┘宝宝好方

3.14

$\Delta \theta$ versus t,$F_{D}$=0.5
q_1=0.5, q_2=0.501, l=g=9.8, \Omega_D=2/3, dt=0.04, \theta(0)=0.2, \omega(0)=0

此处本应有图，显不出来(；°○° ) 我努力了

$\Delta \theta$ versus t,$F_{D}$=1.2
q_1=0.5, q_2=0.501, l=g=9.8, \Omega_D=2/3, dt=0.04, \theta(0)=0.2, \omega(0)=0

此处本应有图，显不出来(ง •̀_•́)ง 还是要加强知识水平。唉！

4.结论

3.12：
得到的图像与预期的结果相符。
3.13：
图像可知在$F_{D}=0.5$时系统运动是有规则的且收敛，曲线的渐近线很好的符合指数规律，指数$\lambda$<0. 但是当$F_{D}=1.2$时情况发生了改变，一开始还是稳定的但是之后就变得无规则不可预测，进入了混沌状态,曲线渐近线近似符合指数规律，指数$\lambda$>0。
3.14：
由图可知，单摆的差异为$q_1=0.5, q_2=0.501$,当$F_D=0.5$时，系统最终稳定在一个常数系统的Lyapunov指数不再存在，当$F_D=0.5$时曲线近似符合指数规律，指数$\lambda$>0,系统仍然是混乱的。

5.致谢

感谢张盛同学的细心教导还有张然带我了解liunx，还有蔡浩老师的指导。

---