一维振子运动的数值求解研究

张定梅

摘要：建立了一维振子运动的牛顿运动学方程，利用MATLAB软件数值研究了振子的动力学特性.通过调节参数实现了简谐振动、阻尼振动、受迫振动及复杂非线性振动.给出了各种振动的运动时间序列，讨论了参数对振子的振动周期及振幅的影响，并使用功率谱和吸引子图分析了振子产生的混沌的特性.与传统的解析求解方法相比，数值求解方法更加快捷、直观，避免了复杂的数学积分运算，并且通过运动轨迹图可快速掌握振子的运动规律.

关键词：一维振子;MATLAB;Runge-Kutta算法;混沌

中图分类号：O322  文献标识码：A  文章编号：1673-260X（2019）09-0019-03

1 引言

机械振动是自然界中普遍存在的运动形式，也是大学物理课程中一项重要的教学内容.对于简单的简谐振动可以利用数学积分的方法给出解析解[1]，但对于存在复杂外部扰动的非线性振动系统，解析求解将变得非常困难，因此数值求解将是一种好的替代方式.关于振动问题数值求解的研究已见有关报道[2-4]，例如，文献[2]利用拉格朗日方法建立了复雜的六弹性振子的二维运动方程，并用maple软件进行了数值求解，给出振子运动时间序列图.文献[3]研究了各种组合弹簧振子中的非线性振动问题，给出了数值解和解析解并进行了对比.文献[4]研究了Duffing振子和Van der Pol振子耦合之后的非线性的动力学行为，利用Simulink仿真给出了不同频率及耦合系数下耦合振子的相图和庞加莱截面图.在大学物理教学中，一维振子运动是一种最基本的教学内容，为了能够系统的、直观的研究一维振子在各种外部扰动下的动力学特性，本文建立了一维振子的牛顿运动学方程，采用四阶Runge-Kutta算法并利用MATLAB软件进行数值求解，通过合理调节参数，实现了简谐振动、阻尼振动、受迫振动及复杂非线性振动，给出了振子的运动时间序列图，并且得到了有趣的混沌现象.

2 四阶Runge-Kutta算法

四阶Runge-Kutta算法是求解常微分方程及方程组常用的数值解法，它基于Taylor展开，要求被求解的函数具有较好的光滑性，并且其截断误差为o（h5），具有较高的求解精度.四阶Runge-Kutta公式的一般表达形式为：

其中：

3 一维振子运动学方程

在实际生活中，一维振子除了受到回复力作用以外，往往还受到阻力和外部扰动，为了全面的描述一维振子的运动情况，本文建立的一维振子运动学方程如下[5-7]：

其中，m为振子的质量，-k1x为线性回复力， 为与速度成正比的阻力，-k3x3为与位移有关的三阶非线性阻力项，F1cos1t+F2sin2t为周期性的外力.为简化运算，m取1kg，并且将二阶常微分方程化为一阶微分方程组求解，具体表达式如下：

4 仿真结果及结果

本文采用四阶Runge-Kutta算法对方程（4）进行数值求解，求解步长为0.01s，采用MATLAB软件进行仿真.

4.1 简谐振动

简谐振动是振子最简单的运动形式，此时的振子只受到回复力的作用，方程（4）的参数取值情况如图题所示.图1给出了振子在不同的回复系数下的振动曲线，如图1（a）所示，当A=1s-2时，此时振子的运动轨迹为一规则正弦曲线，振幅为0.110m，振动频率为0.16Hz.当A=3s-2时（图1（b）），由于回复力的增大，振幅变为0.063m，振动频率变为0.28Hz.当A=6s-2时（图1（c）），振幅减小为0.045m，振动频率增大为0.39Hz.从图中可以直观地看出随着回复系数的增大，振子的振动振幅慢慢缩小，同时振动频率慢慢变大.

4.2 阻尼振动

阻尼振动是指振子受到外界阻力下的运动，本文中考虑阻力与振子的运动速度成正比的情况.在阻尼振荡下，方程（4）的参数取值情况如图题所示.图2（a）给出了阻尼系数B=0.01s-1时振子振动的时间序列，可以看出，由于此时阻尼系数较小，振子的振幅衰减较慢.当继续增加阻尼系数到0.03s-1，如图2（b）所示，振子振动到200s时振幅几乎为零.继续增大阻尼系数到0.06s-1（图2（c）），振子到150s就几乎停止振动了.因此，从图中可以直观地看出，随着阻尼系数的增大，振子的振幅衰减的越来越快.

4.3 受迫振动

受迫振动指的是振子在外部周期力作用下的运动，振子在不受外力作用时，由于自身恢复力的作用会产生一个固有振荡频率，此频率与回复系数有关，本文中回复系数A取2s-2，相应的本征角频率为1.414Hz.图3考虑了三种不同的外部扰动频率下振子的运动情况.图3（a）为D=0.5ms-2，ω1=1Hz情况下的受迫振动振动曲线，此时外部扰动频率与振子的固有频率不相同，但差距较小，振子在72s后趋于规则振荡，振荡幅度为0.52m.当D=0.5ms-2，ω1=1.414Hz时（图3（b）），此时外部扰动频率与振子的固有频率相同，出现了共振，振子的振荡幅度在规则振荡后达到了4.2m.当D=0.5ms-2，ω1=3Hz时（图3（c）），此时外部扰动频率与振子的固有频率差距较大，稳定后的振荡幅度只有0.07m，未形成共振.从图中可以清楚地看出在不同的外部扰动频率下振子振幅及频率的变化.

4.4 非线性振动

图4给出了振子在非线性力作用情况下的振动曲线图及功率谱图.其中功率谱指的是将时域中的时间序列通过傅里叶变换转换到频域进行观察，它代表了该时间序列中某一频率间隔内的功率占总功率的比例，能够直观地看出时间序列中所包含的频率成分.如图4（a）所示，振子的运动轨迹为规则的脉冲包，时间序列（图4（a））中虽然有多个峰值，但整体的振动是重复的，功率谱（图4（a1））中可以看出多个频率成分，其中最强的峰值为振动的主频率（约为0.2GHz），主频率前面的峰值为次谐波.图4（b）所示，此时振子的运动出现混沌行为，时间序列（图4（b））无明显周期行为，并且出现了随机的振动，功率谱（图4（b1））出现了展宽并无明显的峰值.

混沌是一种有趣的物理的现象，它具有自相似性、初值敏感性和内在随机性等特性.为了进一步研究所产生的混沌的性质，图5（a）给出了不同初始条件下的振子的运动轨迹，从图中可以看出，即使所有参数都一样，当初始条件不同，两条轨迹在27s的时候出现了偏离，而且随着时间的增加，偏离越来越大，最后轨迹完全不相同，这就证明了混沌具有初值敏感性.图5（b）给出了该参数情况下的时间序列的吸引子图，图中有多个环路缠绕在一起，并且具有一定的自相似性，表现为奇怪吸引子.

5 总结

本文采用四阶Runge-Kutta算法数值研究了一维振子的运动问题，利用MATLAB软件进行了仿真模拟，通过调节参数实现了简谐振动、阻尼振动、受迫振动及复杂非线性振动，给出了各种振动的运动时间序列图.从图中可以直观地看出，对于简谐振动，回复系数越大，振子运动的振幅越小.对于阻尼振动，随着阻尼系数的增大，振子的振幅衰减的越来越快.

对于受迫振动，当外部扰动的频率与振子本征频率相同时，振子的振幅将大幅度增加.最后对于非线性振动情况，适当调节参数可实现混沌振动.本文简单直观地展示了振子在各种情况下的运动规律，可作为大学物理教学内容的有效补充.

——————————

參考文献：

〔1〕王树平.对称非线性弹簧振子的振动周期的级数解[J].大学物理，2017，36（4）：15-16.

〔2〕包兴明，袁玉全，闫安英，戚作涛对称六弹性振子的二维非线性振动[J].大学物理，2010，35（2）：1-6.

〔3〕廖旭，任学藻.组合线性弹簧振子中的非线性振动[J].大学物理，2008，27（2）：25-28.

〔4〕王晓东，杨绍普，赵志宏.Duffing振子和Van der Pol振子耦合的动力学行为分析[J].石家庄铁道大学学报（自然科学版），2015，28（4）：53-57.

〔5〕卓崇培.非线性物理学[M].天津：天津科学技术出版社，1996.26-61.

〔6〕刘式适，刘式达.物理学中的非线性方程[M].北京大学出版社，2000.6-66.

〔7〕倪致祥，马涛.一种非简谐的微振动模型[J].大学物理，2006，25（9）：14-16.