1、%x（4）=theta-4 guess value%x（5）=l1%x（6）=l2%x（7）=l3%x（8）=l4%x（9）=l6%x（10）=lA guess value%x（11）=lB%Output parameters%y（1）=lA%y（2）=theta-2%y（3）=theta-3%y（4）=theta-4theta2=x（2）;theta3=x（3）;theta4=x（4）;lA=x（10）epsilon=1.0E-6;f=x（6）*cos（theta2）-x（7）*cos（theta3）-x（8）*cos（pi+theta4）+x（5）. *cos（x（1）+pi）-x（9）;
2、 x（6）*sin（theta2）-x（7）*sin（theta3）-x（8）*sin（theta4+pi）+. x（5）*sin（x（1）+pi）; -x（11）*cos（theta4）+lA*cos（x（1）-x（9）; -x（11）*sin（theta4）+lA*sin（x（1）;while norm（f）epsilon J=0 -x（6）*sin（theta2） x（7）*sin（theta3） -x（8）*sin（theta4）; 0 x（6）*cos（theta2） -x（7）*cos（theta3） x（8）*cos（theta4）; cos（x（1） 0 0 x（11）*sin
3、（theta4）; sin（x（1） 0 0 -x（11）*cos（theta4）; dth=inv（J）*（-1.0*f）; lA=lA+dth（1）; theta2=theta2+dth（2）; theta3=theta3+dth（3）; theta4=theta4+dth（4）; f=x（6）*cos（theta2）-x（7）*cos（theta3）-x（8）*cos（pi+theta4）+x（5）.norm（f）;end;y（1）=lA;y（2）=theta2;y（3）=theta3;y（4）=theta4;再进行数据输入，运行程序进行运算。这里我们根据上面分析的1 的极限位置取1 的
4、范围为4055并均分成15个元素：clcclearx1=linspace（40*pi/180,55*pi/180,15）;x=zeros（length（x1）,11）;for n=1:15x（n,:）=x1（:,n） pi/6 8*pi/9 2*pi/3 40 50 75 35 70 75 60;endp=zeros（length（x1）,4）;for k=1:y= rrrposi（x（k,:）;p（k,:）=y; pp = 93.3149 0.7163 2.5455 1.5461 91.3071 0.7045 2.5617 1.5902 89.2387 0.6929 2.5786 1.634
5、7 87.1076 0.6815 2.5963 1.6796 84.9113 0.6703 2.6147 1.7250 82.6463 0.6592 2.6339 1.7709 80.3086 0.6482 2.6539 1.8174 77.8931 0.6372 2.6747 1.8646 75.3930 0.6263 2.6965 1.9126 72.7998 0.6154 2.7192 1.9616 70.1019 0.6043 2.7431 2.0118 67.2833 0.5930 2.7683 2.0635 64.3217 0.5812 2.7950 2.1169 61.1835
6、0.5687 2.8237 2.1728 57.8153 0.5551 2.8549 2.2319输出的P、矩阵的第二列到第四列分别是2 、3 、4 的值，第一列是AF杆的长度L1。第二步进行速度计算：根据速度方程式编写如下rrrvel.m函数：function y=rrrvel（x）%x（2）=theta-2%x（3）=theta-3%x（4）=theta-4%x（5）=dtheta-1%x（6）=l1%x（7）=l2%x（8）=l3%x（9）=l4%x（10）=l6%x（11）=lA%x（12）=lB%Outout parameters%y（1）=V%y（2）=dtheta-2%y（3）=
7、dtheta-3%y（4）=dtheta-4A=-x（7）*sin（x（2） x（8）*sin（x（3） x（9）*sin（pi+x（4） 0; x（7）*cos（x（2） -x（8）*cos（x（3） -x（9）*cos（x（4）+pi） 0; 0 0 x（12）*sin（x（4） cos（x（1）; 0 0 -x（12）*cos（x（4） sin（x（1）;B=x（6）*sin（x（1）+pi）;-x（6）*cos（x（1）+pi）;x（11）*sin（x（1）;-x（11）*cos（x（1）*x（5）;y=inv（A）*B;根据第一步得到的数据进行数据输入，运行程序计算各速度值。程序和所
8、得结果如下：x2=x1 p（:,2） p（:,3） p（:,4） 10*ones（15,1） 40*ones（15,1） 50*ones（15,1）. 75*ones（15,1） 35*ones（15,1） 70*ones（15,1） p（:,1） 60*ones（15,1）;q=zeros（4,15）;for m=1:y2=rrrvel（x2（m,:q（:,m）=y2;qq = 1.0e+003 * Columns 1 through 8 -0.0064 -0.0062 -0.0061 -0.0061 -0.0060 -0.0059 -0.0059 -0.0058 0.0085 0.0089
9、 0.0092 0.0096 0.0101 0.0105 0.0109 0.0114 0.0235 0.0237 0.0239 0.0241 0.0244 0.0247 0.0250 0.0255 -1.0578 -1.0897 -1.1226 -1.1568 -1.1926 -1.2302 -1.2704 -1.3137 Columns 9 through 15 -0.0058 -0.0059 -0.0060 -0.0062 -0.0065 -0.0069 -0.0078 0.0119 0.0125 0.0131 0.0139 0.0148 0.0159 0.0175 0.0259 0.02
10、65 0.0272 0.0281 0.0292 0.0306 0.0327 -1.3610 -1.4136 -1.4734 -1.5431 -1.6273 -1.7337 -1.8767 程序运行得到q矩阵，第一行到第三行分别是a2、a3、a4 的值，第四行是杆AF上滑块运动的速度，即F点的速度。第三步进行加速度计算：编写加速度计算函数rrra.m：function y=rrra（x）%x（1）=th1%x（2）=th2%x（3）=th3%x（4）=th4%x（5）=dth1%x（6）=dth2%x（7）=dth3%x（8）=dth4%x（9）=r1%x（11）=r3%x（12）=r4%x（1
11、3）=r6%x（14）=lA%x（15）=lB%x（16）=V%y（1）=ddth2%y（2）=ddth3%y（3）=ddth4%y（4）=aA=-x（10）*sin（x（2） x（11）*sin（x（3） -x（14）*sin（x（4） 0; x（10）*cos（x（2） -x（11）*cos（x（3） x（14）*cos（x（4） 0; 0 0 x（15）*sin（x（4） cos（x（1）; 0 0 -x（15）*cos（x（4） sin（x（1）;B=-x（10）*x（6）*cos（x（2） x（11）*x（7）*cos（x（3） -x（12）*x（8）*cos（x（4） 0; -x（
12、10）*x（6）*sin（x（2） -x（11）*x（7）*sin（x（3） -x（12）*x（8）*sin（x（4） 0; 0 0 x（15）*x（8）*cos（x（4） -x（5）*sin（x（1）; 0 0 x（15）*x（8）*sin（x（4） x（5）*cos（x（1）;C=x（6）;x（7）;x（8）;x（16）;D=x（9）*x（5）*cos（x（1）;x（9）*x（5）*sin（x（1）;x（14）*x（5）*cos（x（1）+x（16）*sin（x（1）;x（14）*x（5）*sin（x（1）-x（16）*cos（x（1）;y=inv（A）*D-inv（A）*B*C;根据第一
13、步和第二步输入数据，运行程序得到各加速度的值：x3=x1,4） 10*ones（15,1） q（1,:） q（2,: q（3,: .40*ones（15,1） 50*ones（15,1） 75*ones（15,1） 35*ones（15,1） 70*ones（15,1） p（:,1） .60*ones（15,1） q（4,:;f=zeros（4,15）;y3=rrra（x3（m,:f（:,m）=y3; ff =1.0e+005 *Columns 1 through 8 -0.0038 -0.0040 -0.0042 -0.0044 -0.0047 -0.0050 -0.0054 -0.0058
14、0.0060 0.0064 0.0068 0.0072 0.0077 0.0083 0.0089 0.00970.0033 0.0036 0.0039 0.0042 0.0046 0.0050 0.0056 0.0062 -0.3601 -0.3723 -0.3866 -0.4033 -0.4231 -0.4468 -0.4755 -0.5107Columns 9 through 15 -0.0064 -0.0071 -0.0081 -0.0094 -0.0112 -0.0139 -0.0185 0.0106 0.0117 0.0131 0.0150 0.0176 0.0214 0.0276
15、0.0070 0.0080 0.0093 0.0110 0.0134 0.0171 0.0231 -0.5548 -0.6111 -0.6851 -0.7857 -0.9290 -1.1462 -1.5071接下来计算C点在1 =55，w1 =10rad/s时的速度，加速度：Vx=40*10*sin（55）-50*q（1,10）*sin（p（10,2）;Vy=40*10*cos（55）+50*q（1,10）*cos（p（10,2）;ax=100*40*cos（50）-q（1,10）2*50*cos（p（10,2）-f（1,10）*50*sin（p（10,2）;ay=100*40*cos（50）
16、-q（1,10）2*50*cos（p（10,2）-f（1,10）*50*sin（p（10,2）;输出结果：Vx= -230.2208;Vy=-231.1533;ax= 2.3069e+004;ay= 2.3069e+004;表1 各构件的位置、速度和加速度1234L1W2W3W4Va2alradmm rad/smm/sx103 rad/s2 x103mm/s2x1050.6981 0.7168 0.7355 0.7542 0.7729 0.7916 0.8103 0.8290 0.8477 0.8664 0.8851 0.9038 0.9225 0.9412 0.95990.7163 0.70
17、45 0.6929 0.6815 0.6703 0.6592 0.6482 0.6372 0.6263 0.6154 0.6043 0.5930 0.5812 0.5687 0.55512.5455 2.5617 2.5786 2.5963 2.6147 2.6339 2.6539 2.6747 2.6965 2.7192 2.7431 2.7683 2.7950 2.8237 2.85491.5461 1.5902 1.6347 1.6796 1.7250 1.7709 1.8174 1.8646 1.9126 1.9616 2.0118 2.0635 2.1169 2.1728 2.231
18、993.3149 91.3071 89.2387 87.1076 84.9113 82.6463 80.3086 77.8931 75.3930 72.7998 70.1019 67.2833 64.3217 61.1835 57.8153-6.3578 -6.2487 -6.1469 -6.0541 -5.9726 -5.9054 -5.8561 -5.8299 -5.8340 -5.8786 -5.9789 -6.1587 -6.4572 -6.9440 -7.75558.4725 8.8575 9.2469 9.6433 10.0502 10.4720 10.9145 11.3856 1
19、1.8957 12.4601 13.1005 13.8507 14.7646 15.9364 17.547023.5099 23.6948 23.9018 24.1350 24.3994 24.7013 25.0491 25.4538 25.9305 26.5003 27.1937 28.0569 29.1637 30.6408 32.7288-1.0578 -1.0897 -1.1226 -1.1568 -1.1926 -1.2302 -1.2704 -1.3137 -1.3610 -1.4136 -1.4734 -1.5431 -1.6273 -1.7337 -1.8767-0.3790
20、-0.3960 -0.4156 -0.4386 -0.4656 -0.4976 -0.5361 -0.5829 -0.6409 -0.7144 -0.8099 -0.9386 -1.1206 -1.3948 -1.84910.6037 0.6395 0.6788 0.7225 0.7715 0.8274 0.8919 0.9677 1.0587 1.1707 1.3128 1.4999 1.7586 2.1406 2.76000.3334 0.3590 0.3879 0.4210 0.4592 0.5038 0.5567 0.6205 0.6989 0.7976 0.9255 1.0976 1
21、.3402 1.7052 2.3084-0.3601 -0.3723 -0.3866 -0.4033 -0.4231 -0.4468 -0.4755 -0.5107 -0.5548 -0.6111 -0.6851 -0.7857 -0.9290 -1.1462 -1.5071接下来输出图像：角位置图像，程序如下plot（x1,p（:,2）,-,x1,p（:,3）,:,4）,*）title（角位置）;xlabel（theta1/radylabel（theta2、theta3、theta4/rad输出图像如图：图4 2 、3 、4 的变化图像AF长度图像，程序如下：,1）L1/mm输出图像如下：图
22、5 AF长度变化图像角速度w2、w3、w4：plot（x1,q（1,:）,x1,q（2,:,x1,q（3,:W1 rad/sW2、W3、W4 rad/S角速度图6 角速度w2、w3、w4变化图像F点速度：plot（x1,q（4,:）V mm/sF点速度图7 F点速度变化图像角加速度：plot（x1,f（1,:,x1,f（2,:,x1,f（3,:alpha2、alpha3、alpha4 rad/S2角加速度图8 角加速度2、3、4变化图像F点的加速度：plot（x1,f（4,:a mm/s2F点的加速度图9 F点加速度变化图像参考文献【1】 孙恒，陈作模.机械原理【M】.7版.北京：高等教育出版社，2006【2】 曲秀全.基于MATLAB/Simulink平面连杆机构的动态仿真.哈尔滨工业大学出版社，2007【3】 云周工作室.MATLAB数学建模基础教程.人民邮电出版社.2001.7【4】 郭仁生.机械工程设计分析和MATLAB应用.机械工业出版社.2008.6.第二版【5】 同济大学数学系.线性代数.北京：高等教育出版社，2009【6】 同济大学数学系.高等数学. 北京：高等教育出版社，2008【7】 雷培，刘云霞.基于MATLAB的四连杆机构运动分析.机械工程与自动化.2009.4.第2期
