弹簧摆matlab演示及源代码

弹簧摆写出弹簧摆的运动微分方程，并用matlab编写程序，演示弹簧摆的运动弹簧摆的运动微分方程可以通过拉格朗日函数法求得，对于弹簧单摆系统，在极坐标系中其拉格朗日函数为222211()cos()22LTVmgmrrmgrkrlk其中,r为自由度，,,,mkgl分别为摆球质量，弹簧刚度，重力加速度和摆初始长。将拉格朗日函数带入到拉格朗日方程00dLLdtrrdLLdt得到系统的运动微分方程：sin()02sin0kmgrrgrlgkrrg对于这个非线性方程组，利用matlab里面自带的ode45常微分方程组求解函数，对上面的方程组进行数值求解，得到了摆球的运动轨迹，如下图所示0246810121416-1.5-1-0.500.511.522.533.54弹簧摆内共振动画模拟及自由度随时间的变化曲线时间t单位秒对应的自由度摆角随时间变化曲线摆长随时间变化曲线图中给出了弹簧摆运动过程中自由度,r随时间的变化关系，下面将给出这个程序的源代码，演示程序时，首先停顿三秒，然后开始摆动，图形上方的弹簧摆的示意图开始摆动并描绘出轨迹，同时图形下方自由度随时间变化关系也会动态的给出，从这个过程中可以明显看出弹簧摆的内共振现象。程序源码：functionspringmasstheta0=1.5*pi/10;%单摆的初始角度m=20/9.8;k=80;g=9.8;%弹簧质量的参数保证弹簧固有频率是单摆固有频率的2倍L0=1;L=L0+m*g/k;%L0为弹簧原来长度，L为弹簧静止时长度[t,u1]=ode45(@weifen,[0:0.05:15],[L00theta00],[],L,k,m,g);[y1,x1]=pol2cart(u1(:,3),u1(:,1));y1=-y1;%将极坐标换为直角坐标axis([016-1.54]);%text(0,u1(1,1),'摆长l');%text(0,u1(1,3),'摆角\theta');%axisofftitle('弹簧摆内共振动画模拟及自由度随时间的变化曲线','fontsize',15)xlabel('时间t单位秒');ylabel('对应的自由度');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);%弹簧的数据xx=2*xx+7;yy=yy+3.5;%xx(1)=x(1)+15;line([68],[3.53.5],'color','k','linewidth',4)%弹簧单摆系统初始状态fori=1:9line([6+0.2*i6+0.2*i+0.2],[3.53.7],'color','k','linewidth',1);endline([77],[3.52.3],'color','k','linestyle','-.','linewidth',1)bigball=line(xx(1),yy(1),'color','b','marker','.',...'markersize',30,'erasemode','xor');%球ball2=line(xx(1),yy(1),'color','g','linestyle','-',...'linewidth',1.3,'erasemode','none');%轨线spring=line(xx,yy,'color','k','linewidth',2,'erasemode','xor');%弹簧linex=line(t(1),u1(1,1),'color','r','marker','.',...'markersize',10,'erasemode','none');%摆长曲线初始liney=line(t(1),u1(1,3),'color','b','marker','.',...'markersize',10,'erasemode','none');%摆角曲线初始%以下为标注text(6,-0.75,'\uparrow');text(4,-1,'摆角随时间变化曲线');text(9.0,1.9,'摆长随时间变化曲线');text(11.4,1.70,'\downarrow');%在初始状态下停顿三秒pause(2)%以下程序开始实现动画，在弹簧摆开始振动时画出相应的变化曲线fori=2:length(t)yy=-u1(i,1):0.01:0;%弹簧xx=R*sin(yy./u1(i,1)*30*pi);%xx(i)=x(i)+15;[a,r]=cart2pol(xx,yy);a=a+u1(i,3);%set(linex,'XData',t(i),'YData',u1(i,1));%%摆长初始%set(ball2,'XData',2*x1(i)+7,'YData',y1(i)+3.5);%轨线%set(liney,'XData',t(i),'YData',u1(i,3));%摆角初始plot([t(i-1),t(i)],[u1((i-1),1),u1(i,1)],'xr');plot([t(i-1),t(i)],[u1((i-1),3),u1(i,3)],'*b');%legend('摆长曲线','摆角曲线')[xx,yy]=pol2cart(a,r);xx=2*xx+7;yy=yy+3.5;set(bigball,'XData',2*x1(i)+7,'YData',y1(i)+3.5);%球%plot(2*x1(i)+7,y1(i)+3.5,'b');%球%set(ball2,'XData',2*x1(i)+7,'YData',y1(i)+3.5);%轨线plot([2*x1(i-1)+7,2*x1(i)+7],[y1(i-1)+3.5,y1(i)+3.5],'-g');%轨线set(spring,'XData',xx,'YData',yy);%弹簧drawnow;end%以下为弹簧摆的运动微分方程functionvalue=weifen(t,u,l,k,m,g)value=[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)];