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

用Matlab 仿真探究摆角对单摆周期影响
摘要：本文通过Matlab 仿真验证小角摆动是简谐振动，并利用数值法求解微分方程，画出不同摆角下单摆的振动图像，定性分析，得出大角摆动时单摆周期随角度的变化情况。

关键字：单摆周期 摆角大小 0.引言
单摆是生活中常见的一种简单物理模型，物理学中所讨论的单摆是一种理想化的模型,也称数学摆。

它由一根不可伸缩的细线(质量不计),一端固定,另一端悬挂一质量为m 的小球(视为质点)，且摆角小于5度的振动系统。

对于这种理想单摆的周期，不随摆角大小的改变而改变。

但当单摆摆角大于5度时，理想单摆的周期公式不再适用，本文通过建立物理模型，在忽略空气阻力的前提下，用Matlab 进行大角摆动的模拟，画出震动图像，研究大角摆动时摆角对周期的影响。

1.建立物理模型
根据单摆的理想条件，摆线不可伸长且质量忽略不计，空气阻力忽略不计.设摆线长度为l,摆球质量为m ，重力加速度为g,摆球离开平衡位置的角度为θ，对摆球做受力分析如有图所示。

由牛顿第二定律，有
2
2d dt θ
=-
θωθsin sin 2-=l
g
（1-1
其中
l
g =
ω 假定位移很小，θθ=sin ，即小角摆动。

则式1-10l
g
d 2
2
=+
θθ
dt （1-2） 大角摆动时，仍为1-1式，该式为非线性方程，为方便起见，将θ用y 来表示，该式又可以写为下列一阶微分方程组
21y dt dy = ；()12y sin l
g
-dt dy = （1-3） 2.用Matlab 方程求解 2.1小角摆动
用Matlab 求解式1-1，其结果为
y =theta0*cos(1/l^(1/2)*g^(1/2)*t) （2-1）
即
⎪⎪⎭
⎫
⎝⎛=t t g cos 0θθ （2-2）
由式2-1可以看出，当小角摆动时为简谐振动。

其周期为：
T 0=
g
l
22π
ω
π
=⨯ （2-3）
取摆长l=1，重力加速度g=9.8，摆角οθ5=，利用Matlab 求解式1-1，取步长为0.1的点作图，与应用式1-2求出解析解并绘图，将两图象放在一起比较如下图1所示
1
2
3
4
5
6
7
8
9
10
图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所示
12345678910
-1.5-1
-0.5
0.5
1
1.5
图2 大角摆动
观察图2中不同摆角的振动图像，可以看出摆角12
π
θ=时，周期最小，
摆角3
πθ=时，周期最大。

从而得出大角摆动时，单摆周期并不像小角摆动时
一样，不随摆角变化，而是大角摆动周期随着摆角的增大而增大。

结语：对于类似单摆的运动问题，我们需要求解非线性微分方程，然而对于非线性微分方程又很难求出解析解，所以可以应用Matlab求出数值解，再加以分析，从而得出单摆的运动规律。

本文通过计算、作图验证了小角摆动是时单摆的振动周期公式是正确的，得出小角摆动周期不随摆角的改变而改变；对于大角摆动，本文则通过数值解并绘图，定性得到大角摆动时，单摆周期随着摆角的增大而增大。

参考文献：
[1]钞曦旭，MATLAB及其在大学物理课程中的应用[M]，西安：陕西师范大学出版社，2006
[2]陈怀琛，MATLAB及其在理工课程中的应用指南（第三版）[M]，西安：陕西电子科技大学出版社，2007
附录
clear,clc
y=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,clc
y=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-');
hold on
l=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--'); hold on
function dydt=fangcheng(t,y,flag);
dydt=[y(2);-9.8*sin(y(1))]。