matlab四杆机构源程序

1.角位移M函数：
function y=jweiyi(x)
% Input parameters
% x(1)=theta-1
% x(2)=theta-2 guess value
% x(3)=theta-3 guess value
% x(4)=r1
% x(5)=r2
% x(6)=r3
% x(7)=r4
% Output parameters
% y(1)=theta-2
% y(2)=theta-3
theta2=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)];
%
while norm(f)>epsilon
J=[-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;
2. r1旋转360°时，θ2，θ3的M文件程序：
r(1)=50;
r(2)=150;
r(3)=300;
r(4)=400;
dr=pi/180;
th(1)=0;
th(2)=58.4119*dr;
th(3)=154.7912*dr;
y=jweiyi([th(1),th(2),th(3),r(1),r(2),r(3),r(4)])
dth=1*dr;
for i=1:360
y=jweiyi([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);
end
figure(1)
plot(th23(:,1),th23(:,2),th23(:,1),th23(:,3),th23(:,1),th23(:,1)) axis([0,360,0,360])
grid on
title('角位移线图')
xlabel('曲柄转角th(1)')
ylabel('从动件转角th(2),th(3)')
text(300,50,'角th(2)')
text(300,150,'角th(3)')
text(200,200,'角th(1)')
3.角速度的M函数：
function y=rrrvel(x)
% Input parameters
%
% 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
%
% Output parameters
%
% 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;
4. r1旋转360°时ω2，ω3的M文件程序：
th(1)=0;
dr=pi/180;
dth(1)=36;
dth(2)=58.4119*dr;
dth(3)=154.7912*dr;
dth1=1*dr;
r(1)=50;
r(2)=150;
r(3)=300;
r(4)=400;
for j=1:360
x=rrrvel([th(1),dth(2),dth(3),dth(1),r(1),r(2),r(3)]);
y=jweiyi([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];
end
figure(2)
plot(dth23(:,1),dth23(:,2),dth23(:,1),dth23(:,3))
axis([0,360,-800,2000])
grid on
title('角速度线图')
xlabel('曲柄旋转角度')
ylabel('从动件角速度dth(2),dth(3)(rad/s)')
text(250,200,'连杆2角速度')
text(50,100,'摇杆3角速度')
5. r1旋转360°时，α2，α3程序M文件：
r(1)=50;
r(2)=150;
r(3)=300;
r(4)=400;
th(1)=0;
dr=pi/180;
dth(2)=58.4119*dr;
dth(3)=154.7912*dr;
dth(1)=36;
dth2=1*dr;
for k=1:360
x=rrrvel([th(1),dth(2),dth(3),dth(1),r(1),r(2),r(3)]);
y=jweiyi([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;w 2]+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;
end
figure(3)
plot(jiao/dr,a2,jiao/dr,a3)
axis([0,360,-5000,5000])
grid on
title('角加速度线图')
xlabel('曲柄旋转转角')
ylabel('从动件角加速度（rad/s^2）')
text(200,40,'连杆2的角加速度')
text(150,-40,'摇杆3的角加速度')。