matlab仿真弹簧波,[转载]matlab小算法——弹簧的动态仿真（二维）

%%http://control.blog.sina.com.cn/admin/article/article_add.php

theta0=pi/3;

m=1;k=80;g=9.8;

L0=1;L=L0+m*g/k;

[t,u1]=ode45('thbfun',[0:0.005:15],[L0 0 theta0

0],[],L,k,m,g);

[y1,x1]=pol2cart(u1(:,3),u1(:,1));

y1=-y1;

figure

ymax=max(abs(y1));

axis([-1.2 1.2 -1.2*ymax 0.2]);

axis off

title('弹簧摆','fontsize',14);

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);

line([-1 1],[0 0],'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.0 0.51

0.6],'linestyle','-','linewidth',1.3,'erasemode','none');

spring=line(xx,yy,'color','g','linewidth',2,'erasemode','xor');

pause(0.5)

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

end

function F=thbfun(t,u,flag,l,k,m,g)

F=[u(2);

u(1)*u(4)^2+g*cos(u(3))-k/m*(u(1)-1+m*g/k);

u(4);

-2*u(2)*u(4)/u(1)-g*sin(u(3))/u(1)];