作业五Oscillation(matplotlib)

作业 单摆 Euler-Cromer方法

作者：许晗 2013301020084 武汉大学物理学院

程序目的

• 研究单摆在不同近似条件和外界条件下的运动。
• 每次计算前可以选择用哪些近似。（只是图上标记就有点难受了...）
• 用matpoltlib画图。
• 学习Blender建模。(还没有)
• 最后用Vpython展示混沌现象。(还没有)

实现方法与计算

1. 从最简单的单摆开始，设单摆长为 $l$，以竖直方向的夹角 $\theta$为广义坐标的运动方程:
   

   $$\frac{d^2\theta}{dt^2}=-\frac{g}{l}sin\theta \tag 1$$
   将一个二次方程变成两个一次方程，用到$\omega = \frac{d\theta}{dt}$

   Euler-Cromer方法:
   

   $$\omega _{i+1}=w_i-\frac{g}{l}\theta_i \Delta t \tag 2$$
   $$\theta_{i+1}=\theta_i+\omega_{i+1}\Delta t$$
   $$t_{i+1}=t_i+\Delta t$$
   重复足够次数。

   • 这种情况下有解析解：
     $$\theta = \theta_0 sin(\Omega t + \phi)$$
     其中 $\Omega = \sqrt{g/l}$。

2. 不那么平庸：考虑上受到的摩擦，经验告诉我们，在多数情况下可以认为大小正比于速度$l(d\theta / dt)$，q为常数
   

   $$F_{drag}=-q\frac{d\theta}{dt} \tag 3$$
   只需在Euler-Cromer计算过程中(2)式计算$\omega$加上
   $$-q\omega_i \Delta t \tag{3.1}$$
   • 阻力很小(underdamped)有近似解:
     $$\theta (t)=\theta_0 e^{-\frac{qt}{2}}sin(\sqrt{\Omega^2-q^2/4}t+\phi)$$
   • 阻力很大(overdamped)有近似解:
     $$\theta (t)=\theta_0 e^{-(\frac{qt}{2}\pm\sqrt{\Omega^2-q^2/4}t)}$$
   • 阻力在临界值(critically damped)有近似解：
     $$\theta (t)=(\theta_0 + Ct)e^{-qt/2}$$

3. 再特殊一点：外界提供能量，即加外力。而一个周期变化的力是常见的，取正弦变化：
   

   $$F_p/m = F_D sin(\Omega_D t) \tag 4$$
   只需在Euler-Cromer计算过程中(2)式再加一项
   $$+F_D sin(\Omega_D t)\Delta t \tag {4.1}$$
   • 但是这仍有解析解：
     $$\theta (t) = \theta_0sin(\Omega_Dt+\phi)$$
     其中振幅$\theta_0$
     $$\theta_0 = \frac{F_D}{\sqrt{(\Omega^2-\Omega_{D}^2)^2+(q\Omega_D)^2}}$$
     用Euler-Cromer方法数值计算的结果(课本3.7题), $\Omega = \sqrt{9.8}$

   • 比例系数$q = 0$，外力振幅$F_D=0.2$。
     $\Omega_D$ = [1.0, 2.0, 3.0]
     

   • 外力振幅$F_D=0.2$, 外力频率$\Omega_D = 1.0$。
     比例系数$q$ = [q = [0, 1.0, 2.0]]
     

致谢

参考了岳绍圣同学的作图。

链接：
混沌现象（单摆）