《理论力学》课程教学资源（文献资料）Matlab求解理论力学问..._单摆和椭圆摆的运动及周期_高云峰.pdf_大学文库

第43卷第4期力学与实践2021年8月Matlab求解理论力学问题系列（三)单摆和椭圆摆的运动及周期高云峰1)（清华大学航天航空学院，北京100084)一般来说，在理论力学的动力学教学中，绝大部分问题都只需要学生会受力分析、列写动力学方程，并不要求进一步求出具体的运动[1]。原因之一是绝10大部分问题没有解析解，用传统数学方法不好解：原10因之二是学时限制，没有更多时间教学生如何求解。随着计算软件的发展，目前已经可以很轻松地mid求出动力学方程的解，包括运动和力随时间的变化关系，也可以很容易在屏幕上进行动画演示。本篇先从动力学中最简单的单摆问题开始介绍：mg单摆的动力学方程如何求解？数值解的精度如何？大图2单摆的受力图角度情况下周期如何变化？然后介绍椭圆摆的相关问题。虽然方程（1)在线性化时有解析解，但是大角度时比较复杂。但在Matlab中，可以直接调用ode451单摆的相关问题[2]函数统一求解这类常微分方程，其格式是[3]案例1：假设单摆中小球A质量为m，尺寸不计，绳子OA长度为l，不计质量（图1)。试比较不options=odeset(RelTol',numberl,同初始角度对运动的影响，以及如何证明计算的结'AbsTol', number2);(2)果可靠？[t, ] =ode45('file-name', [start.time :单摆运动时，根据其受力图（图2)，可以在切向step-time:end_timel,y0, options);(3)方向列写系统的动力学方程，为其中黑体是固定的格式，斜体是变量或自己命名的mlo = -mg sino - nlo(1)函数。其中options是积分的选项，可以控制积分的精度，RelTol'表示积分中的相对误差，AbsTol'是积分中的绝对误差。在Matlab数值积分计算中，通常number1和number2可以选为10-8~10-10（精度太高会花费更多计算时间，也没有必要），如果省略则默认为10-3。ode45函数则是求解常微分方程的OA一种常用函数，其中[t，习是积分的时间和结果；动力学方程在子程序file_name中描述：yo是初始条图1单摆模型本文于2020-06-01收到。1)E-mail:gaoyunfeng@tsinghua.edu.cn引用格式：高云峰.Matlab求解理论力学问题系列（三）单摆和椭圆摆的运动及周期.力学与实践，2021，43(4)：593-598Gao Yunfeng. The movement and the period of simple pendulum and elliptical pendulum.Mechanics in Engineering.2021, 43(4): 593-598(C)1994-2023ChinaAcademic JournalElectronicPublishingHouse.Allrightsreserved.http://www.cnki.net

第 43 卷第 4 期力学与实践 2021 年 8 月 Matlab 求解理论力学问题系列 (三) 单摆和椭圆摆的运动及周期高云峰1) (清华大学航天航空学院, 北京 100084) 一般来说，在理论力学的动力学教学中，绝大部分问题都只需要学生会受力分析、列写动力学方程，并不要求进一步求出具体的运动 [1]。原因之一是绝大部分问题没有解析解，用传统数学方法不好解；原因之二是学时限制，没有更多时间教学生如何求解。随着计算软件的发展，目前已经可以很轻松地求出动力学方程的解，包括运动和力随时间的变化关系，也可以很容易在屏幕上进行动画演示。本篇先从动力学中最简单的单摆问题开始介绍：单摆的动力学方程如何求解？数值解的精度如何？大角度情况下周期如何变化？然后介绍椭圆摆的相关问题。 1 单摆的相关问题 [2] 案例 1：假设单摆中小球 A 质量为 m，尺寸不计，绳子 OA 长度为 l，不计质量 (图 1)。试比较不同初始角度对运动的影响，以及如何证明计算的结果可靠？单摆运动时，根据其受力图 (图 2)，可以在切向方向列写系统的动力学方程，为 ml¨θ = −mg sin θ − nl ˙θ (1) θ O A 图 1 单摆模型 T mg mlθ lθ 2 lθ A . . . θ 图 2 单摆的受力图虽然方程 (1) 在线性化时有解析解，但是大角度时比较复杂。但在 Matlab 中，可以直接调用 ode45 函数统一求解这类常微分方程，其格式是[3] options = odeset (′RelTol′ , number1, ′AbsT ol′ , number2); (2) [t, y] = ode45 (′ f ile name′ , [start time : step time : end time], y0, options); (3) 其中黑体是固定的格式，斜体是变量或自己命名的函数。其中options是积分的选项，可以控制积分的精度，′RelTol′ 表示积分中的相对误差，′AbsTol′ 是积分中的绝对误差。在 Matlab 数值积分计算中，通常 number1 和 number2 可以选为 10−8 ∼ 10−10 (精度太高会花费更多计算时间，也没有必要), 如果省略则默认为 10−3。ode45函数则是求解常微分方程的一种常用函数，其中 [t, y] 是积分的时间和结果；动力学方程在子程序 f ile name 中描述；y0 是初始条本文于 2020–06–01 收到。 1) E-mail: gaoyunfeng@tsinghua.edu.cn 引用格式: 高云峰. Matlab 求解理论力学问题系列 (三) 单摆和椭圆摆的运动及周期. 力学与实践, 2021, 43(4): 593-598 Gao Yunfeng. The movement and the period of simple pendulum and elliptical pendulum. Mechanics in Engineering, 2021, 43(4): 593-598

力学与实践5942021年第43卷件：[start_time：step-time：end_timel表示积分时按子程序用function开头，注意在Matlab中子等步长step.time从开始积分到结束（等步长积分是程序文件名与函数名相同：zeros(2.1)表示生成一个为了动画演示时每一顿时长相同。2×1的列阵，里面元素都是0：子程序中的动力学注意ode45求解标准的一阶常微分方程组，因此方程直接按照式（3）写出。数值计算中初始角度可以变化，其他参数不变，要把动力学中的二阶微分方程转换为一阶微分方程组。把方程(1)这样处理：设91=9,y2=，得到一如下所示阶微分方程组m=1kg, 1=1m, g=9.8 m/s, (0)=0 (5)91=92(4)图5是方程（1）在不同初始角度下的解（暂时y2=-gsinyi/l-ny2/m没有考虑空气阻力)，可以看到解的周期与初始参数初始条件为y1(0)=0(0)，92(0)=(0)。有关，这也是非线性方程的特点。具体计算时可以把单摆动力学方程求解的源代码见图3，子程序源不同的初值积分结果赋值给yl，y2，y3，y4，在画图代码见图4。plot命令后使用holdon命令可以把不同的曲线叠加在一起；title是给图命名；legend相当于图例，可55539555888856581585955555585第单摆的动力学计算转以给不同的曲线命名以示区别，见图6。55565555555956655555516555555555555355585cle:清楚屏幕clearall :清除内存100globalmssRlengthon50参质量、长度、证和阻尼系数55=1:1enath0=1t=9.8.n=030°alltine-10;hh-0.01;积分时间、步长60050初始条件切始条件90°y0(1)=30*pi/180y0(2)=0:积分工作()/斯options.=odeset(RelTol',le-8,"AbsTol',le-10);“积分的误差选项0[t,y]-ode45Crg_kt_pendulum，[o:hh:alltine],yo,options);$积分结巢放在iy中丽图plot(t,y(:,D)+180/pl，r)）：画曲线gridon;出现网格-50titlec角度变化示意理）图的名称xlabel（时间/s）轴的标注ylabel（角度/deg）轴的标注-100102345图3单摆动力学求解的源代码时间/s电659596669595996969696969695966596669699696669669图5不同角度单摆的解微分方程的动力学方程品5959品9%%9696969665品6696566%function ydot-rg_kt_pendulum(t.y)$子程序plot(tl,yl(:,1)*180/pi,r-):holdonglobal massglengthon$全局变量plot(t2,y2(:,1)*180/pi,r-):ho1donydot-zeros(2.1);给2+1列阵腻0plot(t3,y3(:,1)*180/pi,'r:");holdontheta"y(1):dtheta_dt=y(2):%让变量有物理含义plot(t4,y4(:,1)*180/pi,"r-."):holdon一阶形式的动力学方程ydot(1)=dtheta_dt:全动力学方程legend（5度，30度，60度，90度）%曲线名ydot(2)=-g*sin(theta)/lengtho-1*dtheta_dt/mass：动力学方程2%出现网格grid on:title（角度变化示意图）%图的名称图4动力学子程序源代码xlabel（时间/s）%x轴的标注ylabel（角度/deg）%y轴的标注主程序源代码中global是表示全局变量，在主图6多个曲线画在同一图上程序中赋值后在其他子程序中可以直接调用；plot是画线段的函数，如果很密集，各段微小的线段就构成数值计算的可靠性2了曲线；y(,1)表示y数组中的第1列，实际上就是每一步计算得到的单摆角度：xlabel和ylabel是给数值计算通常是有误差的，包括数值的截断误图中坐标轴加上标注，科学论文中一般需要知道坐差、算法本身的误差等。不过随着计算技术的发标是什么含义，什么单位。展，可显示的最小值越来越小，例如在Matlab中(C)1994-2023ChinaAcademicJournalElectronicPublishingHouse.All rightsreserved.http://www.cnki.net

594 力学与实践 2021 年第 43 卷件；[start time : step time : end time] 表示积分时按等步长 step time 从开始积分到结束 (等步长积分是为了动画演示时每一帧时长相同)。注意ode45求解标准的一阶常微分方程组，因此要把动力学中的二阶微分方程转换为一阶微分方程组。把方程 (1) 这样处理：设 y1 = θ, y2 = ˙θ，得到一阶微分方程组 y˙1 = y2 y˙2 = −g sin y1/l − ny2/m } (4) 初始条件为 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) 写出。数值计算中初始角度可以变化，其他参数不变，如下所示 m = 1 kg, l = 1 m, g = 9.8 m/s 2 , ˙θ(0) = 0 (5) 图 5 是方程 (1) 在不同初始角度下的解 (暂时没有考虑空气阻力)，可以看到解的周期与初始参数有关，这也是非线性方程的特点。具体计算时可以把不同的初值积分结果赋值给 y1，y2，y3，y4，在画图 plot 命令后使用 hold on 命令可以把不同的曲线叠加在一起；title 是给图命名；legend 相当于图例，可以给不同的曲线命名以示区别，见图 6。 0 1 2 3 4 5 ᰦ䰤/s -100 -50 50 100 䀂ᓖ/(°) 5° 30° 60° 90° 0 图 5 不同角度单摆的解图 6 多个曲线画在同一图上 2 数值计算的可靠性数值计算通常是有误差的，包括数值的截断误差、算法本身的误差等。不过随着计算技术的发展，可显示的最小值越来越小，例如在 Matlab 中

第4期高云峰：Matlab求解理论力学问题系列（三）单摆和椭圆摆的运动及周期595输入realmin('single')），显示出单精度最小正浮点数角的变化情况，初始摆角均是30°，其中β=n/2m）为1.1754944×10-38，所以普通计算精度完全满足是阻尼系数。图10显示了阻尼对系统机械能的影响，要求。最后能量趋于平坦，为系统静止时的机械能。图10那么方程（1）的数值积分精度如何呢？这首先中一个细节是：阻尼系数较小时，机械能呈现台阶状取决于options中的精度控制，在前面所说的10-8~的下降，这是因为摆角每次到最大幅值时速度接近，10-10情况下，积分计算精度也基本是这一数量级。X10-s怎么证明这一点呢？可以通过方程的守恒量来判断。不考虑空气阻力时，方程（1)有守恒量（机械能）E=ml?g?- mgl cos e(6)2二-2/之守恒量理论上是一个常数，因此可以用于检测数值计算的可靠性。在数值计算中，可以先把每一-5030°步的角度、角速度等量计算出来，然后计算每一步60°的E。考虑到计算有误差，这个“守恒量”应该在极90°小范围内变化才合理。图7是不同初始条件下方程234015的积分结果，看起来很平坦，从而证明积分的结果很时间/s可靠。图8能量之差230β=1/408=1/2205°-5130910..6090°(0)/-10-20-10-30012345246810120时间/s时间/s图7不同条件下的积分常数图9阻尼对摆角的影响8.0如果想了解积分常数变化的细节，利用能量之B=1/4差是有效的方法。能量之差定义为B=1/2B=5/2-8.5(7)△E=E(t)-E(O)00如果没有事先指定等比例（用axisequal命令表区示r和y轴等比例)，Matlab在画图时坐标轴会自o动根据数据范围进行缩放，因此能量之差的细节变-9.5化就可以反映出来（图8），根据能量之差，可以看出系统机械能的变化范围与积分的允许误差范围同阶，-10.0246810120初始角度小时计算误差更小。时间/s上述分析表明数值计算是可靠的，然后再考虑系统有阻尼的情况。图9显示了不同阻尼情况下摆图10阻尼对机械能的影响(C)1994-2023ChinaAcademic JournalElectronicPublishingHouse.All rightsreserved.http://www.cnki.net

第 4 期高云峰：Matlab 求解理论力学问题系列 (三) 单摆和椭圆摆的运动及周期 595 输入 realmin(′ single′ )，显示出单精度最小正浮点数为 1.175 494 4 × 10−38，所以普通计算精度完全满足要求。那么方程 (1) 的数值积分精度如何呢？这首先取决于options中的精度控制，在前面所说的 10−8 ∼ 10−10 情况下，积分计算精度也基本是这一数量级。怎么证明这一点呢？可以通过方程的守恒量来判断。不考虑空气阻力时，方程 (1) 有守恒量 (机械能) E = 1 2 ml2 ˙θ 2 − mgl cos θ (6) 守恒量理论上是一个常数，因此可以用于检测数值计算的可靠性。在数值计算中，可以先把每一步的角度、角速度等量计算出来，然后计算每一步的 E。考虑到计算有误差，这个 “守恒量” 应该在极小范围内变化才合理。图 7 是不同初始条件下方程的积分结果，看起来很平坦，从而证明积分的结果很可靠。 -10 -8 -6 -4 -2 0 2 ᵪỠ㜭/J 5° 30° 60° 90° 0 1 2 3 4 5 ᰦ䰤/s 图 7 不同条件下的积分常数如果想了解积分常数变化的细节，利用能量之差是有效的方法。能量之差定义为 ∆E = E(t) − E(0) (7) 如果没有事先指定等比例 (用 axis equal 命令表示 x 和 y 轴等比例)，Matlab 在画图时坐标轴会自动根据数据范围进行缩放，因此能量之差的细节变化就可以反映出来 (图 8)，根据能量之差，可以看出系统机械能的变化范围与积分的允许误差范围同阶，初始角度小时计算误差更小。上述分析表明数值计算是可靠的，然后再考虑系统有阻尼的情况。图 9 显示了不同阻尼情况下摆角的变化情况，初始摆角均是 30◦，其中 β = n/(2m) 是阻尼系数。图 10 显示了阻尼对系统机械能的影响，最后能量趋于平坦，为系统静止时的机械能。图 10 中一个细节是：阻尼系数较小时，机械能呈现台阶状的下降，这是因为摆角每次到最大幅值时速度接近， -8 -6 -4 -2 0 2 㜭䟿ѻᐞ/J 0 1 2 3 4 5 ᰦ䰤/s 5° 30° 60° 90° Τ10-8 图 8 能量之差 0 2 4 6 8 1210 ᰦ䰤/s 䀂ᓖ/(°) 30 20 10 0 -10 -20 -30 β// β// β// 图 9 阻尼对摆角的影响 0 2 4 6 8 10 12 ᰦ䰤/s ᵪỠ㜭/J -8.0 -8.5 -9.0 -9.5 -10.0 β// β// β// 图 10 阻尼对机械能的影响

596力学与实践2021年第43卷此时空气阻力很小，机械能接近守恒，所以机械能在系统的运动是周期的吗？在式（7）中消去，有下降过程中就形成了台阶。(mA+mBsin?)lp+mlpsinpcos+(10)(mA+mB)gsinp=03椭圆摆的相关问题椭圆摆是动力学中的一个典型问题，我们想了可以看出，摆角的方程（9）在小角度、线性化后是周解一下：椭圆摆的运动是周期的吗？其周期是多少？期的：类比单摆，大角度情况下摆动周期与初始条件有关。下面的计算只改变椭圆摆的初始角度o，其案例2：圆摆（图11）由质量为mA的滑块A和质量为mB的单摆B构成。滑块可沿光滑水平面他参数均不变。滑动，AB杆长为l，质量不计。试建立系统的运动mi =8 kg, m2= 2 kg(11)微分方程，并求水平面对滑块A的约束力。=0,=0=0图13显示了初始条件对摆角的影响，看起来不AO同条件下摆角都是周期函数，但是周期不同。图14T进一步显示了周期和初始条件的关系，同时表明：初始角为小量时与线性化的结果比较接近。100E5030°图11椭圆摆模型5090()/巢系统有2个自由度，建立Ory坐标系，取A点n位移和AB杆的相对转角为广义坐标（图12)，传统分析可得到系统的运动微分方程以及压力的表-50达式(过程略）为(mA +mB)+mBl(cos-?sin) =0-1001283450l+icos+gsin=0时间/sN=(mA+mB)g+mBl(sinp+p2cos0)) (9)图13椭圆摆的摆角曲线2.3Tro数值计算得到的周期2.2线性化得到的周期A2.12.0TmagafD1.9tTaoBa1.8TmBgmsg1.7图12椭圆摆的受力图020406080初始角度/(°)通常理论力学教材或教学到了这一步就算结束图14椭圆摆的摆动周期了。椭圆摆在运动过程中A滑块和AB杆是否为周期运动？是否为正弦运动？压力是如何变化的？类似图15和图16表明位移和压力也是周期函数，这样的问题传统方法都不好回答。但是初始大角度时位移和压力曲线已经明显偏离标(C)1994-2023ChinaAcademicJournalElectronicPublishingHouse.All rightsreserved.http://www.cnki.net

596 力学与实践 2021 年第 43 卷此时空气阻力很小，机械能接近守恒，所以机械能在下降过程中就形成了台阶。 3 椭圆摆的相关问题椭圆摆是动力学中的一个典型问题，我们想了解一下：椭圆摆的运动是周期的吗？其周期是多少？案例 2：椭圆摆 (图 11) 由质量为 mA 的滑块 A 和质量为 mB 的单摆 B 构成。滑块可沿光滑水平面滑动，AB 杆长为 l，质量不计。试建立系统的运动微分方程，并求水平面对滑块 A 的约束力。 A B 图 11 椭圆摆模型系统有 2 个自由度，建立 Oxy 坐标系，取 A 点位移 x 和 AB 杆的相对转角 φ 为广义坐标 (图 12)，传统分析可得到系统的运动微分方程以及压力的表达式 (过程略) 为 (mA + mB)¨x + mBl( ¨φ cos φ − φ˙ 2 sin φ) = 0 lφ¨ + ¨x cos φ + g sin φ = 0 } (8) N = (mA + mB)g + mBl( ¨φ sin φ + ˙φ 2 cos φ) (9) mAg mBg mBg vr ar ar ae ve y O A x T N B B ϕ n τ 图 12 椭圆摆的受力图通常理论力学教材或教学到了这一步就算结束了。椭圆摆在运动过程中 A 滑块和 AB 杆是否为周期运动？是否为正弦运动？压力是如何变化的？类似这样的问题传统方法都不好回答。系统的运动是周期的吗？在式 (7) 中消去 x¨，有 ( mA + mB sin2 φ ) lφ¨ + mBlφ˙ 2 sin φ cos φ + (mA + mB)g sin φ = 0 (10) 可以看出，摆角的方程 (9) 在小角度、线性化后是周期的；类比单摆，大角度情况下摆动周期与初始条件有关。下面的计算只改变椭圆摆的初始角度 φ0，其他参数均不变。 m1 = 8 kg, m2 = 2 kg x0 = 0, x˙ 0 = 0, φ˙ 0 = 0 } (11) 图 13 显示了初始条件对摆角的影响，看起来不同条件下摆角都是周期函数，但是周期不同。图 14 进一步显示了周期和初始条件的关系，同时表明：初始角为小量时与线性化的结果比较接近。 5° 30° 90° 0 1 2 3 4 5 ᰦ䰤/s -100 -50 50 100 䀂ᓖ/(°) 0 图 13 椭圆摆的摆角曲线 0 20 2.3 2.2 2.1 2.0 1.9 1.8 1.7 40 60 80 (°)/䀂ᓖ࿻ࡍ ઘᵏ/s ᮠ٬䇑㇇ᗇࡠⲴઘᵏ 㓯ᙗॆᗇࡠⲴઘᵏ T T T 图 14 椭圆摆的摆动周期图 15 和图 16 表明位移和压力也是周期函数，但是初始大角度时位移和压力曲线已经明显偏离标

第4期高云峰：Matlab求解理论力学问题系列（三）单摆和圆摆的运动及周期597准正弦曲线了。如果没有数值计算，这些细节无法而在其他位置均为0：因此在频域图中可以直接得出得到。测试函数(11)中的所有参数(图18)。2.00.51501.530°0.4A90°1.0/3赠0.37.0.5u/a0.20AV0.50.1VV-1.00.51.01.52.00230145时间/s时间/s图17测试函数的时域曲线图15椭圆摆的位移曲线1.0160Ai 5°0.2A.30°0.814090RI≥12010.410.610.8A.0.2100A20.280021345510150时间/s频率/Hz图16椭圆摆的压力曲线图18测试函数的频域曲线另外还有一个现象，系统的位移、摆角、压力量下面对椭圆摆的摆角和位移进行快速傅里叶变然都是周期函数，但是周期并不相同。虽然椭圆摆换变化（位移通过平移去掉直流分量），为了方便比的摆角看上去像余弦曲线（图13)，但是否为标准的较不同初始条件的结果，把快速傅里叶变换变化后的结果归一化处理：最大值取为1，处理后的结果见余弦曲线？这可以使用快速傅里叶变换来分析。在图19和图20，可以发现：Matlab中，可以直接调用ft（data_name）函数求出(1）可以看到各条曲线都有一个明显的峰值，但数据data_name的频率。其他位置还有连续的不全为0的值，且大角度时高在处理椭圆摆之前，先看看的效果，设测试倍频处还会有小峰值，这说明椭圆摆的摆角不再是函数为标准的余弦函数；e(t) = Ao + )A;sin(2元fit)(12)(2)由于峰值相对明显，所以在时域图中看起来i=1很像余弦曲线：其中Ao=0.5,A1=1,A2=0.4,A3=0.2;f1=3,(3)各曲线最大峰值对应的主频率不同：频率随f2=5，f3=10.6，其时域图（图17)看不出什么规初始角度增加而减少，或周期随初始角度增加而增律，但经过快速傅里叶变换变化后，在频域图中就可加，符合图14的结论，这反映了非线性函数的周期以清楚看出只会在f=fi处有一个孤立的峰值Ai，或频率与初始条件有关。(C)1994-2023ChinaAcademicJournalElectronicPublishingHouse.All rightsreserved.http://www.cnki.net

第 4 期高云峰：Matlab 求解理论力学问题系列 (三) 单摆和椭圆摆的运动及周期 597 准正弦曲线了。如果没有数值计算，这些细节无法得到。 0 1 2 3 4 5 时间/s 0 0.1 0.2 0.3 0.4 0.5 x/m 5° 30° 90° 图 15 椭圆摆的位移曲线 0 1 2 3 4 5 ᰦ䰤/s 80 100 120 140 160 N/࣋঻ 5° 30° 90° 图 16 椭圆摆的压力曲线另外还有一个现象，系统的位移、摆角、压力虽然都是周期函数，但是周期并不相同。虽然椭圆摆的摆角看上去像余弦曲线 (图 13)，但是否为标准的余弦曲线？这可以使用快速傅里叶变换来分析。在 Matlab 中，可以直接调用fft (data name) 函数求出数据 data name 的频率。在处理椭圆摆之前，先看看fft 的效果，设测试函数为 θ(t) = A0 + ∑ 3 i=1 Ai sin(2πfit) (12) 其中 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)。 0 0.5 1.0 1.5 2.0 ᰦ䰤/s -1.0 -0.5 0 0.5 1.0 1.5 2.0 ᑵᓖ/m 图 17 测试函数的时域曲线 0 5 10.4 0.2 0.1 0 10.6 10.8 ᑵᓖ/m 仁⦷/Hz 10 15 1.0 0.8 0.6 0.4 0.2 0 A A A A 图 18 测试函数的频域曲线下面对椭圆摆的摆角和位移进行快速傅里叶变换变化 (位移通过平移去掉直流分量)，为了方便比较不同初始条件的结果，把快速傅里叶变换变化后的结果归一化处理：最大值取为 1，处理后的结果见图 19 和图 20，可以发现： (1) 可以看到各条曲线都有一个明显的峰值，但其他位置还有连续的不全为 0 的值，且大角度时高倍频处还会有小峰值，这说明椭圆摆的摆角不再是标准的余弦函数； (2) 由于峰值相对明显，所以在时域图中看起来很像余弦曲线； (3) 各曲线最大峰值对应的主频率不同：频率随初始角度增加而减少，或周期随初始角度增加而增加，符合图 14 的结论，这反映了非线性函数的周期或频率与初始条件有关