大幅角单摆振动周期的研究

大幅角单摆振动周期研究的综述


周越 赵诚嘉 刘汉发 王丽丽 张文贤 张桢 阚健 等（B机制091）

摘要：运用插值法、Jaco bi椭圆函数
法、
微分方程的数值解法，研究了大幅度单
摆的周期公式。由插值法计算 出的单摆周期与精确周期相比其相对误差较小，与
其他文献得到的周期公式相对误差相比是最小的，具有 较强的实用性。采用
Jacobi椭圆函数分析单摆的运动，结果的误差也是较小的。运用经典4阶Runge-kutta的方法来研究周期公式，并结合Matlab作图。

一、任意摆角下单摆周期公式的推导

如图，摆长为L的单摆从一较大起摆角θ0
开始摆动，某时刻运动到摆角为θ
处。忽略摩擦力和空气阻力的作用，其运动遵循机械能 守恒定律。取运动最低点
为势能零点，则有：

其中
代入⑴式整理得
⑴



两侧积分

由单摆振动的对称性知


⑵式中θ的积分上界为θ
0
，但被积表达式
⑵

限制θ≠θ
0
，又θ=θ
0
处存在一条渐近线，则⑵式为反常积分，无法用数值模拟求解 。

考虑换元积分法，令
则⑵式化为


其中
*

*式为第一类完全椭圆积分，可将其展开成无穷级数形式
二、插值法
1、线性内插法
*式中令
取点
,

则

其中
即
t?
T?4
?
k
T
?Si
4
?
?
C

⑶

（
上为单值函数
）。
，其中，
得相应直线解析式：

其中

其相对误差表示为

2、抛物线内插法
*式中令
取点
,
则
其中
即

T
其中
其相对误差表示为
f(

3、四点曲线内插法
y

(

T
（
上为单值函数
）。
，得相应抛物线解析式：
（T由< br>(
?
⑶式得出
?
1
?
?
）
?

⑷



k

⑸
,

?
?
2
8
,
,
,

*式中令
（
,,
上为单值函数
）。
取点
则
其中
即

其中
其相对误差表示为
4、猜想
通过比较T
1
,T
2
,T
3
相对误差的大小，猜 想T
3
比T
2
精确，T
2
比T
1
精确。

二、Jacobi椭圆函数法
具体方法：为了便于理解，先给出以下Jacobi 椭圆函数的定义，列出本文
用到的几个基本公式与微分方程。
Jacobi椭圆函数：一般的 ，正弦函数也可以用它的反函数来定义。如：
?
?
?
?
2
f (
?
,k
((,,
(,
y?
2
62
T?4
4
3
k
T
?
?
Si
?
T
,
，得相应三次曲线解析式：




⑹


（
T由⑶式得出
）
，

它的反函数可以定义为

式中积分值u是积分上限x的函数。
类似地，Jacobi椭圆正弦函数，它是第一类椭圆积分的反函数，第一类椭圆积分
的定义是

令z=sinθ，x=sinψ，则
（7）
我们定义：
的反 函数为Jacobi椭圆正弦函数sn（u，k），简记为snu。注：k为Jacobi椭圆
函数的模 ，p为x的振幅函数，记做ψ=amu，当k=0时，退化为三角函数。
把Jacobi椭圆正弦函数 、Jacobi椭圆余弦函数和第三类Jacobi椭圆函数分别定
义为：
由（7）及（8） 可知，当振幅ψ=π2,即x=1时，u=K（k）；这里
由定义式可得
就可以求得它们的导数 公式如下：
u?
?
F(
?
,k
0
F(z,ksnu?
u
d
(
?
x,k
（8）


，于是利用此式和sn，cn，dn三者的定义式，

x

< br>由定义式（8）可知，当k=0时，snu退化为sinu，则k=0时，cnu
退化为cosu，k=0时，dnu退化为1
另一方面，当k=1时，由定义式可知
因而x=tanh u，可见k=1时，sn退化为tanh u，
所以，k=1时，cnu退化为sechu，dnu也退化为sechu，说明有
周期性的椭圆函数当k=1时退化为没有周期性的双曲线函数。
最后由x=snu；cnu和dnu的导数公式可知，它们分别满足微分方程：
，
（9）
单摆运动方程可由角动量定理写成
（10）
即
方程（k）的首次积分（能量积分）为
其中总能量H取决于初始条件，且H>0,令参数
则方程可改写为
由（6）式知：若 H<2mgl，即k
2
<1，单摆将作往复运动，摆幅θ
0

（si nθ
0
2=k）随H值的增大而增大，但若H>2mgl，即k
2
>1，则θ
2
>0,单摆将作
绕悬挂点作旋转运动。
1）

2）

当k﹤﹤1时，θ
0
很小，单摆作谐振动。（
当1> k
2
>0时，单摆的振幅θ
0
随k值的增大而增大。
u(x
x'?
??
2
d
1
ml
2
?
2
?
?
?

（11）

（12）

）
2

3）

当k
2
>1时，θ
0
2
>0，单摆作旋转运动。
本文将只讨论（2）情况，令x=sinθ2,可把方程（12）写为
x
2
=ω
0
2
（1- x
2
）（k
2
- x
2
） （13）
1>k>0（大摆幅振动）
由snu的导数公式
x=x
0
sn（ω
0
t，k）能满足方程（13），代入方程（7）得：

可立刻看出，尝试解
比较左、右两边常数项及sn
2
（ω
0
t，k）项、sn
4
（ω
0
t，k）项的系数，可知应
取x
0
=k，即k=sinθ
0
2
其中，θ
0
是单摆的振幅，故解为
其中，
摆动的周期T可以用椭圆正弦函数的周期4K（k）求出，由
ω
0
T=4K（k）即得
它与单摆振动的周期T
0
=2π ω
0
之比是
它随k的增大而增大，表1列的是θ
0
=30°、60°、90°时，k，K和T T
0
值
θ
0
k= sinθ
0
2 K TT
0
=2Kπ
30° （√3-1）2√2≈0.259 1.598 1.02
60° 12=o.500 1.686 1.07
1.854 1.18 90° 1√2≈0.707
由表1得，若振幅θ
0
<30°，则T与T
0
偏差约在2％以内， 这样就能计算比5°
大的摆角问题了，且数值较精确。
三、微分方程的数值解法
对 于大幅角单摆的研究，由于其运动微分方程为非线性方程，这给计算带来
了很大的困难。通常我们可以运 用例如Euler法或经典4阶Runge-kutta的方法
d
x
0
?0
sn
[1
?
4
sin
T
k?s
du
T?
(14)
（15）
（16）

来研究它。
具体方法：我们知道单摆的运动微分方程为
（17）
并可推导出单摆运动周期的精确公式为：

对于这样的积分是很难求得精确解析解。
我们可以运用经典四阶Runge- kutta法来求解。
将式（17）降阶为

< br>利用Matlab下求解一阶微分方程组初值问题数值解的最常用方法中的
ode45（）函数， 该函数采用的是变步长的四阶五阶Runge-kutta- Felhberg算法。
并由经典的Runge-kutta方法理论分析可知，四阶Runge-ku tta方法其收敛阶为
O(h^5)(h为计算过程中所采用的时间步长)，计算选取的适当的小就可以 保证所
需精度的解
具体的程序省略不写，我们可以得出最大摆角与周期的关系。得出的 结果与
式（3）的前12项相比较。当单摆角度较小时，周期随摆角变化缓慢，基本成线
性关系 ；在摆角较大时，周期变化很快。
四、结语
利用线性内插法对第一类完全椭圆积分进行近似处理得到的大幅度单摆振
动周期的公式, 与其 他文献相比本文相对误差最小、精度最高,而且公式简单。
力学中的一些问题，我们希望能够求得精确的 解析解，但往往很难做到，这是我
们就可以运用Matlab求得比较精确的数值解，并对运动的物体有 一个比较直观
的认识。
在史友进教授的帮助下，本论文得以顺利完成，在此我们对史老师表示衷心
的感谢！！

d
?
dy
2
4T?
?
dt


⑵