用MATLAB模拟弹簧摆运动

用MATLAB模拟弹簧摆运动引言：物体在一定位置附近所作的往复运动称为机械振动。这种振动现象在自然界是广泛存在的。例如，摆的运动，一切发声体的运动，机器开动时各部分的微小颤动等都是机械振动。在不同的振动现象中，最基本最简单的振动是简谐振动。一切复杂的振动都可以分解为若干个简谐振动，弹簧振子和小角度单摆就是简单的简谐振动。计算机模拟可以直观表现物理系统的运动状态，对于理解抽象的物理公式和研究系统的运动带来很多方便，所以在科研中有广泛的应用。这里我们借助于MATLAB,在研究物理系统的运动的同时，也用动画模拟他们的运动。由于这些运动模拟式采用真实的从求解运动方程得到的数据。所以利用matlab的数据拟合可以画出精确的图形。题目：设质量为m的摆锤挂在劲度系数为k,原长为l0的轻弹簧上，弹簧的另一端悬挂于固定点O,系统静止自然下垂时弹簧长度为0mgllk，系统可在过O点的竖直平面内自由摆动，试研究摆锤的运动。原理：系统自由度为2。以O为极点，竖直向下的Ox轴为极轴，建立极坐标系，如图所示，r为质点m到O点距离，为Ox轴与弹簧间的夹角，则系统的拉格朗日函数为由拉格朗日方程可求出系统的运动微分方程为上式未作小摆角近似，因此可用以研究弹簧摆的大摆角运动，但是难以求出解析解。由数值计算结果所画出弹簧摆的大摆角运动轨迹，可发现它的运动情况很复杂，不作数值计算是无法想象的。令1yr，2drydt，3y，4dydt，则成为源程序：functionthbglobalLmkgtheta0=pi/10;%初始角度,可由读者设不同的值m=1;k=80;g=9.8;L0=1;L=L0+m*g/k;%L0为弹簧原来长度，L为弹簧静止时长度[t,u1]=ode45(@thbfun,[0:0.005:15],[L00theta00]);[y1,x1]=pol2cart(u1(:,3),u1(:,1));y1=-y1;%将极坐标换为直角坐标figureymax=max(abs(y1));axis([-1.21.2-1.2*ymax0.2]);%设置坐标范围axisofftitle('弹簧摆','fontsize',14)holdon;R=0.055;%设置弹簧半径yy=-L0:0.01:0;xx=R*sin(yy./L0*30*pi);%用正弦曲线表示弹簧[a,r]=cart2pol(xx,yy);%用坐标变换来画初始位置的弹簧a=a+theta0;[xx,yy]=pol2cart(a,r);%弹簧的数据line([-11],[00],'color','r','linewidth',2)%画横杆，球与弹簧ball=line(xx(1),yy(1),'color','r','marker','.','markersize',70,'erasemode','xor');ball2=line(xx(1),yy(1),'color',[0.50.510.6],'linestyle','-','linewidth',1.3,'erasemode','none');spring=line(xx,yy,'color','r','linewidth',2,'erasemode','xor');pause(0.5)fori=1:length(t)%画出每一步的弹簧位置yy=-u1(i,1):0.01:0;xx=R*sin(yy./u1(i,1)*30*pi);[a,r]=cart2pol(xx,yy);a=a+u1(i,3);[xx,yy]=pol2cart(a,r);set(ball,'XData',x1(i),'YData',y1(i));set(ball2,'XData',x1(i),'YData',y1(i));set(spring,'XData',xx,'YData',yy);drawnow;endfunctionF=thbfun(t,u)globalLmkgF=[u(2);u(1)*u(4)^2+g*cos(u(3))-k/m*(u(1)-L+m*g/k);u(4);-2*u(2)*u(4)/u(1)-g*sin(u(3))/u(1)];程序说明：程序中解微分方程的过程比较简单，主要的技巧是在做模拟动画上。有两点值得注意，一是模拟弹簧的运动，另一点是在模拟弹簧运动的同时画出轨迹图。主程序的文件名是thb，子程序文件为thbfun，两个程序共用的变量L，m，k，g用指令global来传递。程序运行截图：总结：这次小论文的编写，虽然程序还是相对简单的，但是体现了MATLAB真正用于科学计算的优势。通过这门课的学习，对MATLAB也有了更多的了解，以后可以更好的运用MATLAB，辅助各种研究工作，受益匪浅。参考文献：彭芳麟.计算物理基础.北京：高等教育出版社，2010.1