单摆运动方程及其周期近似解

单摆运动周期和相轨迹特性的研究摘要：本文首先分别用基本形式拉格朗日方程和保守系下的拉格朗日方程求解单摆运动方程，其次通过线性近似法求解在小摆角下的周期近似解，再通过构建“局部常化”的近似处理方法，得到大摆角运动周期的一个新结论。最后，用数值模拟（四阶龙格-库塔法）求解无阻尼无驱动单摆非线性方程，用origin作图软件绘制出90时，取不同初值时的相轨迹，并分析了其相轨迹特性，验证对小角度单摆几乎只有摆动，对大角度单摆既有摆动又有转动。关键词：单摆运动周期非线性局部常化椭圆积分数值模拟相轨迹引言：非线性引起复杂性，复杂性产生的根源即“原来是禁锢在笼子里的非线性老虎被释放了”。对线性模型简单、容易分析，且线性微分方程可求其解析解，而非线性模型复杂、不容易分析，非线性方程不容易求其解析解，我们利用这两种性质可对某一具体问题进行不同方式的分析，得出一部分规律。单摆模型是简单与复杂的综合体，对该模型可用：线性化法、近似解析法、数值法和向空间法进行求解，分析。本文求出单摆微分方程后，首先通过线性近似法求解在小摆角下的周期近似解，再通过构建“局部常化”的近似处理方法，得到大摆角运动周期的一个新结论，有用数值法和相空间法，验证了单摆运动在取不同初值时的运动形态，即为摆动或摆动加转动特性，对单摆特性研究有一定价值。除了对无阻尼无驱动单摆系统研究，我们可将该分析方法用与其他几类单摆模型。正文：1单摆运动方程的求解单摆运动问题是一个古老而又十分有趣的问题。对于摆长为L，最大摆角为0的单摆系统，由于只有重力做功，因此满足机械能守恒。分别用基本形式拉格朗日方程和保守系下的拉格朗日方程来求解如下：（1）基本形式拉格朗日方程为1：=12.dTTQdtqq（，，）（1）自由度为1，取广义坐标为，有：广义力为：qrgmqrFQiiii)(2121mglsinθFjlsinθilcosθθrjlcosθilsinθr222iAAlmmvT2220122dTTQdtqqTTmlmldTmldt代入基本拉格朗日方程，得2sin0mlmglsin0gl（2）保守系下的拉格朗日方程为0=12.dLLdtqq（，，）（2）自由度为1，取广义坐标为，有：222222222111222cos1cos2sinTmvmlmlVmglLTVmlmglLmldLdmlmldtdtLmgl代入到（1）式中，有0sinlg令lg20，则有0sin20（3）（3）式是一个非线性微分方程，而大多数非线性微分方程都很难找到其解析解，这给动力系统的分析带来了很大的困难。再者，非线性系统能产生“混沌”现象，其解析解通常也是非常复杂的。文献1和文献2分别从机械能守恒定律和相图关系求出了精确的单摆运动周期公式：近似为sin，则（2）式可简化为200w，对其两边乘以2d，然后积分：202(/)22sin0ddtdwddt得220()2sindwcdt即220()2sindwcdt其中c为常数积分。当摆动到最大角度0时，0ddt。所以2002sincw，因此，002(coscos)dwdt分离变量并积分002(coscos)dwdt如果t=0时，0，并设T是单摆的振动周期，则t=T/4时，0，所以00022042sinsin22Tdw（4）令0sinsinsin22，-1sin1。两边积分，得：1cossincos222dd02sincos2cos2dd则（4）右边被积函数写为：0222022202sincoscos2cos1sin2cossinsin2222cos1sin1sinsin222ddddd积分限位0时00;时，2上式求解有：2022002.41sinsin2TdT202020sin2sin12dTT（5）其中0是单摆的最大摆角。式（5）适合于任意摆角下的单摆运动，但这个公式是用完全椭圆公式表示的，过于复杂，应用时需要查椭圆积分表，因而实用性不强。本文先通过线性近似求出小角度下的单摆等时公式，其次通过构建“局部常化”的近似处理方法，给出在0≦90时近似度较好的一个单摆运动周期解。2小摆角下的周期近似解在单摆运动系统中，如果摆角很小（一般0≦5）时，可将单摆运动近似成一种简谐振动。对于运动方程（3），可以做如下近似：由于摆角很小，所以可以将非线性因子sin作一级近似，将方程（3）转化成线性方程。即有：sin，所以（3）式可写成：020（6）该方程的解为：tAtAtA00201coscossin所以，lgT2200现在，我们将这种近似周期公式0T与精确的周期公式进行比较，其结果如表1示：（T由椭圆积分表查7得）表1：0T与T的相对误差对比0346102030TT00.99980.99970.99930.99810.99250.98280405060708090TT00.96960.95260.93180.90740.87910.8427从表中可以看到，当最大摆角较小时，由0T计算得到的周期相对误差较小，而当最大摆角较大时，相对误差较大。3“局部常化”的近似周期解上面，我们将sin作了一级近似，将非线性微分方程转换成了线性微分方程，但通过表1我们看到，这种近似是很粗略的，当大0较大时产生的误差比较大。现在，我们需通过构建“局部常化”的近似方法来给出近似度较好的单摆运动周期解公式。对于方程（3），我们可以采用局部常化4，将sin作如下变换：0cos2cos2cos222cos2sin2sin（7）在这里，我们将变量2cos视为常量0cos。其中0,0。现在，将（7）式代入（3）式得0222dtd（8）其中，000coscosLg。式（8）是对动力学方（3）的有一个修正，将一个分线性问题转换成一个线性问题。由（6）式我们可以很容易的得到单摆的周期公式0001coscos2TgLT（9）式中的是一个修正常数，修正的方发事将式（6）与标准式（3）在取不同角度时进行比较。通过取的不同值进行尝试比较，我们发现修正值为=0.496时1T有比较高的精度。下面我们将这种近似公式1T与精确的周期公式进行比较，结果如下表2所示.表2：1T与T的相对误差对比0346102030TT11.00001.00001.00001.00001.00000.99980405060708090TT10.99980.99981.00001.00081.00211.0043表2的数据表明，当900时，公式1T的近似度很好，相对公式0T来说近似度明显的提高了。4数值模拟分析相轨迹运动微分方程：{vA)0()0(sin20数值编程：令t=x,y,z则微分方程转换为：{vxzAxyywzzy)()(sin0020用四阶龙格-库塔法5程序如下：Functionf(yAsDouble)AsDoubleg=9.8l=0.6f=-g*Sin(y)/lEndFunctionPrivateSubCommand1_Click()DimgAsDouble,lAsDouble,nAsInteger,hAsDouble,aAsDouble,bAsDoubleDimtAsDouble,t0AsDouble,t1AsDouble,y0AsDouble,y1AsDouble,z0AsDouble,z1AsDoubleDiml1AsDouble,l2AsDouble,l3AsDouble,l4AsDoublea=0:b=2n=800h=(b-a)/nOpen时间.txtForOutputAs#1Open角度.txtForOutputAs#2Open角速度.txtForOutputAs#3y0=20:z0=6:t0=aFort=aTobStepht1=t0+hl1=f(y0)l2=f(y0)l3=f(y0+(h^2*l1)/4)l4=f(y0+h^2*l2/2)y1=y0+h*z0+h^2*(l1+l2+l3)/6z1=z0+h*(l1+2*l2+2*l3+l4)/6Text1.Text=Text1.Text&t&&Format(y1,##.######)&&Format(z1,##.######)&vbCrLfWrite#1,tWrite#2,y1Write#3,z1y0=y1z0=z1NexttClose#1,#2,#3EndSub按理论：（1）取n=800,t=0-2s,g=9.8kg^2/s,l=0.6m。初值：200y,z0=0、1、2、3、4、5rad/s同一横坐标z1，同一纵坐标y1模拟相轨迹如下：-8-6-4-20246816182022###Y1AxisTitleZ1AxisTitleB205-8-6-4-20246816182022###Y1AxisTitleZ1AxisTitleB204-8-6-4-20246816182022###Y1AxisTitleZ1AxisTitleB203-8-6-4-20246816182022###Y1AxisTitleZ1AxisTitleB202-8-6-4-20246816182022###Y1AxisTitleZ1AxisTitleB201-8-6-4-20246816182022###Y1AxisTitleZ1AxisTitleB200分析：由图可知初值在角度y一定的情况下，随着初角速度z的增大，相轨迹椭圆越来越大，并且在在初值为5时，椭圆出现缺口。由于ywzdzdysin20，所以Eyzcos212，即相轨迹方程。总能量为0cos0212zyE初值y不变，初值z增大后，摆能量增大，对应椭圆面积较大，图示符合理论，当初始能量越大，单摆由摆动转为振动，如果继续取值。椭圆缺口会越来越大，最后形状变化无规律可寻。（2）取n=800,t=0-2s,g=9.8kg^2/s,l=0.6m。初值：z0=0rad/s，y0=1、20、30、50、65、75度，同一横坐标z1，模拟相轨迹如下：-8-6-4-202468-2.0-1.5-1.0-0.50.00.51.01.52.0Y1AxisTitleZ1AxisTitleB10-8-6-4-20246816182022###Y1AxisTitleZ1AxisTitleB200-8-6-4-202468293031323334###Y1AxisTitleZ1AxisTitleB300-8-6-4-20246850.050.150.250.350.450.550.6###Y1AxisTitleZ1AxisTitleB500-8-6-4-20246860626466Y1AxisTitleZ1AxisTitleB650-8-6-4-20246875.075.275.475.675.8Y1AxisTitleZ1AxisTitleB750分析：由图可知在初角速度z一定的情况下，随着角度y的增大，相轨迹椭圆变化较明显。初值为1度时，椭圆近似为圆；当初值小于30度时，椭圆纵横向都增大，符合能量公式；当初值为5