基于Matlab模拟单摆在外力矩作用下的动画

clear;clc;close all l=4; m=1; g=9.81; zeta=0.2; I=m*l^2; w0=sqrt(g/l); c=zeta*2*w0*I; dt=0.2/w0; theta(1)=0.6; dtheta(1)=0; Ek=0.5*m*(l*dtheta(1))^2; Ep=l*(1-cos(theta(1)))*m*g; figure h1=subplot(2,1,1); set(h1,'XLim',[-8 8],'YLim',[-5 0]) set(gca,'drawmode','fast') line_handle1=plot([0 l*sin(theta(1))],[0 -l*cos(theta(1))],'-o'); axis equal title('Pendulum Simulation') blank=[' ']; parameters={['Length' blank(1:14) num2str(l)];['Mass' blank(1:16) num2str(m)]; ... ['Damping ratio' blank(1:4) num2str(zeta)];['Gravity' blank(1:14) num2str(g)]}; text(3,-1.5,parameters) text_show={['Time' blank(1:10) '0'];['Theta' blank(1:9) num2str(theta(1))]; ... ['dTheta/dt' blank(1:4) num2str(dtheta(1))];['Energy' blank(1:7) num2str(Ek+Ep)]}; text_handle=text(-7.5,-3.5,text_show); h2=subplot(2,1,2); set(h2,'XLim',[-0.8 0.8],'YLim',[-0.8 0.8]) set(gca,'drawmode','fast') line_handle2=plot(theta(1),dtheta(1)); axis equal title('theta-dtheta curve') xlabel('theta (rad)') ylabel('dtheta (rad/s)') N=80; for k=1:N Tw=5*sin(2*(k-1)*dt); ddtheta(k)=(Tw/I)-(g/l)*sin(theta(k))-(c/(m*l^2))*dtheta(k); dtheta(k+1)=dtheta(k)+dt*ddtheta(k); theta(k+1)=theta(k)+dt*dtheta(k); Ek=0.5*m*(l*dtheta(k+1))^2; Ep=l*(1-cos(theta(k+1)))*m*g; set(line_handle1,'xdata',[0 l*sin(theta(k+1))],'ydata',[0 -l*cos(theta(k+1))]) set(h1,'XLim',[-8 8],'YLim',[-5 0]) text_show={['Time' blank(1:10) num2str(k*dt)];['Theta' blank(1:9) num2str(theta(k+1))]; ... ['dTheta/dt' blank(1:4) num2str(dtheta(k+1))];['Energy' blank(1:7) num2str(Ek+Ep)]}; set(text_handle,'string',text_show) set(line_handle2,'xdata',theta,'ydata',dtheta) set(h2,'XLim',[-0.8 0.8],'YLim',[-0.8 0.8]) drawnow pause(0.001) end msgbox('MATLAB编程答疑，请加QQ: 1575304183','MATLAB答疑','help') web http://url.cn/NSFcAs -browser