10.2 分析单摆的运动
1.单摆运动的微分方程和数值算法
（1）微分方程
最基本的单摆是由绝对挠性且长度不变、质量可以忽略不计的一根绳子或一根竿和一个锤组成的。在绳或竿一端的下方系着锤，绳或竿的另一端固定，锤来回移动，假设单摆绳子的长度为l，锤的质量为m，摆动角度为θ，则考虑空气阻尼〔所受的阻尼为
（q为阻尼系数）〕时，可由牛顿第二定理得到描述单摆运动的微分方程：
其中θ0和v0分别是初始摆角和初始驱动速度。
当q=0时，方程（10.11）是无阻尼单摆运动方程。在无阻尼且摆角很小（0＜θ≤5°）的情况下，可以假设θ≈sinθ，此时单摆运动方程变为二阶常系数线性微分方程，可以确定其解析解。但是，在一般情况下，方程（10.11）没有解析解，因此要用适当的数值方法求其数值解以进行分析。
（2）数值算法
设
=v，则方程（10.11）变为如下二元一阶非线性常微分方程组：
其中，θ是位置，v是单摆的摆动速度。
确定常微分方程（10.12）的数值解时，可以用四阶龙格-库塔法，对应方程（10.12）的四阶龙格-库塔计算公式为
其中h是时间步长。
2.单摆动态可视化和运动分析方程模板的制作
下面的所有制作过程以图10.7为参考。
第一步，输入已知条件。在单元格A2、B2、C2、D2、E2、F2、T2和H2中分别输入重力加速度、单摆线长度、阻尼系数、初始位置、初始驱动速度、开始时间、模拟终止时间、时间步长参数。在单元格I2中输入计算总次数的公式“=（G2-F2）/H2”；在单元格J2中输入实现循环和计算当前计算步骤的公式“=IF（J2＜I2/2，J2+1，0）”。
图10.7 单摆有关的数据(https://www.daowen.com)
第二步，计算当时位置和速度。在单元格A5、B5和C5中分别输入“=E2”“=D2”和“=E2”；用式（10.13a）及式（10.13b）在单元格D5和E5中分别输入“=C14”和“=-$C$2*C14-$A$2/$B$2*SIN（B14）”；用式（10.13c）及式（10.13d）在单元格F5和G5中分别输入“=C5+$H$2/2*E5”和“=-$A$2/$B$2*SIN（B5+$H$2/2*D5）-$C$2*（C5+$H$2/2*E5）”；用式（10.13e）及式（10.13f）在单元格H5和I5中分别输入“=C5+$H$2/2*G5”和“=-$A$2/$B$2*SIN（B5+$H$2/2*F5）-$C$2*（C5+$H$2/2*G5）”；用式（10.13g）及式（10.13h）在单元格J5和K5中分别输入“=C7+I7*$H$2”和“=-$A$2/$B$2*SIN（B7+$H$2*H7）-$C$2*（C7+$H$2*I7）”；在单元格A6和A7中分别输入前一时刻和现在时刻的计算公式“=IF（J2＜=1，A5+H2，A7+H2）”和“=IF（J2＜=1，A5+2*H2，A6+H2）”；在单元格B6和C6中用式（10.13i）及（10.13j）分别输入前一时刻位置和速度的计算公式“=IF（J2＜=1，B5+H2/6*（D5+2*F5+2*H5+J5），B7+H2/6*（D7+2*F7+2*H7+J7））”和“=IF（J2＜=1，C5+H2/6*（E5+2*G5+2*I5+K5），C7+H2/6*（E7+2*G7+2*I7+K7））”；在单元格B7和C8中分别输入现在时刻位置和速度的计算公式“=B6+$H$2/6*（D6+2*F6+2*H6+J6）”和“=C6+$H$2/6*（E6+2*G6+2*I6+K6）”；选中单元格区域D5：K5，把鼠标指针移动到单元格K5的右下角，当其变成“
”符号时，按住鼠标左键向下方拖动两行。
第三步，计算单摆线的终点坐标和确定画图用的点坐标。在单元格A10、B10、C10和D10中分别输入“=-B2/2”“=B2”“=0”和“=B2”；在单元格A11、B11、C11和D11中分别输入“=B2/2”“=B2”“=B2*SIN（B7）”和“=B2-B2*COS（B7）”。
第四步，绘制单摆动画仿真图。选中单元格区域A10：B11，依次单击“插入”、图表下拉按钮，“更多散点图（
）...”，选择“XY散点图→带平滑线的散点图”，单击“确定”按钮，便可绘制出点（-5，10）与点（5，10）连接的横线图像。然后在一个直角坐标系内使用画多个函数图像的方法（参考图3.3～3.5）分别绘制单摆竿（单元格区域C10：D11）及单摆锤（单元格区域C11：D11）的图像，并且对图像进行装饰，得到图10.8所示的动画仿真图。
图10.8 动画仿真演示图
第五步，计算位置和速度。在单元格A14、B14和C14中分别输入“=F2”“=D2”和“E2”；选中单元格区域D5：K5，单击右键，选择“复制”，再选中单元格区域D14：K14，单击右键，选择“粘贴”；在单元格A15中输入“=A14+$H$2”；选中单元格区域B6：C6，单击右键，选择“复制”，再选中单元格区域B15：C15，单击右键，选择“粘贴”；拖动填充法自动填充单元格区域A15：K114，如图10.9所示。
图10.9 位置和速度计算结果
第六步，绘制位置和速度图。选中单元格区域A14：B114，依次单击“插入”、图表下拉按钮，“更多散点图（
）...”选择“XY散点图→带平滑线的散点图”，单击“确定”按钮，得到图10.10所示的单摆位置曲线。选中刚绘制的位置曲线图，单击右键，选择“选择数据”，弹出“选择数据源”对话框后按单击“添加”按钮，此时弹出“编辑数据系列”对话框，在系列名称内填“速度曲线”，在x轴系列值内选择填“=Sheet1！$A$14：$A$114”，在y轴系列值内选择填“=Sheet1！$C$14：$C$114”，确定添加速度曲线。利用同样的方法绘制当前位置图，在系列名称内填“当前位置”，在x轴系列值内选择填“=Sheet1！$A$7”，在y轴系列值内选择填“=Sheet1！$B$7”确定添加当前位置。然后绘制当前速度图，在系列名称内填“当前速度”，在x轴系列值内选择填“=Sheet1！$A$7”，在y轴系列值内选择填“=Sheet1！$C$7”确定添加当前速度。最后绘制相图，在系列名称内填“相图”，在x轴系列值内选择填“=Sheet1！$C$14：$C$114”，在y轴系列值内选择填“=Sheet1！$B$14：$B$114”确定添加相图对图表进行装饰，得到图10.11所示的图。
图10.10 单摆位置曲线
图10.11 有阻尼单摆的位置曲线、速度曲线和相图
第七步，在“文件”菜单下选择“选项”，在“选项”对话框中选择“公式”副对话框。勾选“选择启用迭代计算
）”复选框，最多迭代计算次数设为1，最大误差设为0.001，单击“确定”按钮。这时，如果我们按功能键F9，则图10.8所示的单摆开始摆动，单摆在图10.11中的当前位置并以当前速度沿着位置曲线和速度曲线移动时，很容易观察出单摆在各个时刻的速度和位置。
这个模板巧用Excel绘制图像和实现动画功能“数”、“形”及“动”连贯形象直观地演示了单摆运动的动态模拟仿真；实现了交互性，建立了不同参数情况下研究单摆运动和模拟的平台，得出了单摆无阻尼、单摆有阻尼、摆角为任何角度值时速度和位置的数值解，以及速度和位置的曲线和相应的相图。在这个模板内可以改变重力加速度、单摆线长度、阻尼系数、初始位置、初始驱动速度、开始时间、模拟终止时间、时间步长等参数，最终可得出相应的演示模板。例如，在这个模板内阻尼系数改为零，就可以得到无阻尼单摆动画演示模板，图10.12是无阻尼单摆位置曲线、速度曲线和相图。
图10.12 有阻尼单摆的位置曲线、速度曲线和相图
