第一篇：Bazinga：大角度单摆的运动近似(1)：小角度近似

第二篇：Bazinga：大角度单摆的运动近似(2):周期的解与椭圆积分

上次我们已经成功解决了周期的表达式以及周期近似，接下来我们把注意力集中到微分方程 [ddot heta = - frac{g}{L}sin heta ] 的解 heta(t) 上来.此时我们不使用小角度近似，并试图求该微分方程的数值解.在往后的几篇文章中我会把自己的级数解法更新上来.在本文会先用一种我临时受到启发的解法来解决该问题.灵感出自《费曼物理讲义》.

首先对该物理问题所需要的数学基础做一个铺垫.

图片来自www.bing.com

本文的核心数学思想是线性近似.对于一个已知函数 y=f(x) ,由自变数在 x=a 处增加引起的增量 [Delta y = f(a + Delta x)] ,此增量可以表示为

[Delta y = f(a)Delta x + o(Delta x)]

在 [Delta x o 0] 时高阶无穷小量 [o(Delta x) o 0] ，这时可以定义微分了.不过这里我们采用近似的记法，即当 [Delta x] 很小时， [Delta y approx f(a)Delta x]

这个近似模型可以很好的运用在运动学中，比如在本文的情况下 [omega = frac{{d heta }}{{dt}}] , [alpha = frac{{domega }}{{dt}}]

所以可以近似的说，当我们的步长 [Delta t] 足够小的时候，

[ heta ({t_n}) = heta ({t_{n - 1}}) + omega ({t_n})Delta t]

[omega ({t_n}) = omega ({t_{n - 1}}) + alpha ({t_{n - 1}})Delta t]

注意在这里我们定义 [{t_n} = nDelta t] ,这说明了对于任意时刻的角度位移，只要知道它上一刻的角度位移和此刻的角速度就可以得到；而此刻的角速度又由上一刻的角速度和上一刻的角加速度给出.此时牛顿第二定律的微分方程给了我们便利，因为它成功的把角加速度和角度位移联系了起来：

[alpha ({t_n}) = frac{{{d^2}}}{{d{t^2}}} heta ({t_n}) = - frac{g}{L}sin heta ({t_n})]

这也就是说

[left{ egin{array}{l} heta ({t_n}) = heta ({t_{n - 1}}) + omega ({t_n})Delta t\ omega ({t_n}) = omega ({t_{n - 1}}) - frac{g}{L}sin heta ({t_{n - 1}})Delta t end{array} ight.]

故只要我们知道了前一时刻的角加速度和角度位移，就可以表示出这一时刻的角度位移.此时我们为了化简计算，令常数 [frac{g}{L} = 1] .同时为了获得较高的精确度，我们让 [Delta t = 0.01] .

我们还需要给出初始条件.考虑到我们研究的是单摆在大角度情况下的运动，那么我们令 t=0 时 [ heta ({t_0}) = frac{pi }{2}] ,而角速度自然是 [omega ({t_0}) = 0] .

为了画出角度位移随时间变化的图像，我让我同学 @天琴心弦 （编程大佬）帮我写了一段Python代码，从而帮助我画出这个函数:

import numpy as np
import matplotlib.pyplot as plt

telement = list()
telement.append(np.pi/2)
for i in range(1000):
telement.append(0)
oelement = list()
oelement.append(0)
for i in range(1000):
oelement.append(0)

def theta(n):
if n == 0:
return telement[n]
else:
print("This is "+ str(n))
if telement[n] == 0:
telement[n] = telement[n - 1] + 0.01 * omega(n)

return telement[n]

def omega(n):
if n == 0:
return oelement[n]
else:
if oelement[n] == 0:
oelement[n] = oelement[n - 1] - np.sin(telement[n - 1]) * 0.01

return oelement[n]

x = list()
for i in range(1000):
x.append(i*0.01)
y = list()
for i in range(1000):
y.append(theta(i))

plt.plot(x, y)
plt.show()

运行之后大概是这个样子：

我让时间t以0.01的步长运行了1000次，所以得出来的应该是在前10s内角度的变化，图像如下：

在前10s内角度的变化

看上去，该函数是一个三角函数，所以当摆角很小的时候，单摆运动可以近似的看作简谐运动.那么当摆角逐渐变大时，会发生什么情况？

释放角度为3pi/4(135度)
释放角度为0.99pi(178.2度)
释放角度为pi(180度)

当释放角度刚好为 pi 时程序会出现问题，这是因为 sin(pi)=0 ，这就导致角速度一直不会有变化，故角度就不会有变化.所以在最后的图象上角度总是一个常量 pi .

从上面几张图片可以找到一些规律：

1.随著释放角度的增大周期也会跟著增大

2.角度的变化始终是周期性的

3.当释放角度越来越大时，简谐运动的近似会越来越不精确.

通过上面的程序我们已经初步知道了如何用数值表示这个函数.下次我会讲讲如何求这个函数的解析解.

参考

[1] 《费恩曼物理学讲义》 第一卷 第9章 牛顿的运动学方程


推荐阅读：
相关文章