北京科技大学数学大作业物理组

《Matlab数学实验》大作业
研究课题轻杆单摆振动的周期和运动规律
学院数理学院
班级物理1101
组长周慧斌
组员童鑫周涛
2012年12月
一、应用题目及背景
单摆是质点振动系统的一种，是最简单的摆。

绕一个悬点来回摆动的物体，都称为摆，但其周期一般和物体的形状、大小及密度的分布有关。

但若把尺寸很小的质块悬于一端固定的长度且不能伸长的轻杆（质量可忽略不计）上，把质块拉离平衡位置，放手后质块往复振动，此即轻杆单摆的振动问题。

当轻杆和过悬点铅垂线所成角度小于10°时，可视为质点的振动周期 T只和当地的重力加速度g有关，而和质块的质量、形状和振幅的大小都无关系，其运动状态可用简谐振动公式表示，称为单摆或数学摆。

如果振动的角度大于 10°，则振动的周期将随振幅的增加而变大，就不成为单摆了。

如摆球的尺寸相当大，轻杆的质量不能忽略，就成为复摆（物理摆），周期就和摆球的尺寸有关了。

本文将从严格的物理学角度，应用Matlab软件对各种角度下的轻杆类单摆的周期和运动规律进行模拟分析。

此项研究对物理学其它模型的研究，各种科研活动中精密仪器的研制具有重要指导意义。

此项研究解决的主要问题如下：
1.设想一长为l的轻杆，连接一个质量为m的摆球，形成一个单摆，不计摩擦。

建立单摆的周期与角振幅关系的数学模型，并用matlab 进行曲线模拟.
2.编程演示单摆振动的动画，比较单摆振动和简谐振动的规律。

3.当单摆角振幅的度数为1°到7 °时(间隔为1 °)，将单摆运动的角位置和角速度与简谐振动进行比较。

当单摆角振幅的度数为
30 °到150 °时(间隔为30 °)，另加179 ° ，同样进行比较。

二、数学模型的建立
如图1所示，设角位置为θ，摆锤的运动方程为
即：
在小角度的情况下，sinθ ≈ θ，可得
ω0为圆频率： 振动周期为：
图1 可见：在小角振动的情况下，单摆的周期与角振幅无关，这称为单摆的等时性。

摆锤的角速度为ω = d θ/dt ，因此
可得
积分得
22d sin d g t l θθ=-22d sin d ml mg t θθ=-0ω=002π
2T ω
==2202d 0d t θωθ+=21cos 2g C l
ωθ=+mg d sin d g l ωωθθ=-22d d d d d d d d d d t t t θωωθωωθθ
===
当t = 0时，ω = 0，θ = θm ，可得C = -gcos θm/l 。

因此角速度大小为
单摆的周期为：
对于任何角振幅θm ，通过数值积分和符号积分都能计算周期。

利用半角公式可得：
设 2sin
k m θ= 并设ksinx = sin(θ/2)，因此
可得
第一类完全椭圆积分定义为
则周期为
其中
d d t θω==m m 000
πT T ==
⎰m 001πT T θ=
⎰002π2T ω=
=π/2π/2000012ππT T T ==⎰⎰π/2002πT =⎰π/2
0K()k =⎰02K()π.T T k =002π2T ω==π/2002πT T =⎰02K(),πT k =m sin 2
k θ=
将周期的椭圆积分公式按二项式展开得
其中(2n – 1)!! = 1·3…·(2n – 1)。

利用定积分公式
可得用无穷级数表示的单摆周期
如果只取常数项，可得单摆小角度的周期T0。

如果取前两个正弦项，则得
利用级数计算周期究竟要取多少项，则根据精度确定。

三、算法及程序代码
1.单摆的周期与角振幅关系
【算法】MATLAB 定义的第一类完全椭圆积分为
周期可表示为 π/2
20102(21)!![1(sin )]d π2!n n n n T T k x x n ∞=-=+∑⎰π/220(21)!!πsin d 2!2n n n x x n -=⎰2m 01(21)!!{1[sin ]}2!2
n n n n T T n θ∞=-=+∑224m m 0222113(1sin sin ...)22242
T T θθ=++
+π/20K(m)=⎰02K()π.T T m =
其中2
sin m 2m θ=. 对于一定的角振幅，根据精度在单摆周期的无穷级数中去有限项，计算周期的近似值。

【程序代码】
clear
theta=0:179;
thm=theta*pi/180;
T0=ellipke(sin(thm/2).^2)*2/pi; figure
plot(theta,T0,'r-')
fs=16;
title('单摆的周期与角振幅的关系','FontSize',fs)
xlabel('\it\theta\rm_m/\circ','FontSize',fs)
ylabel('\itT/T\rm_0','FontSize',fs)
grid on
legend('第一类完全椭圆积分')
text(0,2,'\itT\rm_0=2\pi(\itl/g\rm)^{1/2}','FontSize',fs)
2.轻杆单摆运动的动画模拟
【程序代码】
clear
thetam=input('输入单摆角振幅：');
thm=thetam*pi/180;
T=ellipke(sin(thm/2)^2)*2/pi;
t=linspace(0,T*2*pi);
options.RelTol=1e-6;
[t,TH]=ode45('fun',t,[thm,0],options);
th=TH(:,1);
figure
plot([0;0],[0;1.05],'-.','LineWidth',2)
axis equal off ij
fs=16;
title('单摆的振动','FontSize',fs)
hold on
plot(0,0,'o')
plot(exp(i*(linspace(-thm,thm)+pi/2)),'m--','LineWidth',2) x=sin(thm);
y=cos(thm);
plot([0,-x],[0,y],'r:','LineWidth',2)
pole=plot([0,x],[0,y],'r','LineWidth',3);
ball=plot(x,y,'c.','MarkerSize',50);
txt{1}=['\it\theta\rm_m=',num2str(thetam),'\circ'];
txt{2}=['\itT/T\rm_0=',num2str(T)];
txt{3}='\itT\rm_0=2\pi(\itl/g\rm)^{1/2}';
text(-x/2,y/2,txt,'FontSize',fs)
pause
while get(gcf,'CurrentCharacter')~=char(27)
for i=1:length(t)
x=sin(th(i));
y=cos(th(i));
set(ball,'XData',x,'YData',y)
set(pole,'XData',[0 x],'YData',[0 y])
drawnow
pause(0.01)
end
end
3.轻杆单摆运动规律的模拟
【程序代码】
clear
w0t=linspace(0,4*2*pi);
a=1:5;
n=length(a);
options.RelTol=1e-6;W=[]; THETA=[];
t0=[];
for i=1:n
th0=a(i)*pi/180;
[w0t,TH]=ode45('p5_5_2fun',w0t,[th0;0],options); THETA=[THETA,TH(:,1)*180/pi];
W=[W,TH(:,2)];
t0=[t0,ellipke(sin(th0/2)^2)*2/pi];
end
t=w0t/2/pi;
figure
plot(t,THETA(:,1),'o-',t,THETA(:,2),'d-',t,THETA(:,3),'s-', ...t,THETA(:,4),'^-',t,THETA(:,5),'v-')
grid on
fs=16;
xlabel('时间\itt/T\rm_0','FontSize',fs)
ylabel('角度\it\theta\rm/\circ','FontSize',fs)
title('单摆的角位移(同色点为简谐振动的标准点)','FontSize',fs)
leg=[repmat('\it\theta\rm_m=',n,1),num2str(real(a'))];%角振leg=[leg,repmat('\circ,\itT/T\rm_0=',n,1),num2str(real(t0') )];
legend(leg,4)
text(0,a(end),'\itT\rm_0=2\pi(\itl/g\rm)^{1/2}','FontSize',
fs)
[T,WT]=meshgrid(t0,w0t);
THm=ones(size(t))*a;
THh=THm.*cos(WT./T);
hold on
plot(t,THh,'.')
figure
plot(t,W)
plot(t,W(:,1),'o-',t,W(:,2),'d-',t,W(:,3),'s-',t,W(:,4),'^-',...t,W(:,5),'v-')
grid on
title('单摆的角速度(同色点为简谐振动的标准点)','FontSize',fs)
xlabel('时间\itt/T\rm_0','FontSize',fs)
ylabel('角速度\it\omega/\omega\rm_0','FontSize',fs)
legend(leg,4)
text(0,max(W(:)),'\it\omega\rm_0=2\pi/\itT\rm_0','FontSize' ,fs)
Wm=ones(size(w0t))*a*pi/180;
Wh=-Wm./T.*sin(WT./T);
hold on
plot(t,Wh,'.')
四.执行结果及分析
1.单摆的周期与角振幅的关系
程序执行结果如图2.a和图2.b所示。

图2. a
数值积分和符号积分与第一类完全椭圆积分公式计算的结果完全吻合。

当角振幅在20º以内时，单摆的周期几乎不变；当角振幅在20º到40º之间时，单摆的周期稍有增加；当角振幅大于40º时，单摆的周期显著增加；当角振幅接近180º时，单摆的周期急剧增加。

在图2.b中，实线表示用第一类完全椭圆积分公式计算的周期的精确值，点表示用级数计算的周期的近似值，容差取6
10 .当角振幅等于5º时，只要在周期的级数中取一个正弦项即可达到精度。

当角振幅等于150º时，则需要取148个正弦项才能达到精度。

当角振幅在
155º到165º之间时，取150个正弦项虽然不能达到精度，但是周期的近似值与精确值基本吻合。

当角振幅等于90º时，则需要取15个正弦项才能达到精度。

可见：在通常振幅的情况下，可用级数求和的方法计算单摆的周期，但是在很大振幅的情况下，就需要用积分的方法或完全椭圆积分函数才能保证周期的精度。

图2. b
2.轻杆单摆的动画演示效果
为了演示单摆的振动，需要求微分方程中角度的数值解。

摆锤的坐标为x = lsinθ，y = lcosθ，x轴取向右的方向为正，y轴取向下的方向为正。

（1）当角振幅为60º时，单摆的初始状态如图3所示，单摆的周期为 1.0732T0。

图3
（2）当角振幅小于5º时，单摆振动周期约等于小角振动的周期，如图4。

图4
（3）当角振幅为90º时，单摆的周期为1.18T0
，如图5。

图5
（4）当角振幅为179º时，单摆的初始状态如图6所示，单摆的周期
为3.90T0。

图6
3. 轻杆单摆运动规律及角位移曲线
（1）在角振幅较小的情况下，单摆的周期近似为小角单摆的周期，其角位移完全可以用余弦函数，如图7所示。

图7
（2）在角振幅较小的情况下，其角速度完全可以用正弦函数表示，如图8所示。

图8
（3）当角振幅在90º以内时，单摆的角位移与简谐运动的标准点基本上是重合的，因此可用余弦函数近似表示；当角振幅等于150º时，单摆的角位移与简谐运动的标准点有所偏离，如图9所示。

图9
（4）当角振幅在90º以内时，单摆的角速度曲线与大多数正弦点(少量极值附近的点除外)重合；当角振幅等于150º时，单摆的角速度与正弦曲线偏离较多；当角振幅接近180º时，角速度与正弦曲线偏离很大，峰值附近的曲线尖而窄，零值附近的曲线变得十分“平直”。

图10
五.参考资料网站
/view/60720.htm
/view/cd18810a4a7302768e9939f3.html /question/237542980.html
/HotWord/4251328.htm
附录：组员成绩评定
童鑫：负责程序和算法编写，对研究工作做出了突破性的贡献。

最后评定结果：8分
周涛：负责文档整理与校正，资料搜集等工作。

最后评定分数：8分
周慧斌：负责问题分析和数学模型的建立，公式解析等工作。