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: 1530497909','MATLAB答疑','help')
web http://url.cn/NSFcAs -browser
