用Matlab仿真探究摆角对单摆周期影响

用Matlab仿真探究摆角对单摆周期影响摘要：本文通过Matlab仿真验证小角摆动是简谐振动，并利用数值法求解微分方程，画出不同摆角下单摆的振动图像，定性分析，得出大角摆动时单摆周期随角度的变化情况。关键字：单摆周期摆角大小0.引言单摆是生活中常见的一种简单物理模型，物理学中所讨论的单摆是一种理想化的模型,也称数学摆。它由一根不可伸缩的细线(质量不计),一端固定,另一端悬挂一质量为m的小球(视为质点)，且摆角小于5度的振动系统。对于这种理想单摆的周期，不随摆角大小的改变而改变。但当单摆摆角大于5度时，理想单摆的周期公式不再适用，本文通过建立物理模型，在忽略空气阻力的前提下，用Matlab进行大角摆动的模拟，画出震动图像，研究大角摆动时摆角对周期的影响。1.建立物理模型根据单摆的理想条件，摆线不可伸长且质量忽略不计，空气阻力忽略不计.设摆线长度为l,摆球质量为m，重力加速度为g,摆球离开平衡位置的角度为，对摆球做受力分析如有图所示。由牛顿第二定律，有22ddt=-sinsin2lg（1-1）其中lg假定位移很小，sin，即小角摆动。则式1-1可表示为0lgd22dt（1-2）大角摆动时，仍为1-1式，该式为非线性方程，为方便起见，将用y来表示，该式又可以写为下列一阶微分方程组Tmg21ydtdy；12ysinlg-dtdy（1-3）2.用Matlab方程求解2.1小角摆动用Matlab求解式1-1，其结果为y=theta0*cos(1/l^(1/2)*g^(1/2)*t)（2-1）即ttgcos0（2-2）由式2-1可以看出，当小角摆动时为简谐振动。其周期为：T0=gl22（2-3）取摆长l=1，重力加速度g=9.8，摆角5，利用Matlab求解式1-1，取步长为0.1的点作图，与应用式1-2求出解析解并绘图，将两图象放在一起比较如下图1所示012345678910-0.1-0.08-0.06-0.04-0.0200.020.040.060.080.1图1小角摆动图1中曲线为求解式1-2所得图像，而离散点为求解式1-1所得，观察图像可得两方程的解几乎吻合，可以说明当较小时（5），两方程的解几乎相等，故周期公式此时较为准确，即单摆周期不随摆角变化而变化。上述结论仅仅适用于摆角很小时（5），当摆角很大时，sin不再成立，式1-1与式1-2的解不再相近，故周期公式（2-3）不再成立。下面我们继续讨论摆角比较大时的单摆运动规律。2.2大角摆动用Matlab求解式1-1，前面已经提到为方便起见用式1-3表示。利用数值法求解微分方程，用Matlab分别绘制12，6，3时的振动图像，如下图2所示012345678910-1.5-1-0.500.511.5pi/12pi/6pi/3图2大角摆动观察图2中不同摆角的振动图像，可以看出摆角12时，周期最小，摆角3时，周期最大。从而得出大角摆动时，单摆周期并不像小角摆动时一样，不随摆角变化，而是大角摆动周期随着摆角的增大而增大。结语：对于类似单摆的运动问题，我们需要求解非线性微分方程，然而对于非线性微分方程又很难求出解析解，所以可以应用Matlab求出数值解，再加以分析，从而得出单摆的运动规律。本文通过计算、作图验证了小角摆动是时单摆的振动周期公式是正确的，得出小角摆动周期不随摆角的改变而改变；对于大角摆动，本文则通过数值解并绘图，定性得到大角摆动时，单摆周期随着摆角的增大而增大。参考文献：[1]钞曦旭，MATLAB及其在大学物理课程中的应用[M]，西安：陕西师范大学出版社，2006[2]陈怀琛，MATLAB及其在理工课程中的应用指南（第三版）[M]，西安：陕西电子科技大学出版社，2007附录clear,clcy=dsolve('D2y+g/l*y=0','y(0)=theta0,Dy(0)=0','t')l=1;g=9.8;theta0=5/180*pi;t=(0:0.01:2)*pi;y=eval(y);plot(t,y,t,0)clear,clcy=dsolve('D2y+g/l*y=0','y(0)=theta0,Dy(0)=0','t')l=1;g=9.8;theta0=5/180*pi;t=(0:0.1:10);y=eval(y);[t,u1]=ode45('fangcheng',[0:0.1:10],[5/180*pi,0],[]);plot(t,u1(:,1),'r*',t,y,'b-',t,0,'k-');holdonl=1;g=9.8;[t1,u1]=ode45('fangcheng',[0:0.01:10],[pi/12,0],[]);[t2,u2]=ode45('fangcheng',[0:0.01:10],[pi/6,0],[]);[t3,u3]=ode45('fangcheng',[0:0.01:10],[pi/3,0],[]);plot(t1,u1(:,1),'r-',t2,u2(:,1),'b:',t3,u3(:,1),'k--');holdonfunctiondydt=fangcheng(t,y,flag);dydt=[y(2);-9.8*sin(y(1))]