单摆在MATLAB中的应用

MATLAB语言课程论文运用MATLAB语言解决单摆的问题姓名：李毛毛学号:12010245339 专业：电子信息工程班级：2010级指导老师：汤全武 学院：物理电气信息学院完成日期：2011/12/25运用MATLAB语言解决单摆的相关问题（李毛毛 电子信息工程 2010级1班）【摘要】：借助MATLAB计算软件, 制作一个动态单摆的模型，使人能直观地看到单摆的运动情况；运用MATLAB语言对单摆水平分力和竖直分力进行研究，数据可以直观的看出力的变化范围；最后还运用MATLAB语言对单摆的角速度和单摆的二阶性进行了分析。建立单摆摆角的二阶微分方程模型，并化为单摆的摆角和角速度的一阶微分方程组，在相平面上分析单摆在圆周上的运动时，摆角和角速度的关系。【关 键 词】:单摆 周期 摆长 单摆的二阶性 角速度 MATLAB语言 单摆是一种理想的物理模型，它由理想化的摆球和摆线组成摆线由质量不计、不可伸缩的细线提供；摆球密度较大，而且球的半径比摆线的长度小得多，这样才可以将摆球看做质点，由摆线和摆球构成单摆在满足偏角10的条件下，单摆的周期（1）单摆问题是高中物理及大学普通物理实验教学中的一个基础问题。单摆在摆角比较小时，其运动规律近似为准简谐振动。但是当摆角比较大时，即单摆在大摆角情况下运动时，这种近似已不在成立，其运动方程满足非线性微分方程。因此，对摆角大小的限制成为该实验中必须满足的条件。不同的实验条件下，最大摆角的取值不同，其中包括。因此，这就为在实验过程中对摆角的统一取值造成困难，给实验带来较大的误差。同时，学生对单摆在大摆角情况下运动时其运动周期及运动规律的理解也存在困难。MATLAB是用于算法开发，数据可视化、数据分析以及数值计算的高级技术计算语言。本文将借助于MATLAB这个强大的数学软件，计算了任意摆角下单摆运动周期的精确解，消除了摆角问题带来的误差。同时利用该软件，仿真出大摆角时单摆的运动情况，为单摆测重力加速度实验中大摆角问题的讲解提供了较好的手段。一、单摆的介绍1.单摆的定义用一根绝对挠性且长度不变、质量可忽略不计的线悬挂一个质点，在重力作用下在铅垂平面内作周期运动，就成为单摆。2摆长的定义 摆长指的是从悬点到摆球重心的距离。3平衡位置 当摆球静止在平衡位置O点时，细线竖直下垂，摆球所受重力G和悬线的拉力F平衡，O点就是摆球的平衡位置。4回复力 将摆球拉离O点，从A位置释放，摆球在细线的拉力F和重力G的作用下,沿AOA圆弧在平衡点O附近来回振动，当摆球运动到任一位置P时，此时摆球受重力G，拉力F作用，由于摆球沿圆弧运动，所以将重力分解成切线方向分力G1；和沿半径方向G2，重力分力G1不论是在O左侧还是右侧始终指向平衡位置，而且正是在G1作用下摆球才能回到平衡位置5条件 摆角106单摆二阶性的介绍(1) 建立单摆摆角的二阶微分方程模型 设单摆摆线与铅垂线夹角为，则摆锤切向速度为，把重力mg分解成与两个分量 由牛顿第二定律知： （2）:即 ： （3）其中，这是单摆摆角的二阶非线性微分方程模型.(2) 建立单摆的摆角和角速度的一阶微分方程组模型上述二阶方程模型等价于下面的一阶微分方程组模型： （4） （其中为角速度。由此可以研究摆角和角速度的相位关系，从而观察单摆的运动情况。第二步：摆角和角速度的相位关系及其相位图： (5)积分之，得摆角和角速度的相位关系： （6）是非负常数, 这里未考虑空气阻力。上（6）式表达了单摆的能量守恒律，左端第一项是动能项，第二项是势能项，是总能量（不考虑摆锤质量、摆长等常数因子）。 （7）号表示摆线摆动方向是逆时针或顺时针方向。 总能量是单摆开始摆动时给定的。令，这是单摆的最大势能。由上述公式我们知道，角速度关于摆角是周期变化的，周期为，即单摆的动能随着摆角的变化而周期变化的。二、利用MATLAB处理单摆的相关问题1运用Matlab语言制作单摆 %挂摆横梁 %by dynamic %see also %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=pi/4; %初始角度 x0=l*sin(theta0);%初始x坐标 y0=-l*cos(theta0);%初始y坐标 axis(-0.75,0.75,-1.25,0); axis off %创建摆锤 %擦除模式为xor head=line(x0,y0,color,r,linestyle,.,erasemode,xor,markersize,40); %创建摆杆 body=line(0;x0,-0.05;y0,color,b,linestyle,-,erasemode,xor); %摆的运动 t=0;%时间变量 dt=0.01;%时间增量 while 1 t=t+dt; theta=theta0*cos(sqrt(g/l)*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= 12. 研究单摆周期与摆长之间的关系程序如下：l=1:0.1:2; %摆长的变化范围g=9.8; %当地的重力加速度t=2*pi*sqrt(l/g) %周期的计算公式运行结果如下：t =Columns 1 through 7 2.0071 2.1051 2.1987 2.2884 2.3748 2.4582 2.5388Columns 8 through 112.6169 2.6928 2.7666 2.83853. 对单摆水平和垂直方向的分力的研究假设单摆的各个量如下：初始角速度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 721.3961 20.7192 19.9853 19.1966 18.3552 17.4634 16.5237 Columns 8 through 1415.5386 14.5108 13.4432 12.3387 11.2004 10.0313 8.8346 Columns 15 through 217.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-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 7019.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 915.2290 6.4868 7.7268 8.9457 10.1399 11.3063 12.4417 Columns 92 through 9813.5429 14.6070 15.6309 16.6119 17.5474 18.4346 19.2713 Columns 99 through 10520.0551 20.7838 21.4554 22.0682 22.6203 23.1104 23.5370 Columns 106 through 11223.8991 24.1955 24.4255 24.5884 24.6839 24.7116 24.6714 Columns 113 through 11924.5635 24.3882 24.1460 23.8374 23.4635 23.0251 22.5235 Columns 120 through 12621.9600 21.3363 20.6541 19.9151 19.1214 18.2753 17.3790 Columns 127 through 13316.4350 15.4459 14.4143 13.3432 12.2355 11.0942 9.9224 Columns 134 through 1408.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.2082fy = Columns 1 through 77.2355 6.1318 5.0651 4.0383 3.0542 2.1156 1.2249 Columns 8 through 140.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 422.1577 3.0985 4.0846 5.1133 6.1817 7.2870 8.4261 Columns 43 through 499.5959 10.7931 12.0146 13.2568 14.5164 15.7900 17.0741 Columns 50 through 5618.3651 19.6595 20.9537 22.2442 23.5275 24.7999 26.0581 Columns 57 through 6327.2986 28.5179 29.7128 30.8799 32.0160 33.1180 34.1829 Columns 64 through 7035.2078 36.1899 37.1264 38.0148 38.8527 39.6377 40.3676 Columns 71 through 7741.0406 41.6548 42.2083 42.6999 43.1280 43.4915 43.7894 Columns 78 through 8444.0209 44.1854 44.2824 44.3117 44.2731 44.1668 43.9930 Columns 85 through 9143.7523 43.4453 43.0728 42.6358 42.1357 41.5737 40.9513 Columns 92 through 9840.2704 39.5327 38.7403 37.8953 37.0001 36.0572 35.0691 Columns 99 through 10534.0385 32.9683 31.8613 30.7208 29.5497 28.3512 27.1288 Columns 106 through 11225.8857 24.6253 23.3512 22.0667 20.7755 19.4810 18.1869 Columns 113 through 11916.8966 15.6138 14.3419 13.0845 11.8449 10.6266 9.4330 Columns 120 through 1268.2672 7.1326 6.0322 4.9691 3.9461 2.9661 2.0317 Columns 127 through 1331.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 Columns 155 through 1611.3447 2.2422 3.1872 4.1774 5.2098 6.2818 7.3903 Columns 162 through 1688.5324 9.7048 10.9044 12.1279 13.3719 14.6329 15.9076 Columns 169 through 17517.1925 18.4839 19.7785 21.0725 22.3625 23.6449 24.9162 Columns 176 through 18226.1729 27.4116 28.6288 29.8212 30.9856 32.1187 33.2175 Columns 183 through 18934.2788 35.2999 36.2779 37.2101 38.0939 38.9270 39.7071 Columns 190 through 19640.4319 41.0995 41.7082 42.2561 42.7419 43.1641 43.5216 Columns 197 through 20143.8135 44.0389 44.1972 44.2880 44.31104.讨论角速度与摆角之间的关系方程,其中是摆角，是角速度。的范围是1.5,2。编的程序如下: q = 0.5; %对q赋值w0 = 2.0/3; %角速度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;endplot(bdata,y2data,*-) %绘制曲线xlabel(b)ylabel(角速度)运行结果如图3所示:图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;endAX,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) m/s);%设置y轴的标注set(get(AX(2),Ylabel),string,角位移 theta(t) rad/s);xlabel(时间 t/s);legend(H1,线速度 v(t),1);%说明曲线1名称legend(H2,角位移 theta(t),2); %说明曲线2名称结果如图4所示：图4单摆的二阶性三、结论 1.利用Matlab语言制作单摆，可以看到一个动态的单摆运动，使人更加清晰，便捷的看到一幅动态的画面，便于理解。 2.从运行结果可以看出：随着摆长的增加周期也随之增大。即：摆长与周期成正比例关系。 3.单摆的水平分量和竖直分量随着摆角的变化而变化，当单摆角为0度时，竖直分量为0；当摆角为90度时，水平分量为0。 4.角速度与摆角有一定的关系，当给定一个摆角时，就会有一个角速度与之对应。结果应该是角速度随b线性增加。 5.运用Matlab语言研究单摆的二阶性，就是对两个状态变量进行微分的，可以使繁琐的运算量变得简洁明了，运行结果清晰可见。四、课程体会 本课程开始讲的是一些matlab的入门技术，除了一些基本操作与介绍之外，还初步认识了简单指令的编制，认识了一些matlab的特殊符号，例如运算用到的加减乘除。（+*/）对我来说还是比较容易上手的。 其次，我们接触到两种数据类型，一个是double，还有就是char,另外还介绍了赋值语句，数学计算，常用函数，输入输出语句和数据文件。 最后，我们接触到基本的matlab选择结构，还有控制这个结构的关系运算符和逻辑运算符。也就是if结构。这个结构对运算起到很大作用，跟elseif配合使用的话，就可以将很图方法。还有，我们学习如何控制画图附加功能，例如线的宽度和符号的颜色。这些属性可由指定的“propertyname”和值Value决定，“propertyname”和值Value将出现在plot命令的数据后。在这个学期的学习中， 学了一些初步的矩阵运算以及画图方法。在书上的练习中也得到了实践，各种矩阵的运算也可以运用到现在所学的电路计算中。活学活用，的确为我们提供了不少的方便。但现今所学到的matlab技术还只是很初步，但是我已经了解到了matlab的实用性，所以还是会在今后继续学习这门课程的。除此之外还学了利用axis命令的画从这个学期开始，我们学习了一门新课程，开始我对matlab这门课程不是很了解，只知道它是一种计算机程序，可以帮助我们运算处理一些数据，但并不知道这种程序与其他计算机编程程序有什么区别与特点。运用Matlab语言可以使繁琐的绘图变得简单，绘图本来是一件繁琐的工作，需要确定很多参数，但是MATLAB高层绘图函数代替用户做了这些工作，给定了参数的默认值，这样就是使用户省去了一些操作细节，用起来比较方便。通过学习matlab这门课程后，我知道，matlab是一种专业的计算机程序，主要用于工程科学的矩阵数学运算。但在之后的几年里，它又渐渐发展成为一种非常灵活的计算体系，可以解决更多技术上的问题。在解决工程技术问题方面，matlab比其他任何计算机语言都简单高效，对于我们学电气信息工程专业的学生，这门课程给我们日后本专业的工作提供了很大的方便。参考文献：1. 刘卫国. MATLAB程序设计与应用M.北京：高等教育出版社，2002.2. 籍延坤;焦志伟. 单摆作简谐振动的最大摆角的分析 期刊论文 -大学物理实验2001(04).3. 叶慧群单摆周期近似解法综述EJ浙江师范大学学报：自然科学版，2004，29(3)：34384. 刘卫国.科学计算与MATLAB语言M.北京：中国铁道出版社，2000.5. 张志涌.精通MATLAB6.5M.北京：北京航空航天大学出版社，2003.6. 易德文，盛忠志利用单摆测重力加速度时的最大摆角的再讨论J物理实验，2003；23(7)：44457. 高闯，魏薇，肖发新单摆角度问题分析EJ高等函授学报，2004，17(1)：35368. 贾利群物理实验学M_北京：东方出版社，1999：103