资源描述：
《用Matlab软件求单摆的运动.pdf》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、2005年赣南师范学院学报l.3第三期JournalofGannanTeachersCollegeJune.2005用Matlab软件求单摆的运动X罗颖,罗兴垅(赣南师范学院物理与电子信息科学系,江西赣州341000)摘要:直接从单摆的运动微分方程出发,应用Matlab软件,求出了线性单摆阻尼振动的解析解和非线性单摆无阻尼振动、阻尼振动的数值解;画出了单摆的振动曲线和相图.关键词:单摆;阻尼;非线性;数值解;Matlab软件中图分类号:O439文献标识码:A文章编号:1004-8332(2005)03-0090-031线性单摆的
2、无阻尼振动一根不会伸缩、长度为l的细线,上端固定(或一根刚性轻杠,上端与无摩擦的铰链相连),下端悬挂一质量为m的小球就构成一单摆.当摆角H很小时,sinHU0,单摆作简谐振动,令s=lH,其运动方程为s=Acos(X0t+<0)(1)式中X0=g/l,称为单摆的固有圆频率,其周期为T0=2P/X0=2Pl/g(2)式(1)对时间t求一阶导数,得小球的速度v=-AX0sin(X0t+<0)(3)由式(1)与式(3)得22sv+=1(4)AAX0上式表明单摆的速度-位移曲线图(通常称为相图)为一椭圆.2线性单摆的阻尼振动当小球摆动的
3、速度较小时,小球受到一个与速度方向相反的阻力f=-Cds/dt,C为阻力系数,它与物体的形状以及周围媒质的性质有关.根据牛顿第二定律有2dsds22+2B+X0s=0(5)dtdt式中B=C/2m,称为阻尼因数.对于一定振动系统,根据比值(BBX0)小于、等于、大于1,则称单摆处于欠阻尼、临界阻尼和过阻尼振动状态.设振动系统的周期T0=2s,初始条件s(0)=0.10m,v(0)=0.用Matlab软件中的函数/dsolve0与/ezplot0可求出单摆阻尼振动方程并描绘振动曲线:(Ñ)当BBX0=0.3BP时,单摆作欠阻尼振动
4、.仅需运行以下程序:clear;clf;symsst;*s=dsolve('D2s+pi^2s=0','s(0)=0.10,Ds(0)=0'),v=diff(s,t,l),ezplot(s,[0,6,-0.1,0.1],1),ezplot(s,v[0,6],2)即可得到单摆的阻尼振动(见图1)和相图(见图2),振动方程(略).从图1、图2可以看出,线性单摆的阻尼振动周期T大于无阻尼的振动周期(T0=2s),随着振幅的指数衰减,相轨迹向内卷缩,呈螺旋状,最后趋于中心的不动点(稳定平衡位置).该点是动力学系统的一个吸引子(attra
5、ctor).(Ò)当BBX0=1时,单摆作临界阻尼振动.运行程序:X收稿日期:2004-05-21修回日期:2005-04-12作者简介:罗颖(1980)),女,江西信丰人,赣南师范学院物理与电子信息科学系教师,主要从事物理实验的教学与研究.第3期罗颖,罗兴垅用Matlab软件求单摆的运动91***s=dsolve('D2s+2piDs+pi^2s=0','s(0)=0.10,Ds(0)=0')ezplot(s,[0,12,0,0.1],1)图1单摆的欠阻尼振动曲线图2单摆的欠阻尼振动的相图可得单摆的临界阻尼振动曲线见图3曲线1
6、,振动方程(略).(Ó)当BBX0=9BP时,单摆作过阻尼振动.运行程序:***s=dsolve('D2s+29Ds+pi^2s=0','s(0)=0.10,Ds(0)=0')ezplot(s,[0,12,0,0.1],2)可得单摆的过阻尼振动曲线见图3曲线2,振动方程(略).从图3可看出,单摆作临界阻尼、过阻尼振动时,分别经时间t临=1.25T0=2.5s,t临=4.5T0=9s,其振幅趋于零,即单摆静止.3非线性单摆的无阻尼振动根据牛顿第二定律,单摆的运动微分方程为2dH2图3单摆的临界阻尼与2+X0sinH=0(6)dt过
7、阻尼振动曲线PdH上式是一个非线性的二阶微分方程.设X0=Prad/s,初始条件为H
8、t=0=,
9、t=0=0.下面介绍用Matlab2dt[1]HH0软件的函数/ode450求数值解的方法(非线性单摆无阻尼振动的解的解为sin=ksn(X0t,k),k=sin,222T0H0为单摆的振幅,sn为Jacobi椭圆正弦函数,振动周期T=K(k),K(k)为第一类完全椭圆积分,当H0PP=时,T=1.180T0).2HdH(Ñ)首先把高阶微分方程改写成一阶微分方程组令y1=H,y2=,y=dHdtdt于是原二阶微分方程可改写成如下一阶
10、微分方程组dy1Pdty2y1(0)=,初始条件为=2dy22-X0siny1y2(0)0dt(Ò)然后根据上述一阶微分方程组编写M函数文件(即ODE文件)%fodelxdb2.mfunctiondydt=fodelxdb2(t,y)*w=pi;dydt=[y(
