牛头刨床运动分析 程序一： clear;clc; %初始条件 L1=125;L6=275;L3=600;L4=150;L6p=575;w1=1;theta1=0:1:360; %theta3 方位角 theta3=atan((L6+L1.*sind(theta1))./(L1.*cosd(theta1))); theta3=theta3.*180/pi; for i=1:361 if theta3(i)<0 theta3(i)=90+abs(abs(theta3(i))-90); end end %S3 长度 S3=(L1.*cosd(theta1))./cosd(theta3); %theta4 方位角 theta4=asin((L6p-L3.*sind(theta3))./L4); theta4=theta4.*180./pi; for i=1:361 if theta4(i)<0 theta4(i)=180+abs(theta4(i)); B=[-w1.*L1.*sind(theta1(i)); w1.*L1.*cosd(theta1(i)); 0; 0]; C=inv(A)*B; vS3(i)=C(1); w3(i)=C(2); w4(i)=C(3); ve(i)=C(4); end %% %aS3,a3,a4,aSe 的求解 for i=1:361 else theta4(i)=180-theta4(i); end end %Se 长度 Se=L3.*cosd(theta3)+L4.*cosd(theta4); %% %vS3,w3,w4,vSe 的求解 for i=1:361 A=[cosd(theta3(i)), -S3(i).*sind(theta3(i)), 0, 0; sind(theta3(i)), S3(i).*cosd(theta3(i)), 0, 0; 0 ,-L3.*sind(theta3(i)) ,-L4.*sind(theta4(i)), -1; 0, L3.*cosd(theta3(i)) ,L4.*cosd(theta4(i)), 0]; 

D=[-w3(i).*sind(theta3(i)),-vS3(i).*sind(theta3(i))-S3(i).*w3(i).*cos d(theta3(i)),0,0; w3(i).*cosd(theta3(i)),vS3(i).*cosd(theta3(i))-S3(i).*w3(i).*sind(the ta3(i)),0,0; 0,-L3.*w3(i).*cosd(theta3(i)),-L4.*w4(i).*cosd(theta4(i)),0; 0,-L3.*w3(i).*sind(theta3(i)),-L4.*w4(i).*sind(theta4(i)),0;]; E=[vS3(i);w3(i);w4(i);ve(i)]; F=w1.*[-w1.*L1.*cosd(theta1(i)); -w1.*L1.*sind(theta1(i)); 0; 0]; G=inv(A)*(-D*E+F); aS3(i)=G(1); a3(i)=G(2); a4(i)=G(3); ae(i)=G(4); end %% figure(1) plot(theta1,theta3,'g');hold on; h1=plotyy(theta1,theta4,theta1,Se);grid on; xlabel('theta1'); ylabel(h1(1),'theta3、theta4/度'); ylabel(h1(2),'Se/mm'); legend('theta3','theta4','Se') title('位置线图'); figure(2) plot(theta1,w3,'g');hold on; h2=plotyy(theta1,w4,theta1,ve);grid on; xlabel('theta1'); ylabel(h2(1),'w3、w4/(rad/s)'); ylabel(h2(2),'ve/(m/s)'); legend('w3','w4','ve') title('速度线图'); figure(3) plot(theta1,a3,'g');hold on; h3=plotyy(theta1,a4,theta1,ae);grid on; xlabel('theta1'); ylabel(h3(1),'a3、a4/(rad/s^2)'); ylabel(h3(2),'ae/(m/s^2)'); legend('a3','a4','ae') title('加速度线图'); 

程序二： clear all;clc; w1=1;l1=0.125;l3=0.6;l6=0.275;l61=0.575;l4=0.15; for m=1:3601 o1(m)=pi*(m-1)/1800; o31(m)=atan((l6+l1*sin(o1(m)))/(l1*cos(o1(m)))); if o31(m)>=0 o3(m)=o31(m); else o3(m)=pi+o31(m); end; s3(m)=(l1*cos(o1(m)))/cos(o3(m)); o4(m)=pi-asin((l61-l3*sin(o3(m)))/l4); se(m)=l3*cos(o3(m))+l4*cos(o4(m)); if o1(m)==pi/2 o3(m)=pi/2; s3(m)=l1+l6; end if o1(m)==3*pi/2 o3(m)=pi/2; s3(m)=l6-l1; end A1=[cos(o3(m)),-s3(m)*sin(o3(m)),0,0;sin(o3(m)),s3(m)*cos(o3(m)),0,0; 0,-l3*sin(o3(m)),-l4*sin(o4(m)),-1;0,l3*cos(o3(m)),l4*cos(o4(m)),0]; B1=w1*[-l1*sin(o1(m));l1*cos(o1(m));0;0];D1=A1\B1;E1(:,m)=D1;ds(m)=D1 (1);w3(m)=D1(2);w4(m)=D1(3);ve(m)=D1(4); A2=[cos(o3(m)),-s3(m)*sin(o3(m)),0,0;sin(o3(m)),s3(m)*cos(o3(m)),0,0; 0,-l3*sin(o3(m)),-l4*sin(o4(m)),-1;0,l3*cos(o3(m)),l4*cos(o4(m)),0]; B2=-[-w3(m)*sin(o3(m)),(-ds(m)*sin(o3(m))-s3(m)*w3(m)*cos(o3(m))),0,0; w3(m)*cos(o3(m)),(ds(m)*cos(o3(m))-s3(m)*w3(m)*sin(o3(m))),0,0;0,-l3* w3(m)*cos(o3(m)),-l4*w4(m)*cos(o4(m)),0;0,-l3*w3(m)*sin(o3(m)),-l4*w4 (m)*sin(o4(m)),0]*[ds(m);w3(m);w4(m);ve(m)]; C2=w1*[-l1*w1*cos(o1(m));-l1*w1*sin(o1(m));0;0];B=B2+C2;D2=A2\B;E2(:, m)=D2;dds(m)=D2(1);a3(m)=D2(2);a4(m)=D2(3);ae(m)=D2(4); end; % o11=o1*180/pi;o33=o3*180/pi; % y=polyfit(o11,o33,4); for x=1:3601 x1=pi*(x-1)/1800; y1(x)=polyval(y,x1); end % o11=o1*180/pi;y=[o3*180/pi;o4*180/pi];w=[w3;w4];a=[a3;a4];figure; % subplot(221);h1=plotyy(o11,y,o11, se); axis equal; % title('位置线图 

');xlabel('\it\theta1');ylabel('\it\theta3,\theta4,Se'); % subplot(222);h2=plotyy(o11,w,o11,ve); % title('速度线图'); % xlabel('\it\theta1');ylabel('\it\omega3,\omega4,Ve'); % subplot(212);h3=plotyy(o11,a,o11,ae); % title('加速度线图'); % xlabel('\it\theta1');ylabel('\it\alpha3,\alpha4,\alphaE'); % F=[o11;o3./pi*180;o4./pi*180;se;w3;w4;ve;a3;a4;ae]';G=F(1:100:3601,:) ;