弹簧摆matlab演示及源代码

弹簧摆

写出弹簧摆的运动微分方程，并用matlab编写程序，演示弹簧摆的运动 弹簧摆的运动微分方程可以通过拉格朗日函数法求得，对于弹簧单摆系统，在极坐标系中其拉格朗日函数为

LTV 11mg2m(r2r22)mgrcosk(rl)22k

其中r,为自由度，m,k,g,l分别为摆球质量，弹簧刚度，重力加速度和摆初始长。

将拉格朗日函数带入到拉格朗日方程

ddtddtLL0rrLL0

得到系统的运动微分方程：

kmgrrgsin(rl)0gkr2rgsin0

对于这个非线性方程组，利用matlab里面自带的ode45常微分方程组求解函数，对上面的方程组进行数值求解，得到了摆球的运动轨迹，如下图所示

43.532.52弹簧摆内共振动画模拟及自由度随时间的变化曲线对应的自由度摆长随时间变化曲线1.510.50-0.5-1-1.502摆角随时间变化曲线46810时间t 单位秒121416

图中给出了弹簧摆运动过程中自由度r,随时间的变化关系，下面将给出这个程序的源代码，演示程序时，首先停顿三秒，然后开始摆动，图形上方的弹簧摆的示意图开始摆动并描绘出轨迹，同时图形下方自由度随时间变化关系也会动态的给出，从这个过程中可以明显看出弹簧摆的内共振现象。

程序源码：

function springmass

theta0=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],[L0 0 theta0 0],[],L,k,m,g);

[y1,x1] = pol2cart(u1(:,3),u1(:,1)); y1 = -y1;%将极坐标换为直角坐标

axis([0 16 -1.5 4]);

% text(0,u1(1,1),'摆长l');

% text(0,u1(1,3),'摆角\heta');

%axis off

title('弹簧摆内共振动画模拟及自由度随时间的变化曲线','fontsize',15)

xlabel('时间t 单位秒');

ylabel('对应的自由度');

hold on;

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([6 8],[3.5 3.5],'color','k','linewidth',4)%弹簧单摆系统初始状态

for i=1:9

line([6+0.2*i 6+0.2*i+0.2],[3.5 3.7],'color','k','linewidth',1);

end

line([7 7],[3.5 2.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,'\ow');

text(4,-1,'摆角随时间变化曲线');

text(9.0,1.9,'摆长随时间变化曲线');

text(11.4,1.70,'\\downarrow');

%在初始状态下停顿三秒

pause(2)

%以下程序开始实现动画，在弹簧摆开始振动时画出相应的变化曲线

for i = 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

%以下为弹簧摆的运动微分方程

function value=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)];