单摆在MATLAB中的应用

《MATLAB语言》课程论文





运用MATLAB语言解决单摆的问题






姓名：李毛毛
学号:
专业：电子信息工程
班级：2010级
指导老师：汤全武
学院：物理电气信息学院
完成日期：20111225

运用MATLAB语言解决单摆的相关问题
（李毛毛 电子信息工程 2010级1班）
【摘要】
：
借助MATLAB计算软件, 制作一个动态单摆的模 型，使人能直观地看到单摆的
运动情况；运用MATLAB语言对单摆水平分力和竖直分力进行研究，数 据可以直观的看出
力的变化范围；最后还运用MATLAB语言对单摆的角速度和单摆的二阶性进行了分 析。建
立单摆摆角的二阶微分方程模型，并化为单摆的摆角和角速度的一阶微分方程组，在相平面
上分析单摆在圆周上的运动时，摆角和角速度的关系。
【关 键 词】
:单摆 周期 摆长 单摆的二阶性 角速度
MATLAB语言

单摆是一种理想的物理模型，它由理想化的摆球和摆线 组成．摆线由
质量不计、不可伸缩的细线提供；摆球密度较大，而且球的半径比摆线的
长度小得 多，这样才可以将摆球看做质点，由摆线和摆球构成单摆．在满
足偏角
?
<10°的条 件下，单摆的周期


t?2
?
l
g
（1） < br>单摆问题是高中物理及大学普通物理实验教学中的一个基础问题。单
摆在摆角比较小时，其运动规 律近似为准简谐振动。但是当摆角比较大时，
即单摆在大摆角情况下运动时，这种近似已不在成立，其运 动方程满足非
线性微分方程。因此，对摆角大小的限制成为该实验中必须满足的条件。
不同的实 验条件下，最大摆角的取值不同，其中包括
3?,5?,10?
。因此，这就
为在实验 过程中对摆角的统一取值造成困难，给实验带来较大的误差。同
时，学生对单摆在大摆角情况下运动时其 运动周期及运动规律的理解也存
在困难。MATLAB是用于算法开发，数据可视化、数据分析以及数值 计算的
高级技术计算语言。本文将借助于MATLAB这个强大的数学软件，计算了任
意摆角下 单摆运动周期的精确解，消除了摆角问题带来的误差。同时利用
该软件，仿真出大摆角时单摆的运动情况 ，为单摆测重力加速度实验中大
摆角问题的讲解提供了较好的手段。
一、单摆的介绍

1.单摆的定义
用一根绝对挠性且长度不变、质量可忽略不计的线悬挂一个质点，在重力作用下在铅垂平面内作周期运动，就成为单摆。

2

2．摆长的定义

摆长指的是从悬点到摆球重心的距离。
3．平衡位置
当摆球静止在平衡位置O点时，细线竖直下垂，摆球所受重力G和悬 线的拉
力F平衡，O点就是摆球的平衡位置。
4．回复力
将摆球拉离O 点，从A位置释放，摆球在细线的拉力F和重力G的作用下,
沿AOA'圆弧在平衡点O附近来回振动， 当摆球运动到任一位置P时，此时摆球
受重力G，拉力F作用，由于摆球沿圆弧运动，所以将重力分解成 切线方向分力
G1；和沿半径方向G2，重力分力G1不论是在O左侧还是右侧始终指向平衡位置，而且正是在G1作用下摆球才能回到平衡位置
5．条件
摆角θ≤10°
6．单摆二阶性的介绍

(1) 建立单摆摆角的二阶微分方程模型
设单摆摆线与铅垂线夹角为
?
(t)
，则摆锤切向速度为
v?l
d< br>?
dt
，把重力
mg
分解成
MP
与
MQ两个分量
由牛顿第二定律知：
m
dv
??mgsin
?
dt
（2）
d
2
?
?ksin
?
?0
2
:即 ：
dt
（3）
k?
g
l
，这是单摆摆角的二阶非线性微分方程模型. 其中
(2) 建立单摆的摆角和角速度的一阶微分方程组模型
上述二阶方程模型等价于下面的一阶微分方程组模型：
?
d
?
?
?
?
?
dt
?
?
d
?
??ksin
?
?
?
dt


（4）
3

（

其中
?
为角速度。
由此可以研究摆角和角速度的相位关系，从而观察单摆的运动情况。
第二步：摆角和角速度的相位关系及其相位图：
d
?
ksin
?
??
d
??
(5)

积分之，得摆角和角速度的相位关系：
1
2
?
?k(1?cos
?
)?E
2

（6）
E
是非负常数, 这里未考虑空气阻力。
上（6）式表达了单摆的 能量守恒律，左端第一项是动能项，第二项是势能项，
E
是总能量（不考虑摆锤质量、摆长等常 数因子）。
?
??2[E?k(1?cos
?
)]
（7）
?
号表示摆线摆动方向是逆时针或顺时针方向。
总能量
E
是单摆开始摆动时给定的。令
E
0
?2k
，这是单摆的最大势能。由
上述公式我们知道，角速度关于摆角是周期变化的，周期为
2
?
，即单摆的动 能
随着摆角的变化而周期变化的。
二、利用MATLAB处理单摆的相关问题
1．运用Matlab语言制作单摆
%挂摆横梁
%by dynamic
%see also http:
%2008.6.9
h=figure('numbertitle','off',' name','擦除动画演示(挂摆横梁)——
Matlabsky')
%绘制横梁
plot([-0.2;0.2],[0;0],'-k','linewidth',20); %画初始位置的单摆
g=0.98; %重力加速度,可以调节摆的摆速
l=1; %摆长
theta0=pi4; %初始角度
x0=l*sin(theta0);%初始x坐标
y0=-l*cos(theta0);%初始y坐标
axis([-0.75,0.75,-1.25,0]);

4

axis off %创建摆锤
%擦除模式为xor

head=line(x0,y0,'color','r','linestyle','.','e rasemode','xor','ma
rkersize',40); %创建摆杆

body=line([0;x0],[-0.05;y0],'color','b','lines tyle','-','erasemod
e','xor'); %摆的运动
t=0;%时间变量
dt=0.01;%时间增量
while 1
t=t+dt;
theta=theta0*cos(sqrt(gl)*t);%单摆角度与时间的关系
x=l*sin(theta); y=-l*cos(theta);
if ~ishandle(h),return,end
set(head,'xdata',x,'ydata',y);%改变擦除对象的坐标数据
set(body,'xdata',[0;x],'ydata',[-0.05;y]);
drawnow;%刷新屏幕

运行结果如图2所示 :

图2单摆的模型
h=

5

1
2. 研究单摆周期与摆长之间的关系
程序如下：
l=1:0.1:2;
%摆长的变化范围

g=9.8;
%当地的重力加速度

t=2*pi*sqrt(lg)
%周期的计算公式

运行结果如下：
t =
Columns 1 through 7
2.0071 2.1051 2.1987 2.2884 2.3748 2.4582 2.5388
Columns 8 through 11
2.6169 2.6928 2.7666 2.8385
3. 对单摆水平和垂直方向的分力的研究
假设单摆的各个量如下：初始角速 度300度秒，摆与水平夹角30度，杆长0.45
米,小球质量2kg。
程序如下：
theta0=-0.524;%最初的角度
w=-5.24;% 角速度
m=2.0; %摆球质量
L=0.45; %摆长
g=9.8; %当地重力加速度
t=0:0.01:2 %时间的变化范围
theta=theta0+w*t; %t时刻所对应的角度
fx=m*w*w*L*cos(theta) %垂直方向的分量
fy=m*w*w*L*sin(theta)+m*g %水平方向的分量
运行结果如下：
fx =
Columns 1 through 7
21.3961 20.7192 19.9853 19.1966 18.3552 17.4634 16.5237
Columns 8 through 14
15.5386 14.5108 13.4432 12.3387 11.2004 10.0313 8.8346
Columns 15 through 21
7.6137 6.3719 5.1127 3.8393 2.5555 1.2646 -0.0297
Columns 22 through 28
-1.3240 -2.6146 -3.8981 -5.1708 -6.4294 -7.6703 -8.8902
Columns 29 through 35
-10.0856 -11.2534 -12.3903 -13.4931 -14.5589 -15.5848 -16.5678
Columns 36 through 42
-17.5054 -18.3950 -19.2340 -20.0202 -20.7515 -21.4258 -22.0413
Columns 43 through 49
-22.5963 -23.0893 -23.5188 -23.8839 -24.1833 -24.4164 -24.5824
Columns 50 through 56

6

-24.6810 -24.7118 -24.6747 -24.5700 -24.3977 -24.1586 -23.8530
Columns 57 through 63
-23.4821 -23.0466 -22.5479 -21.9873 -21.3663 -20.6867 -19.9503
Columns 64 through 70
19.1591 -18.3153 -17.4212 -16.4794 -15.4923 -14.4626 -13.3933
Columns 71 through 77
-12.2872 -11.1473 -9.9769 -8.7791 -7.5571 -6.3145 -5.0544
Columns 78 through 84
-3.7806 -2.4963 -1.2052 0.0892 1.3834 2.6738 3.9568
Columns 85 through 91
5.2290 6.4868 7.7268 8.9457 10.1399 11.3063 12.4417
Columns 92 through 98
13.5429 14.6070 15.6309 16.6119 17.5474 18.4346 19.2713
Columns 99 through 105
20.0551 20.7838 21.4554 22.0682 22.6203 23.1104 23.5370
Columns 106 through 112
23.8991 24.1955 24.4255 24.5884 24.6839 24.7116 24.6714
Columns 113 through 119
24.5635 24.3882 24.1460 23.8374 23.4635 23.0251 22.5235
Columns 120 through 126
21.9600 21.3363 20.6541 19.9151 19.1214 18.2753 17.3790
Columns 127 through 133
16.4350 15.4459 14.4143 13.3432 12.2355 11.0942 9.9224
Columns 134 through 140
8.7234 7.5005 6.2569 4.9962 3.7218 2.4371 1.1458
Columns 141 through 147
-0.1487 -1.4428 -2.7329 -4.0155 -5.2871 -6.5442 -7.7833
Columns 148 through 154
-9.0011 -10.1941 -11.3592 -12.4931 -13.5926 -14.6549 -15.6769
Columns 155 through 161
-16.6559 -17.5892 -18.4742 -19.3085 -20.0898 -20.8159 -21.4849
Columns 162 through 168
-22.0949 -22.6442 -23.1314 -23.5551 -23.9141 -24.2075 -24.4344
Columns 169 through 175
-24.5943 -24.6866 -24.7112 -24.6679 -24.5569 -24.3786 -24.1332
Columns 176 through 182
-23.8217 -23.4447 -23.0034 -22.4989 -21.9327 -21.3063 -20.6213
Columns 183 through 189
-19.8798 -19.0837 -18.2352 -17.3367 -16.3905 -15.3994 -14.3660
Columns 190 through 196
-13.2931 -12.1838 -11.0410 -9.8679 -8.6677 -7.4438 -6.1993
Columns 197 through 201
-4.9379 -3.6629 -2.3779 -1.0863 0.2082
fy =

7

Columns 1 through 7
7.2355 6.1318 5.0651 4.0383 3.0542 2.1156 1.2249
Columns 8 through 14
0.3847 -0.4028 -1.1353 -1.8110 -2.4279 -2.9843 -3.4787
Columns 15 through 21
-3.9097 -4.2762 -4.5772 -4.8118 -4.9794 -5.0795 -5.1118
Columns 22 through 28
-5.0763 -4.9731 -4.8025 -4.5648 -4.2608 -3.8913 -3.4573
Columns 29 through 35
-2.9600 -2.4008 -1.7812 -1.1029 -0.3678 0.4222 1.2647
Columns 36 through 42
2.1577 3.0985 4.0846 5.1133 6.1817 7.2870 8.4261
Columns 43 through 49
9.5959 10.7931 12.0146 13.2568 14.5164 15.7900 17.0741
Columns 50 through 56
18.3651 19.6595 20.9537 22.2442 23.5275 24.7999 26.0581
Columns 57 through 63
27.2986 28.5179 29.7128 30.8799 32.0160 33.1180 34.1829
Columns 64 through 70
35.2078 36.1899 37.1264 38.0148 38.8527 39.6377 40.3676
Columns 71 through 77
41.0406 41.6548 42.2083 42.6999 43.1280 43.4915 43.7894
Columns 78 through 84
44.0209 44.1854 44.2824 44.3117 44.2731 44.1668 43.9930
Columns 85 through 91
43.7523 43.4453 43.0728 42.6358 42.1357 41.5737 40.9513
Columns 92 through 98
40.2704 39.5327 38.7403 37.8953 37.0001 36.0572 35.0691
Columns 99 through 105
34.0385 32.9683 31.8613 30.7208 29.5497 28.3512 27.1288
Columns 106 through 112
25.8857 24.6253 23.3512 22.0667 20.7755 19.4810 18.1869
Columns 113 through 119
16.8966 15.6138 14.3419 13.0845 11.8449 10.6266 9.4330
Columns 120 through 126
8.2672 7.1326 6.0322 4.9691 3.9461 2.9661 2.0317
Columns 127 through 133
1.1456 0.3101 -0.4724 -1.1998 -1.8701 -2.4815 -3.0323
Columns 134 through 140
-3.5209 -3.9461 -4.3066 -4.6015 -4.8300 -4.9914 -5.0853
Columns 141 through 147
-5.1114 -5.0697 -4.9603 -4.7834 -4.5396 -4.2296 -3.8541
Columns 148 through 154
-3.4142 -2.9112 -2.3464 -1.7213 -1.0377 -0.2975 0.4974

8

Columns 155 through 161
1.3447 2.2422 3.1872 4.1774 5.2098 6.2818 7.3903
Columns 162 through 168
8.5324 9.7048 10.9044 12.1279 13.3719 14.6329 15.9076
Columns 169 through 175
17.1925 18.4839 19.7785 21.0725 22.3625 23.6449 24.9162
Columns 176 through 182
26.1729 27.4116 28.6288 29.8212 30.9856 32.1187 33.2175
Columns 183 through 189
34.2788 35.2999 36.2779 37.2101 38.0939 38.9270 39.7071
Columns 190 through 196
40.4319 41.0995 41.7082 42.2561 42.7419 43.1641 43.5216
Columns 197 through 201
43.8135 44.0389 44.1972 44.2880 44.3110
4.讨论角速度与摆角之间的关系
方程

dy
1
dy?y
2
,
2
?q?y
2
?sin
?
y
1
?
?b?cos
?
w
0
?t
?
,其中
y
1
是摆角，
y
2
是角速
dtdt
度。
b
的范围是[1.5,2]。
编的程序如下:
q = 0.5; %对q赋值
w0 = 2.03; %角速度
bdata = []; %b赋值
y2data = []; %角速度赋初值
for b = 1.5:0.01:2 %b的变化范围
fun = @(t,y)[y(2);-q*y(2)-sin(y(1))+b*cos(w0*t)]; %对角速度微分
ts = linspace(0,5,100); %对ts赋值
[t,y] = ode45(fun,ts,[0,2]); % 在此区域生成网格坐标
y1 = y(:,1); %摆角
y2 = y(:,2); %角速度
y2i = y2(abs(y1-1.59)<0.048); %b的约束条件
bdata = [bdata;b*ones(size(y2i))];
y2data = [y2data;y2i];
end
plot(bdata,y2data,'*-') %绘制曲线
xlabel('b')
ylabel('角速度')
运行结果如图3所示:


9

图3 角速度与常数b的关系
5.单摆的二阶性
程序如下：
%单摆系统，要考虑空气阻力，以及出示状态，经过数学分析，该系统是一个二阶系统
%什么是二阶系统，就是有两个状态变量（可以进行微分的）
g=9.8;%当地重力加速度
T=15;%该单摆的周期
dt=0.0001;%时间的变化间隔
t=0:dt:T;%给时间t赋值
L=1;%摆长
m=10;%摆锤
K=5;
v0=0;%初速度赋值
theta0=3.1;
v=zeros(size(t));%存储变量初始化
theta=zeros(size(t));%变量初始化
v(1)=v0;
theta(1)=theta0;%赋初值
for k=1:length(t)%开始仿真 二阶方程组 描述系统运行
v(k+1)=v(k)+(g*sin(theta(k))-K.m.*v(k)).*dt;
theta(k+1)=theta(k)-1.L.*v(k).*dt;
end
[AX ,H1,H2]=plotyy(t,v(1:length(t)),t,theta(1:length(t )),'plot');
set(H1,'LineStyle','-');%设置作图类型，前者实线后者虚线
set(H2,'LineStyle','-.'); %设置作图类型，前者实线后者虚线
set(get(AX(1),'Ylabel'),'string','线速度 v(t) ms');
%设置y轴的标注
set(get(AX(2),'Ylabel'),'string','角位移 theta(t) rads');

10

xlabel('时间 ts');
legend(H1,'线速度 v(t)',1);%说明曲线1名称
legend(H2,'角位移 theta(t)',2); %说明曲线2名称
结果如图4所示：

图4单摆的二阶性
三、结论
1. 利用Matlab语言制作单摆，可以看到一个动态的单摆运动，使人更加清
晰，便捷的看到一幅动态的 画面，便于理解。
2.从运行结果可以看出：随着摆长的增加周期也随之增大。即：摆长与周期
成正比例关系。
3.单摆的水平分量和竖直分量随着摆角的变化而变化，当单摆角为0度时，
竖直分量为0；当摆角为9 0度时，水平分量为0。
4.角速度与摆角有一定的关系，当给定一个摆角时，就会有一个角 速度与之
对应。结果应该是角速度随b线性增加。
5.运用Matlab语言研究单 摆的二阶性，就是对两个状态变量进行微分的，
可以使繁琐的运算量变得简洁明了，运行结果清晰可见。
四、课程体会

本课程开始讲的是一些matlab的入门技术，除了一些基本 操作与介绍之外，
还初步认识了简单指令的编制，认识了一些matlab的特殊符号，例如运算用到< br>的加减乘除。（+—*）对我来说还是比较容易上手的。
其次，我们接触到两种数据类 型，一个是double，还有就是char,另外还介
绍了赋值语句，数学计算，常用函数，输入输出 语句和数据文件。
最后，我们接触到基本的matlab选择结构，还有控制这个结构的关系 运算
符和逻辑运算符。也就是if结构。这个结构对运算起到很大作用，跟elseif

11

配合使用的话，就可以将很图方法。还有，我们学习如何控制 画图附加功能，例
如线的宽度和符号的颜色。这些属性可由指定的“propertyname”和值V alue
决定，“propertyname”和值Value将出现在plot命令的数据后。

在这个学期的学习中， 学了一些初步的矩阵运算以及画图方法。在书上的
练习中也得到了 实践，各种矩阵的运算也可以运用到现在所学的电路计算中。活
学活用，的确为我们提供了不少的方便。 但现今所学到的matlab技术还只是很
初步，但是我已经了解到了matlab的实用性，所以还是 会在今后继续学习这门
课程的。
除此之外还学了利用axis命令的画从这个学期开始，我们 学习了一门新课
程，开始我对matlab这门课程不是很了解，只知道它是一种计算机程序，可以帮助我们运算处理一些数据，但并不知道这种程序与其他计算机编程程序有什么
区别与特点。 运用Matlab语言可以使繁琐的绘图变得简单，绘图本来是一件繁琐的工作，
需要确定很多参数 ，但是MATLAB高层绘图函数代替用户做了这些工作，给定了
参数的默认值，这样就是使用户省去了 一些操作细节，用起来比较方便。
通过学习matlab这门课程后，我知道，matlab是一种专 业的计算机程序，
主要用于工程科学的矩阵数学运算。但在之后的几年里，它又渐渐发展成为一种
非常灵活的计算体系，可以解决更多技术上的问题。在解决工程技术问题方面，
matlab比其他任 何计算机语言都简单高效，对于我们学电气信息工程专业的学
生，这门课程给我们日后本专业的工作提供 了很大的方便。
参考文献：
1. 刘卫国. MATLAB程序设计与应用[M].北京：高等教育出版社，2002.
2. 籍延坤;焦志伟. 单摆作简谐振动的最大摆角的分析 [期刊论文] -大学
物理实验2001(04).
3. 叶慧群．单摆周期近似解法综述EJ]．浙江师范大学学报：自然科学版，
2004，29(3)：34 —38．
4. 刘卫国.科学计算与MATLAB语言[M].北京：中国铁道出版社，2000.
5. 张志涌.精通MATLAB6.5[M].北京：北京航空航天大学出版社，2003.
6. 易德文，盛忠志．利用单摆测重力加速度时的最大摆角的再讨论[J]．物
理实验，20 03；23(7)：44—45．
7. ]高闯，魏薇，肖发新．单摆角度问题分析EJ]．高等函授 学报，2004，
17(1)：35—36．
8. 贾利群．物理实验学[M]_北京：东方出版社，1999：103．

12