搜索
Matlab（程序双摆运动的研究的）双摆运动有下面的微分方程组来描述x1’=x3;x2’=x4x3’=A(1)x4’=A(2)其中，A表示一个列向量，由式（10-1）计算。(m1+m2)L1m2L2cos(x1-x2)-1m2L2x4*x4sin(x2-x1)-（m1+m2）gsinx1A=M1L1cos（x1-x2）m2L2m2L1x3*x3sin（x1-x2）-m2gsinx2其中，m1和m2分别表示上摆和下摆的质量，L1和L2分别表示上摆和下摆的长度，g表示重力加速度，x1和x2分别表示上摆和下摆的角位移，x3和x4分别表示上摆和下摆的角速度。对上面的微分方程组可以使用ode45函数来求解。其对应的MATLAB代码为clearall;m1=1;m2=1;L1=1;L2=1;g=9.8;Da=inline(['[x(3);x(4);',...'inv([(m1+m2)*L1,m2*L2*cos(x(1)-x(2));',...'m1*L1*cos(x(1)-x(2)),m1*L2])*'...'[m2*L2*x(4)^2*sin(x(2)-x(1))-(m1+m2)*g*sin(x(1));',...'m2*L1*x(3)^2*sin(x(1)-x(2))-m2*g*sin(x(2))]]'],'t','x',...'flag','m1','m2','L1','L2','g');set(gcf,'DoubleBuffer','on');[t,x]=ode45(Da,[0,20];[0.8,0.8,0,0],[],m1,m2,L1,L2,g);axis([-(L1+L2),(L1+L2),-(L1+L2)*1.8,1]);axissquare;holdon;gh1=plot([0,L1*exp(i*(x(1)-pi/2))],'r-');set(gh1,'linewidth',2,'markersize',6,'marker','o');gh2=plot([L1*exp(i*(x(1)-pi/2)),L1*exp(i*(x(1)-pi/2))+L2*exp(i*(x(2)-pi/2))],'b-');set(gh2,'linewidth',2,'markersize',6,'marker','o');fork=2:size(x,1);C1=[0,L1*exp(i*(x(k,1)-pi/2))];C2=[L1*exp(i*(x(k,1)-pi/2)),L1*exp(i*(x(k,1)-pi/2))+L2*exp(i*(x(k,2)-pi/2))];set(gh1,'xdata',real(C1),'ydata',imag(C1));set(gh2,'xdata',real(C2),'ydata',imag(C2));title(['t=',num2str(t(k))],'fontsize',12);pause(0.1);endfigure;subplot(2,3,1);plot(t,x(:,1));title('t-\theta_1');xlabel('t');ylabel('\theta_1');subplot(2,3,2);plot(t,x(:,2));title('t-\theta_2');xlabel('t');ylabel('\theta_2');subplot(2,3,3);plot(t,x(:,3));title('t-\omega_1');xlabel('t');ylabel('\omega_1');subplot(2,3,4);plot(t,x(:,4));title('t-\omega_2');xlabel('t');ylabel('\omega_2');subplot(2,3,5);plot(x(:,1),x(:,3));title('\theta_1-\omega_1');xlabel('\theta_1');ylabel('\omega_1');subplot(2,3,6);plot(x(:,2),x(:,4));title('\theta_2-\omega_2');