MATLAB机构运动仿真

机械原理第一次作业（matlab7.0）：求：r1旋转360°时，θ2，θ3，ω2，ω3，α2，α3和C点的加速度.设r1=1000,r2=2000,r3=1800,r4=25001、角位移的M函数：functiony=diyiwen(x)%Inputparameters%x(1)=theta-1%x(2)=theta-2guessvalue%x(3)=theta-3guessvalue%x(4)=r1%x(5)=r2%x(6)=r3%x(7)=r4%Outputparameters%y(1)=theta-2%y(2)=theta-3theta2=x(2);theta3=x(3);%epsilon=1.0E-6;%f=[x(4)*cos(x(1))+x(5)*cos(theta2)-x(7)-x(6)*cos(theta3);x(4)*sin(x(1))+x(5)*sin(theta2)-x(6)*sin(theta3)];%whilenorm(f)>epsilonJ=[-x(5)*sin(theta2)x(6)*sin(theta3);x(5)*cos(theta2)-x(6)*cos(theta3)];dth=inv(J)*(-1.0*f);theta2=theta2+dth(1);theta3=theta3+dth(2);f=[x(4)*cos(x(1))+x(5)*cos(theta2)-x(7)-x(6)*cos(theta3);x(4)*sin(x(1))+x(5)*sin(theta2)-x(6)*sin(theta3)];norm(f);end;y(1)=theta2;y(2)=theta3;r1旋转360°时，θ2，θ3的M文件程序：r(1)=1000;r(2)=2000;r(3)=1800;r(4)=2500;dr=pi/180;th(1)=0;th(2)=59.8897*dr;th(3)=106.0172*dr;y=diyiwen([th(1),th(2),th(3),r(1),r(2),r(3),r(4)])dth=1*dr;fori=1:360y=diyiwen([th(1),th(2),th(3),r(1),r(2),r(3),r(4)]);th23(i,:)=[th(1)/dr,th(2)/dr,th(3)/dr];th(1)=th(1)+dth;th(2)=y(1);th(3)=y(2);endfigure(1)plot(th23(:,1),th23(:,2),th23(:,1),th23(:,3),th23(:,1),th23(:,1))axis([0,360,0,360])gridontitle('角位移线图')xlabel('曲柄转角th(1)')ylabel('从动件转角th(2),th(3)')text(300,50,'角th(2)')text(300,150,'角th(3)')text(200,200,'角th(1)')2、角速度的M函数：functiony=dierwen(x)%Iputparameters%%x(1)=theta-1%x(2)=theta-2%x(3)=theta-3%x(4)=dtheta-1%x(5)=r1%x(6)=r2%x(7)=r3%%Outputparameters%%y(1)=dtheta-2%y(2)=dtheta-3%A=[-x(6)*sin(x(2))x(7)*sin(x(3));x(6)*cos(x(2))-x(7)*cos(x(3))];B=[x(5)*sin(x(1));-x(5)*cos(x(1))]*x(4);y=inv(A)*B;ω2，ω3r1旋转360°时ω2，ω3的M文件程序：th(1)=0;dr=pi/180;dth(1)=10;dth(2)=59.8897*dr;dth(3)=106.0172*dr;dth1=1*dr;r(1)=1000;r(2)=2000;r(3)=1800;r(4)=2500;forj=1:360x=dierwen([th(1),dth(2),dth(3),dth(1),r(1),r(2),r(3)]);y=diyiwen([th(1),dth(2),dth(3),r(1),r(2),r(3),r(4)]);dth(2)=y(1);dth(3)=y(2);th(1)=th(1)+dth1;dth23(j,:)=[th(1)/dr,x(1)/dr,x(2)/dr];endfigure(2)plot(dth23(:,1),dth23(:,2),dth23(:,1),dth23(:,3))axis([0,360,-400,400])gridontitle('角速度线图')xlabel('曲柄旋转角度')ylabel('从动件角速度dth(2),dth(3)(rad/s)')text(250,200,'连杆2角速度')text(50,100,'摇杆3角速度')3、r1旋转360°时，α2，α3程序M文件：r(1)=1000;r(2)=2000;r(3)=1800;r(4)=2500;th(1)=0;dr=pi/180;dth(2)=59.8897*dr;dth(3)=106.0172*dr;dth(1)=10;dth2=1*dr;fork=1:360x=dierwen([th(1),dth(2),dth(3),dth(1),r(1),r(2),r(3)]);y=diyiwen([th(1),dth(2),dth(3),r(1),r(2),r(3),r(4)]);dth(2)=y(1);dth(3)=y(2);w1=x(1);w2=x(2);C=[-r(2)*sin(dth(2)),r(3)*sin(dth(3));r(2)*cos(dth(2)),-r(3)*cos(dth(3))];D=-[-w1*r(2)*cos(dth(2)),w2*r(3)*cos(dth(3));-w2*r(2)*sin(dth(2)),w2*r(3)*sin(dth(3))]…*[w1;w2]+10*[10*r(1)*cos(th(1));10*r(1)*sin(th(1))];a=inv(C)*D;a2(k)=a(1);a3(k)=a(2);jiao(k)=th(1);th(1)=th(1)+dth2;endfigure(3)plot(jiao/dr,a2,jiao/dr,a3)axis([0,360,-100,100])gridontitle('角加速度线图')xlabel('曲柄旋转转角')ylabel('从动件角加速度（rad/s^2）')text(200,40,'连杆2的角加速度')text(150,-40,'摇杆3的角加速度')4、C点的加速度函数：C点的M函数文件：functiony=cdjsd(x)%x(1)=r(1);%x(2)=r(2);%x(3)=th(1);%x(4)=th(2);%x(5)=w1;%x(6)=w2;%x(7)=a(1);y=[-x(1)*sin(x(3)),-x(2)*sin(x(4));x(1)*cos(x(3)),x(2)*cos(x(4))]*[0;x(7)]+[-x(1)*x(5)*cos(x(3)),-x(2)*x(6)*cos(x(4));-x(1)*x(5)*sin(x(3)),-x(2)*x(6)*sin(x(3))]*[x(5)^2;x(6)^2];r1旋转360°时，C点的加速度的M文件程序：r(1)=1000;r(2)=2000;r(3)=1800;r(4)=2500;th(1)=0;dr=pi/180;dth(2)=59.8897*dr;dth(3)=106.0172*dr;dth(1)=10;dth2=1*dr;fork=1:360x=dierwen([th(1),dth(2),dth(3),dth(1),r(1),r(2),r(3)]);y=diyiwen([th(1),dth(2),dth(3),r(1),r(2),r(3),r(4)]);dth(2)=y(1);dth(3)=y(2);w1=x(1);w2=x(2);C=[-r(2)*sin(dth(2)),r(3)*sin(dth(3));r(2)*cos(dth(2)),-r(3)*cos(dth(3))];D=-[-w1*r(2)*cos(dth(2)),w2*r(3)*cos(dth(3));-w2*r(2)*sin(dth(2)),w2*r(3)*sin(dth(3))]*[w1;w2]+10*[10*r(1)*cos(th(1));10*r(1)*sin(th(1))];a=inv(C)*D;a2=a(1);ac=cdjsd([r(1),r(2),th(1),dth(2),w1,w2,a2]);th(1)=th(1)+dth2;acx(k)=