高云峰

(清华大学航天航空学院,北京100084)

一般来说，在理论力学的动力学教学中，绝大部分问题都只需要学生会受力分析、列写动力学方程，并不要求进一步求出具体的运动[1]。原因之一是绝大部分问题没有解析解，用传统数学方法不好解；原因之二是学时限制，没有更多时间教学生如何求解。

随着计算软件的发展，目前已经可以很轻松地求出动力学方程的解，包括运动和力随时间的变化关系，也可以很容易在屏幕上进行动画演示。

本篇先从动力学中最简单的单摆问题开始介绍：单摆的动力学方程如何求解？数值解的精度如何？大角度情况下周期如何变化？然后介绍椭圆摆的相关问题。

1 单摆的相关问题[2]

案例1：假设单摆中小球A质量为m，尺寸不计，绳子OA长度为l，不计质量(图1)。试比较不同初始角度对运动的影响，以及如何证明计算的结果可靠？

图1 单摆模型

单摆运动时，根据其受力图(图2)，可以在切向方向列写系统的动力学方程，为

图2 单摆的受力图

虽然方程(1)在线性化时有解析解，但是大角度时比较复杂。但在Matlab中，可以直接调用ode45函数统一求解这类常微分方程，其格式是[3]

其中黑体是固定的格式，斜体是变量或自己命名的函数。其中options是积分的选项，可以控制积分的精度，′RelTol′表示积分中的相对误差，′AbsTol′是积分中的绝对误差。在Matlab数值积分计算中，通常number1和number2可以选为10−8∼10−10(精度太高会花费更多计算时间，也没有必要),如果省略则默认为10−3。ode45函数则是求解常微分方程的一种常用函数，其中[t,y]是积分的时间和结果；动力学方程在子程序filename中描述；y0是初始条件；[starttime:steptime:endtime]表示积分时按等步长steptime从开始积分到结束(等步长积分是为了动画演示时每一帧时长相同)。

注意ode45求解标准的一阶常微分方程组，因此要把动力学中的二阶微分方程转换为一阶微分方程组。把方程(1)这样处理：设y1=θ,y2=，得到一阶微分方程组

初始条件为y1(0)=θ(0)，y2(0)=˙θ(0)。

单摆动力学方程求解的源代码见图3，子程序源代码见图4。

图3 单摆动力学求解的源代码

图4 动力学子程序源代码

主程序源代码中global是表示全局变量，在主程序中赋值后在其他子程序中可以直接调用；plot是画线段的函数，如果很密集，各段微小的线段就构成了曲线；y(:,1)表示y数组中的第1列，实际上就是每一步计算得到的单摆角度；xlabel和ylabel是给图中坐标轴加上标注，科学论文中一般需要知道坐标是什么含义，什么单位。

子程序用function开头，注意在Matlab中子程序文件名与函数名相同；zeros(2,1)表示生成一个2×1的列阵，里面元素都是0；子程序中的动力学方程直接按照式(3)写出。

数值计算中初始角度可以变化，其他参数不变，如下所示

图5是方程(1)在不同初始角度下的解(暂时没有考虑空气阻力)，可以看到解的周期与初始参数有关，这也是非线性方程的特点。具体计算时可以把不同的初值积分结果赋值给y1，y2，y3，y4，在画图plot命令后使用hold on命令可以把不同的曲线叠加在一起；title是给图命名；legend相当于图例，可以给不同的曲线命名以示区别，见图6。

图5 不同角度单摆的解

图6 多个曲线画在同一图上

2 数值计算的可靠性

数值计算通常是有误差的，包括数值的截断误差、算法本身的误差等。不过随着计算技术的发展，可显示的最小值越来越小，例如在Matlab中输入realmin(′single′)，显示出单精度最小正浮点数为1.175 494 4×10−38，所以普通计算精度完全满足要求。

那么方程(1)的数值积分精度如何呢？这首先取决于options中的精度控制，在前面所说的10−8∼10−10情况下，积分计算精度也基本是这一数量级。怎么证明这一点呢？可以通过方程的守恒量来判断。

不考虑空气阻力时，方程(1)有守恒量(机械能)

守恒量理论上是一个常数，因此可以用于检测数值计算的可靠性。在数值计算中，可以先把每一步的角度、角速度等量计算出来，然后计算每一步的E。考虑到计算有误差，这个“守恒量”应该在极小范围内变化才合理。图7是不同初始条件下方程的积分结果，看起来很平坦，从而证明积分的结果很可靠。

图7 不同条件下的积分常数

如果想了解积分常数变化的细节，利用能量之差是有效的方法。能量之差定义为

如果没有事先指定等比例(用axis equal命令表示x和y轴等比例)，Matlab在画图时坐标轴会自动根据数据范围进行缩放，因此能量之差的细节变化就可以反映出来(图8)，根据能量之差，可以看出系统机械能的变化范围与积分的允许误差范围同阶，初始角度小时计算误差更小。

图8 能量之差

上述分析表明数值计算是可靠的，然后再考虑系统有阻尼的情况。图9显示了不同阻尼情况下摆角的变化情况，初始摆角均是30°，其中β=n/(2m)是阻尼系数。图10显示了阻尼对系统机械能的影响，最后能量趋于平坦，为系统静止时的机械能。图10中一个细节是：阻尼系数较小时，机械能呈现台阶状的下降，这是因为摆角每次到最大幅值时速度接近，此时空气阻力很小，机械能接近守恒，所以机械能在下降过程中就形成了台阶。

图9 阻尼对摆角的影响

图10 阻尼对机械能的影响

3 椭圆摆的相关问题

椭圆摆是动力学中的一个典型问题，我们想了解一下：椭圆摆的运动是周期的吗？其周期是多少？

案例2：椭圆摆(图11)由质量为mA的滑块A和质量为mB的单摆B构成。滑块可沿光滑水平面滑动，AB杆长为l，质量不计。试建立系统的运动微分方程，并求水平面对滑块A的约束力。

图11 椭圆摆模型

系统有2个自由度，建立Oxy坐标系，取A点位移x和AB杆的相对转角φ为广义坐标(图12)，传统分析可得到系统的运动微分方程以及压力的表达式(过程略)为

通常理论力学教材或教学到了这一步就算结束了。椭圆摆在运动过程中A滑块和AB杆是否为周期运动？是否为正弦运动？压力是如何变化的？类似这样的问题传统方法都不好回答。

系统的运动是周期的吗？在式(7)中消去¨x，有

可以看出，摆角的方程(9)在小角度、线性化后是周期的；类比单摆，大角度情况下摆动周期与初始条件有关。下面的计算只改变椭圆摆的初始角度φ0，其他参数均不变。

图13显示了初始条件对摆角的影响，看起来不同条件下摆角都是周期函数，但是周期不同。图14进一步显示了周期和初始条件的关系，同时表明：初始角为小量时与线性化的结果比较接近。

图13 椭圆摆的摆角曲线

图14 椭圆摆的摆动周期

图15和图16表明位移和压力也是周期函数，但是初始大角度时位移和压力曲线已经明显偏离标准正弦曲线了。如果没有数值计算，这些细节无法得到。

图15 椭圆摆的位移曲线

图16 椭圆摆的压力曲线

另外还有一个现象，系统的位移、摆角、压力虽然都是周期函数，但是周期并不相同。虽然椭圆摆的摆角看上去像余弦曲线(图13)，但是否为标准的余弦曲线？这可以使用快速傅里叶变换来分析。在Matlab中，可以直接调用fft(dataname)函数求出数据dataname的频率。

在处理椭圆摆之前，先看看fft的效果，设测试函数为

其中A0=0.5,A1=1,A2=0.4,A3=0.2;f1=3,f2=5,f3=10.6，其时域图(图17)看不出什么规律，但经过快速傅里叶变换变化后，在频域图中就可以清楚看出只会在f=fi处有一个孤立的峰值Ai，而在其他位置均为0：因此在频域图中可以直接得出测试函数(11)中的所有参数(图18)。

图17 测试函数的时域曲线

图18 测试函数的频域曲线

下面对椭圆摆的摆角和位移进行快速傅里叶变换变化(位移通过平移去掉直流分量)，为了方便比较不同初始条件的结果，把快速傅里叶变换变化后的结果归一化处理：最大值取为1，处理后的结果见图19和图20，可以发现：

图19 角度曲线的频谱

图20 位移曲线的频谱

(1)可以看到各条曲线都有一个明显的峰值，但其他位置还有连续的不全为0的值，且大角度时高倍频处还会有小峰值，这说明椭圆摆的摆角不再是标准的余弦函数；

(2)由于峰值相对明显，所以在时域图中看起来很像余弦曲线；

(3)各曲线最大峰值对应的主频率不同：频率随初始角度增加而减少，或周期随初始角度增加而增加，符合图14的结论，这反映了非线性函数的周期或频率与初始条件有关。

4 小结

现代科学的发展，使得计算的重要性大为提升，它和理论分析、试验验证同为科学研究的三大重要手段。例如，1963年MIT的气象学家Lorenz在计算大气对流问题时发现了混沌现象、现代飞机的设计涉及大量动力学软件的计算，等等。

在这种背景下，在动力学中引入Matlab，除了会列写动力学方程，还能计算和演示系统整个运动过程，便于发现系统丰富的动力学现象。

本篇着重介绍了动力学中运动微分方程的求解，首先把二阶微分方程转换为一阶微分方程组，然后采用荣格库塔法进行求解。涉及到的Matlab核心函数是odeset(设置积分精度)和ode45(求微分方程)；plot和hold on命令可以实现多个图叠在一起。另外介绍了fft函数(把时域数据转换到频域)，可以分析复杂曲线的频率或周期。

可以看出，借助Matlab可以更加深入地研究动力学问题，解决传统分析方法无法处理的问题。