dynamics.m

matlab中用于解决高等动力学hamilton系统的相关程序

L=0.15;
m=0.2;
I=1/12*m*L^2;
C=[3 L 3/2*L L/2;m*L (I+m*L^2/2) 3/4*m*L^2 1/4*m*L^2;3/2*m*L 1/2*m*L^2 (I+5/4*m*L^2) 1/2*m*L^2; m*L/2 1/4*m*L^2 1/2*m*L^2 (I+1/4*m*L^2)];
Impulse=[1;1*L/2;0;0];
A=inv(C)*Impulse;
[t,u] = ode45(@pulse,[0 10],[-0.15;0;0;0;0]);
figure(1)
for i=1:90
    clf    
    axis([-0.6 0.6 -0.6 0.6]);
    xlabel('X (m)');
    ylabel('Y (m)');
    linewidth=1000;
    line([u(i,1)+L/2*cos(u(i,3));u(i,1)-L/2*cos(u(i,3))],[u(i,2)+L/2*sin(u(i,3));u(i,2)-L/2*sin(u(i,3))],'LineWidth',3);
    line([u(i,1)+L/2*cos(u(i,3));u(i,1)+L/2*cos(u(i,3))+L*cos(u(i,4))],[u(i,2)+L/2*sin(u(i,3));u(i,2)+L/2*sin(u(i,3))+L*sin(u(i,4))],'LineWidth',3);
    line([u(i,1)+L/2*cos(u(i,3))+L*cos(u(i,4));u(i,1)+L/2*cos(u(i,3))+L*cos(u(i,4))+L*cos(u(i,5))],[u(i,2)+L/2*sin(u(i,3))+L*sin(u(i,4));u(i,2)+L/2*sin(u(i,3))+L*sin(u(i,4))+L*sin(u(i,5))],'LineWidth',3);
    pause(0.01)
    MOV(i)=getframe;
end