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

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

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

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

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

如图，摆长为L的单摆从一较大起摆角θ

开始摆动，某时刻运动到摆角为θ

处。忽略摩擦力和空气阻力的作用，其运动遵循机械能守恒定律。取运动最低点为势能零点，则有：

⑴其中

代入⑴式整理得

两侧积分

t =

由单摆振动的对称性知

4T =⑵

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

C

限制θ≠θ0，又θ=θ

处存在一条渐近线，则⑵式为反常积分，无法用数值模拟求解。 考虑换元积分法，令

则⑵式化为

4T =*

其中

k Si =

*式为第一类完全椭圆积分，可将其展开成无穷级数形式

⑶

二、插值法

1、线性内插法 *式中令（π上为单值函数）。

取点

,

π，其中，得相应直线解析式：

则

其中

即

1T =⑷

其中

2T =

其相对误差表示为

2、抛物线内插法 *式中令

f(,k ϕ（π上为单值函数）。

取点

,

(,π,

(,π，得相应抛物线解析式：

8

则

其中

即

⑸

其中

其相对误差表示为

（T 由⑶式得出）

3、四点曲线内插法

*式中令

f(,k ϕ（

2π上为单值函数）。

取点

,

(,4π,

(,6π,

(,2π，得相应三次曲线解析式：

2y =

则

34T =

其中

即

⑹

其中

其相对误差表示为

（T 由⑶式得出）

4、猜想

通过比较T 1,T 2,T 3 相对误差的大小，猜想T 3比T 2精确，T 2比T 1精确。

二、Jacobi 椭圆函数法

具体方法：为了便于理解，先给出以下Jacobi 椭圆函数的定义，列出本文用到的几个基本公式与微分方程。

Jacobi 椭圆函数：一般的，正弦函数也可以用它的反函数来定义。如：

，

它的反函数可以定义为

0x

u =⎰

式中积分值u 是积分上限x 的函数。

类似地，Jacobi 椭圆正弦函数，它是第一类椭圆积分的反函数，第一类椭圆积分的定义是

F(,k ϕ

令z=sin θ，x=sin ψ，则

F(,z k （7）

我们定义：

(,u x k

（8）

的反函数为Jacobi 椭圆正弦函数sn （u ，k ），简记为snu 。注：k 为Jacobi 椭圆函数的模，p 为x 的振幅函数，记做ψ=amu ，当k=0时，退化为三角函数。

把Jacobi 椭圆正弦函数、Jacobi 椭圆余弦函数和第三类Jacobi 椭圆函数分别定义为：

sn u =

由（7）及（8）可知，当振幅ψ=π/2,即x=1时，u=K （k ）；这里

由定义式可得

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

就可以求得它们的导数公式如下：

由定义式（8）可知，当k=0时，snu 退化为sinu ，则k=0时，cnu 退化为cosu ，k=0时，dnu 退化为1 另一方面，当k=1时，由定义式可知

(u x ，

因而x=tanh u ，可见k=1时，sn 退化为tanh u ，

所以，k=1时，cnu 退化为sechu ，dnu 也退化为sechu ，说明有 周期性的椭圆函数当k=1时退化为没有周期性的双曲线函数。 最后由x=snu ；cnu 和dnu 的导数公式可知，它们分别满足微分方程：

2

2'x =（9）

单摆运动方程可由角动量定理写成

2d ml （10）

即

∙∙

方程（k ）的首次积分（能量积分）为

21 （11）

其中总能量H 取决于初始条件，且H>0,令参数

则方程可改写为

∙（12）

由（6）式知：若H<2mgl ，即k 2<1，单摆将作往复运动，摆幅θ0

（sin θ0/2=k ）随H 值的增大而增大，但若H>2mgl ，即k 2>1，则θ2>0,单摆将作绕悬挂点作旋转运动。 1） 当k ﹤﹤1时，θ0很小，单摆作谐振动。（

）

2）

当1> k 2>0时，单摆的振幅θ0随k 值的增大而增大。

3） 当k 2>1时，θ02>0，单摆作旋转运动。

本文将只讨论（2）情况，令x=sin θ/2,可把方程（12）写为 x 2=ω02（1- x 2）（k 2- x 2） （13） 1>k>0（大摆幅振动） 由snu 的导数公式

d sn

可立刻看出，尝试解

x=x 0sn （ω0t ，k ）能满足方程（13），代入方程（7）得：

0[1

比较左、右两边常数项及sn 2（ω0t ，k ）项、sn 4（ω0t ，k ）项的系数，可知应取x 0=k ，即k=sin θ0 /2

其中，θ0是单摆的振幅，故解为θ

(14)

其中，

s =

摆动的周期T 可以用椭圆正弦函数的周期4K （k ）求出，由 ω0T=4K （k ）即得

4（15）

它与单摆振动的周期T 0=2π/ω0之比是

（16）

它随k 的增大而增大，表1列的是θ0 =30°、60°、90°时，k ，K 和T/ T 0值 由表1得，若振幅θ0 <30°，则T 与T 0偏差约在2％以内，这样就能计算比5°大的摆角问题了，且数值较精确。

三、微分方程的数值解法

对于大幅角单摆的研究，由于其运动微分方程为非线性方程，这给计算带来

了很大的困难。通常我们可以运用例如Euler 法或经典4阶Runge-kutta 的方法